© The Authors 2022

1. Introduction

Galaxies are fundamentally shaped by their star formation history. At any given epoch, the star formation rate (SFR) determines the population of massive, short-lived stars that regulate and change the composition and energetics of the interstellar medium (ISM). The SFR sets the rate at which the ISM is dispersed through galactic winds or locked away into long-lived, low-mass stars. It correlates with the nuclear activity of galaxies and the growth of stellar mass throughout cosmic time (Madau et al. 1996; Bouwens et al. 2011). Finally, the SFR is a fundamental input and observational constraint for galaxy evolution models (Chiappini et al. 1997; Kauffmann & Haehnelt 2000; de Rossi et al. 2009), again because it sets both the level of massive star feedback and the rate at which gas and dust are depleted.

In external unresolved galaxies the SFR is most often determined via a wide range of proxy measurements (see e.g., Kennicutt & Evans 2012). In the Milky Way, one can – at least in principle – count up the young stars. Low-mass stars can be discerned as being young if they have not yet reached the main sequence; massive stars can be discerned as young by their position in the colour-magnitude diagram (CMD), or simply by the fact that they are short-lived. However, not all stars can readily be recognised as young. Therefore, any such counting exercise must account for the fraction of the stellar mass function enclosed in such a census.

For nearby (d ≲ 500 pc) molecular clouds in the Milky Way, nearly complete lists of young stellar objects (YSOs) are available, the ages of which (t) can be derived by measuring their CMD positions or infrared excesses. The SFR is therefore given by ⟨Ṁ_∗⟩ = N_YSO ⟨M_YSO⟩/t, where ⟨M_*⟩ is the mean mass of YSOs and is usually assumed to be ⟨M_*⟩ = 0.5 M_⊙ (see e.g., Lada et al. 2010). In the local group, analogous SFR estimates have been based on the CMD position of massive luminous stars. For example, Williams et al. (2011) have applied CMD fitting techniques to reconstruct spatially resolved maps of the stellar age distribution in nearby galaxies.

The ‘current’ SFR of the Milky Way has been estimated by various groups over the last decade, for instance by Robitaille & Whitney (2010), Chomiuk & Povich (2011), and Licquia & Newman (2015). Specifically, Robitaille & Whitney (2010) used the census of YSOs from Robitaille et al. (2008, hereafter R08) to construct a population synthesis model for YSOs in the Galaxy, applied the same observational constraints as for the R08 census, and varied the model SFR such that the number of detected synthetic YSOs matched the observed number. They inferred a global SFR (across the entire Milky Way) of 0.68–1.45 M_⊙ yr⁻¹. Chomiuk & Povich (2011) normalised various estimates for the Galactic SFR to the same initial mass function (IMF) and population synthesis models across methods, and estimated ⟨Ṁ_∗⟩ = 1.9 ± 0.4 M_⊙ yr⁻¹. Licquia & Newman (2015) revisited this work and estimated ⟨Ṁ_∗⟩ = 1.65 ± 0.19 M_⊙ yr⁻¹ by combining previous measurements from the literature using a hierarchical Bayesian statistical method.

A conceptual limitation of such analyses is the identification of a measured age distribution of stars in a Galactic neighbourhood with the local or global star formation history of the Milky Way, for two reasons. First, young stars that are now seen in a modest volume (say, a few 100 pc) co-rotating with the Sun were not necessarily born within that volume. Second, we know from external disc galaxies that the instantaneous SFR varies dramatically as a function of position within a galaxy (see e.g., Fig. 1 in Kreckel et al. 2018, for a recent illustration); therefore, the extrapolation of a local census of young stars to a galaxy-wide SFR may be far off. To acknowledge this crucial distinction, but not completely discard traditional terminology, we use the term effective star formation rate for the quantity that is naively inferred from the age distribution. In addition, the SFR of galaxies may vary on short time scales, so any estimate of Ṁ_∗ ∼ M_∗(< Δt)/Δt will depend on the age range Δt of the young stars considered.

In this study, we set out to model and determine the age distribution over a substantive fraction (6 × 6 kpc around the Sun) of the Galactic disc, using upper main sequence, i.e. OB(A), stars. To do this, we devised and applied a non-parametric model¹ for the age distribution of stars, which predicts their absolute magnitude distribution via isochrones. These predictions are then compared to OB(A) stars across many spatial patches in the Galactic disc, resulting in age-resolved spatial maps of stars younger than 1 Gyr. On that basis, we explored to what extent the modelled age distribution and its effective SFR reflect the actual disc star formation history of the Milky Way. Recently, Ruiz-Lara et al. (2020) made a comprehensive estimate of the age distribution of Galactic disc stars within 2 kpc by modelling Gaia DR2 colour–magnitude diagrams. However, their focus was on 0.5–10 Gyr timescales as opposed to the younger populations that we study here. Our analysis draws in the sample of massive stars from Zari et al. (2021), which we describe in Sect. 2. We lay out our methodology in Sect. 3. Our results, presented in Sect. 4, show a rapid increase in star formation rate compared to previous years, which could be related to the presence of numerous massive star-forming regions in the area considered. We discuss our results and summarise them in Sect. 5.

2. Data: Luminous young disc stars

The sample that we used in this study consists of ∼500 000 OB(A) stars in the ‘filtered’ sample presented in Zari et al. (2021, hereafter Z21). The sample was selected by combining Gaia EDR3 photometry and astrometry and 2MASS photometry by applying the following criteria. First, Z21 selected sources brighter than G = 16 mag and with absolute magnitude in the 2MASS K_s band M_K < 0 mag, which roughly corresponds to a late B-type star. Then, they applied several cuts (described in Sect. 2.1 of Z21, Eqs. (2)–(5)) in colour–colour space to exclude bright red giant-branch and asymptotic giant-branch stars. To clean the sample from spurious astrometric solutions, they removed sources with astrometric fidelity < 0.9 (Rybizki et al. 2022). Finally, they selected sources in the Galactic disc that are and remain near the disc midplane, i.e. with |z|< 300 pc, and that belong to a low-velocity dispersion component in the vertical direction (identified with a Gaussian mixture model; see Sect. 4.2 of Zari et al. (2021). The latter condition was used to exclude stars with large vertical motions that point to a sample contamination by old, hot stars.

The Z21 sample extends to distances of ∼5 kpc from the Sun (Fig. 1, left). However, in this study we focused on a 6 kpc × 6 kpc region of the Galactic disc (centred on the Sun, shown in Fig. 1, right), where we can reasonably assume that our sample, is complete. To derive our spatially resolved SFR map of the Galactic disc, we divided our volume into 256 spaxels of different sizes, each containing the same number of stars (∼1000). We modelled the magnitude distribution of the sources in each spaxel, as described in the next section.

Fig. 1. Density distribution of the filtered sample of massive hot stars from Zari et al. (2021). In the left panel, the coloured lines represent the spiral arm locations derived by Reid et al. (2019). The right panel highlights the 6 × 6 kpc region considered in this study. Such region is divided in ‘spaxels’ of different sizes, each containing the same number of sources. In both panels, the Sun is in X, Y = (0, 0).

3. Modelling the stellar age distribution

Ideally, we would like to model the observed colour-magnitude distribution of our sample stars in various spaxels to constrain their age distribution. However, this would require detailed star-by-star de-reddening, which would be difficult if not impossible (however, see e.g., Ruiz-Lara et al. 2020). Therefore, we resorted to modelling a more restricted set of observables centred on the near-infrared K band, which is less affected by extinction; we modelled the distribution of absolute magnitudes M_K for stars that are luminous (M_K < 0) and (intrinsically) bluer than the red-giant branch. This latter colour criterion can be satisfied even in the presence of extensive reddening. These stars should all be massive young main-sequence or evolved stars (Z21).

Specifically, we chose the absolute K-band magnitude distribution of sample members as a statistic n(M_K | x,α), which we take to be a function of position in the Galaxy x. We assume that at each x there is an underlying age or birthrate distribution, b(τ | x, α), whose temporal dependence is described by a set of parameters α, specified in Eq. (2). The expected observed distribution in absolute magnitudes then depends on stellar evolution, the initial mass function, the birthrate distribution, and the observational selection function:

$$\begin{array}{cc}& n\left(M_K|\mathit{x},\mathit{\alpha }\right)={\displaystyle \int }\mathrm{d}\tau {\displaystyle \int }\mathrm{d}{m}_0{\displaystyle \int }\mathrm{d}{T}_{\mathrm{eff}}\hfill \\ \hfill & \times S_c(M_K,T_{\mathrm{eff}})p(M_K,T_{\mathrm{eff}}|\tau ,m_0)\xi (m_0)b(\tau |\mathit{x},\mathit{\alpha }).\hfill \end{array}$$(1)

In Eq. (1), the term S_c(M_K, T_eff) reflects the observational selection function of the sample as specified in Eq. (3); p(M_K, T_eff|τ, m₀) reflects the probability density of a star to have absolute magnitude M_K and effective temperature T_eff at a given birth mass m₀ and age τ; and ξ(m₀) is the initial mass function in units of mass [M_⊙]⁻¹. The integral over age τ should in principle extend from 0 to t_Hubble. In practice, we consider minimum and maximum ages from τ_min to τ_max, as the very youngest stars are too enshrouded to be in the sample, and stars with τ ≳ 1 Gyr will be too faint in their main-sequence phase. In the next section, we describe the single terms of Eq. (1).

3.1. Functional forms of the model components

We specify the birthrate b(τ | x, α) as a function of time τ at a given position in the Galaxy x, given parameters α. We want to avoid imposing a particularly restrictive functional form, as we have no reason to believe that the age distribution is particularly smooth or spatially homogeneous. Therefore, we parameterise it by step functions in time that can be different in each spatial pixel x:

$$\begin{array}{c}\hfill b(\tau |\mathit{x},\mathit{\alpha })={\displaystyle \sum _{\text{age bin}i}}{\alpha }_i(\mathit{x}){\chi }_{A_i}(\tau ),\end{array}$$(2)

where A_i are age intervals, the rate coefficients α_i are the elements of α and must satisfy α_i > 0, and the indicator function χ_Ai = 1 if τ ∈ A_i and zero otherwise. After some experimentation to explore the trade-off between functional flexibility and the noisiness of individual bins, we adopted five bins:

$$\begin{array}{cc}& A_1=5-10\mathrm{Myr}\hfill \\ \hfill & A_2=10-30\mathrm{Myr}\hfill \\ \hfill & A_3=30-100\mathrm{Myr}\hfill \\ \hfill & A_4=100-300\mathrm{Myr}\hfill \\ \hfill & A_5=300\mathrm{Myr}-1\mathrm{Gyr},\hfill \end{array}$$

spaced approximately logarithmic in age, which corresponds to a linear dependence in turn-off magnitude.

We use a broken power law Kroupa IMF (Kroupa 2001, 2002) whose piecewise components ξ(m)≈m^−Γ are

$$\begin{array}{cc}& \mathrm{\Gamma }=0.3\mathrm{for}0.01<M/M_{\odot }<0.08\hfill \\ \hfill & \mathrm{\Gamma }=1.8\mathrm{for}0.08<M/M_{\odot }<0.5\hfill \\ \hfill & \mathrm{\Gamma }=2.7\mathrm{for}0.5<M/M_{\odot }<1\hfill \\ \hfill & \mathrm{\Gamma }=2.3\mathrm{for}1<M/M_{\odot }<120.\hfill \end{array}$$

This form of the IMF accounts for unresolved binaries, which are important across most of the IMF’s mass spectrum (e.g., Moe & Di Stefano 2017). Previous studies (Robitaille & Whitney 2010; Chomiuk & Povich 2011) have not used IMF forms that make this correction. The choice of the IMF matters in this context, as all such studies match the observed frequency of young stars in a finite mass range and then extrapolate to the SFR across the full range of ξ(m). In Sect. 5, we discuss how our results would have been affected using an alternate, simpler IMF model (e.g., Chomiuk & Povich 2011; Robitaille & Whitney 2010).

The term S_c in Eq. (1) is the selection function of the sample, i.e. a function that returns the (unitless) probability that a star is in the photometric catalogue we model, given its observed properties. In this work, we approximate the selection function described in Sect. 2 and Z21 as a function of M_K and T_eff. We presume that S_c(M_K, T_eff) is either 1 or 0. Since M_K is an observable, the function S_c in M_K is directly linked to the selection of the observed sample. The dependence of S_c on T_eff is somewhat more indirect, as it is not a direct photometric observable, especially in the presence of reddening. In practise, we only need to know T_eff well enough to eliminate RGB stars at that luminosity, which is possible even for reddened stars (Z21). The selection function S_c assures that the integral in Eq. (1) has non-zero contributions only from parts of the isochrones where both M_K and T_eff are within the ‘selected’ range. We approximate the sample selection function as

$$\begin{array}{c}\hfill S_c(M_K,T_{\mathrm{eff}})=\{\begin{array}{c}1M_K<min\{0,M_K^{\max }(\mathit{x})\}\wedge T_{\mathrm{eff}}>7000\mathrm{K}\hfill \\ 0\text{otherwise.}\hfill \end{array}\end{array}$$(3)

The condition $M_K<\min \{0,M_K^{\text{max}}(\mathit{x})\}$ arises from the fact that for the most distant patches, x, the apparent magnitudes, G, of stars with absolute magnitudes of M_K = 0 may be too faint and lie in a magnitude range where Gaia is not necessarily complete. We therefore impose a brighter cut to ensure that we can fully model the selection function. We derived $M_K^{\text{max}}(\mathit{x})$ by smoothing the observed M_K distribution with a Gaussian kernel and taking the mode of the distribution. This reduces the number of stars in each spaxel to 550–800 stars, depending on distance (the more distant spaxels have fewer sources).

Concerning the probability density of M_K, the overall data-model comparison requires us to predict how many hot stars of luminosity M_K we would expect from isochrones for a stellar population of unit mass and age τ: p(M_K|τ) (in units of mag⁻¹). To compute p(M_K|τ), we first applied our selection function S_c(M_K, T_eff) and marginalised (integrate) over all T_eff and m₀. We used PARSEC isochrones (Bressan et al. 2012; Tang et al. 2014; Chen et al. 2015) from 1 Myr to 1 Gyr, equally spaced by 0.1 dex in log(age/yr). Then, for each age, we sampled the Kroupa IMF and we assigned to each star the appropriate absolute magnitude M_K by finding the nearest isochrone point according to the mass of the star. Finally, we counted the number of stars in bins with widths of ΔM_K = 0.25 mag, between M_K, min and M_K, max, and we normalised such number for the total number of stars and the bin width ΔM_K. Therefore, we obtain the quantity n(M_K, τ), which corresponds to the integral over the initial mass and effective temperature m₀, T_eff of Eq. (1):

$$\begin{array}{cc}\hfill n(M_K,\tau )& ={\displaystyle \int }\mathrm{d}{m}_0{\displaystyle \int }\mathrm{d}{T}_{\mathrm{eff}}\hfill \\ \hfill & \times S_c(M_K,T_{\mathrm{eff}})p(M_K,T_{\mathrm{eff}}|\tau ,m_0)\xi (m_0).\hfill \end{array}$$(4)

Since we do not know the individual ages of the stars in bin x, we also needed to integrate over isochrone age τ before comparing with the data:

$$\begin{array}{c}\hfill n(M_K|\mathit{x},\mathit{\alpha })={\displaystyle {\int }_{{\tau }_{\min }}^{{\tau }_{\max }}}n(M_K,\tau )b(\tau |\mathit{x},\mathit{\alpha })\mathrm{d}\tau ,\end{array}$$(5)

which we did numerically. As the isochrones are logarithmically spaced, it proves useful to define the logarithmic age ω = log τ and transforming integration variables to obtain

$$\begin{array}{c}\hfill n(M_K|\mathit{x},\mathit{\alpha })={\displaystyle {\int }_{{\omega }_{\min }}^{{\omega }_{\max }}}n(M_K|{10}^{\omega })b({10}^{\omega })ln(10)e^{ln(10)\omega }\mathrm{d}\omega ,\end{array}$$(6)

where we dropped the arguments x and α from b for conciseness. Now, Eq. (6) represents the predicted M_K distribution at a given position, x, given the model assumptions about b(τ | x, α).

3.2. The data likelihood function

We compared the predictions from Eq. (6) to the observed magnitude distribution {M_K, i} detected in patch x by quantifying and optimising the likelihood of the data. We assume that the stars in our dataset are sampled by an inhomogeneous Poisson point process with a rate function of λ = n(M_K |x, α). This leads the joint probability of the data, given the model parameters α of p(N = N_⋆, {M_K}|α) (the likelihood) to be the product of the two terms specified below. The first term,

$$\begin{array}{c}\hfill \frac{{\left[{\displaystyle \int }n(M_K|\mathit{x},\mathit{\alpha })\mathrm{d}{M}_K\right]}^N}{N!}exp[-\int n(M_K|\mathit{x},\mathit{\alpha })\mathrm{d}{M}_K],\end{array}$$(7)

states that the probability of observing N stars in a spaxel follows a Poissonian distribution and thus sets the normalisation and informs us of the total size of the sample. The second term,

$$\begin{array}{c}\hfill {\displaystyle \prod _i}\frac{n(M_{K,i}|\mathit{x},\mathit{\alpha })}{\int n(M_K|\mathit{x},\mathit{\alpha })\mathrm{d}{M}_K},\end{array}$$(8)

specifies the form of the distribution of the data and states that the observed rate is normalised over the volume in M_K space that could have been observed within the survey selection constraints as expressed by the survey selection function (e.g., Bovy et al. 2012; Rix & Bovy 2013; Rix et al. 2021). By multiplying Eqs. (7) and (8), we obtain

$$\begin{array}{c}\hfill \mathcal{L}(\{M_{K,i}\}|\mathit{x},\mathit{\alpha })=({\displaystyle \prod _{i=1}^{N_{\ast }}}n(M_{K,i}|\mathit{x},\mathit{\alpha }))e^{-N_{\ast ,\mathrm{tot}}},\end{array}$$(9)

where the normalisation factor is

$$\begin{array}{c}\hfill N_{\ast }^{\mathrm{tot}}(\mathit{x})={\displaystyle {\int }_{M_{K,\min }}^{M_{K,\max }=0}}n(M_K|\mathit{x},\mathit{\alpha })\mathrm{d}{M}_K.\end{array}$$(10)

Using uniform priors between 0 and 0.1, we sampled the α_i(x) parameters characterising the birthrate b(τ | x, α) in each spaxel from the posterior with emcee (Foreman-Mackey et al. 2013) and using, respectively, the 50th, 16th, and 84th percentiles as the parameter best estimate and the 1σ levels. We call b(τ | x, α_best) the birthrate obtained by using the best estimates for α_i.

3.3. Tests on mock data

To test our model, we generated a mock magnitude distribution n(M_K|x, α) (see Eq. (6)) using a known birthrate b(τ | x, α_true) = ∑_iα_i, trueχ_Ai(τ) (where we chose α_true = α_best), and the same selection function, set of isochrones, and IMF as used in our model. We then used our fitting procedure described in the previous section to derive the parameters α_i, pred and compare them with the true values α_i, true. To use this test in a physically meaningful range of the parameter space and to anticipate the results (Sect. 4), the α_i, true were chosen to reproduce the observed number of stars in a spaxel at a given position, x.

Figure 2 (top row) shows the true (black lines) and predicted values (orange lines) for α_i, for three realisations of the same mock magnitude distribution. The shaded areas correspond to 1σ and 2σ levels. Figure 2 (bottom row) shows the corresponding true mock magnitude distribution (black histogram) and the predicted magnitude distribution for the best fit α_i value (thick solid orange histogram) and the 1σ values (dashed thin histograms). The predicted parameters are compatible within 1−2σ with the true values, and they are always consistent with each other. The increase in birthrate in the second age interval (10–30 Myr) is always correctly retrieved. The difference between the estimated values of the parameters most likely depends on the different sampling of the bright end of the magnitude distribution (M_K < −4 mag), which is also where the true magnitude distribution differs the most from that estimated using the predicted parameters. This is due to the very low number of stars predicted (and observed) at extremely high luminosity.

Fig. 2. Model fits to mock datasets used as tests. Top row: black lines show the true value of the birthrate αi parameters. The orange lines show the estimated αi parameters for three different realisations of the mock datasets generated from the same initial parameters. The orange shaded areas correspond to the 1- and 2-σ levels. Bottom row: black histograms represent the magnitude distribution of three different mock datasets. The thick orange histograms represent the magnitude distribution obtained by using the best fit αi. The dashed histograms represent the magnitude distribution obtained by using the 16th and 84th percentiles for αi.

4. Results

In this section, we present the results obtained by fitting our model to the data presented in Sect. 2. We computed the values and uncertainties for the coefficients α_i, j that imply b(τ_i | x_j, α_j) for all the spaxels x_j in the map of Fig. 1 (right) by sampling the likelihood of Eq. (9). Example values of the α_i parameters and their uncertainties are given in Table 1 and are available on CDS.

These results are essentially number- or mass-density maps of mono-age stellar populations (in five age bins) across our surrounding 6 kpc  ×  6 kpc area of the Galactic disc (Figs. 3 and 4), as well as the overall radial gradients they show (Figs. 5 and 6). When summing over the entire area, we obtain an estimate of the integrated age distribution of the stars (Fig. 7), which reflects the formation rate of the stars that are now in this portion of the disc.

Fig. 3. Density distribution of stars younger than 500 Myr in the Galacic plane, as derived from Eq. (11), i.e. by using the best estimates for the birthrate coefficients, αi, in each spaxel, j.

Fig. 4. Density distribution of stars with ages of 5 < t < 50 Myr (top left), 50 < t < 100 Myr (top right), 100 < t < 200 Myr (bottom left), and 200 < t < 500 Myr (bottom right). The density distribution changes as a function of time, becoming gradually smoother. Over-arching structures corresponding to spiral arms are visible for all age intervals. The solid black lines indicate the location of the spiral arms from Reid et al. (2019).

Fig. 5. Distribution of star formation rate surface density as a function of galactocentric radius. The thick lines correspond to the average star formation rate surface density for different age intervals (red: 5 < τ < 10 Myr; orange: 50 < τ < 100 Myr; yellow: 400 < Myrτ < 1 Gyr). The shaded areas correspond to the 16th and 84th percentiles. The left Y-axis shows the star formation rate expressed as the total number of stars born; the right Y-axis shows the more conventional amount of stellar mass formed, which requires the adoption of an IMF (here, Kroupa 2001).

Fig. 6. Star formation rate density in different galactocentric radii, as a function of time. The left and right Y-axes show the SFR expressed as the total number of stars born and as the stellar mass formed, respectively. The thick lines correspond to the average star formation rate surface density for three intervals in galactocentric radius. The shades areas correspond to the 16th and 84th percentiles. Ages below 5 Myr are not included in our analysis, and thus the corresponding area in the plot is hatched.

Fig. 7. Effective star formation history of the sample, in units of yr−1 (left y-axis) and M⊙ yr−1 (first right y-axis). The second right y-axis corresponds to the total SFR obtained by multiplying by the ratio f = SFR5 − 11/SFRtot (see Sect. 5). Vertical bars correspond to errors on the parameters. Ages below 5 Myr are not included in our analysis, and thus the corresponding area in the plot is hatched.

4.1. Mono-age stellar density maps of the Galactic disc

Figure 3 shows the number surface density of stars, n(Δτ | x, α), younger than 500 Myr, defined as

$$\begin{array}{c}\hfill n(\mathrm{\Delta }\tau |\mathit{x},\mathit{\alpha })\equiv \frac{1}{\mathrm{Area}(\mathit{x})}{\displaystyle {\int }_{{\tau }_{\min ,i}}^{{\tau }_{\max ,i}}}{b}_{\mathrm{best}}(\tau |\mathit{x},{\mathit{\alpha }}_{\mathrm{best}})\mathrm{d}\tau ,\end{array}$$(11)

where the integral of b(τ | x, α_best) over time gives the number of stars formed between τ_min < τ < τ_max and Area(x) is the spaxel area in kpc². The density distribution in the Galactic plane of sources younger than 500 Myr is qualitatively very similar to that shown in Fig. 1, lending credence to our results. The total number of stars in each spaxel is different from the number of observed stars per spaxel, because with Eq. (11) we predict the total number of stars at all masses, not only those included in the filtered sample from Z21. These maps can be converted into maps of the (hypothetical, initial) stellar surface mass density of stars in that age bin, simply by multiplying by the initial mean mass per star, ⟨M_*⟩ = 0.22 M_⊙ for our Kroupa IMF.

The density distribution of the sample can be further investigated by splitting the stars into age bins, with ages of 5 < τ < 50 Myr (Fig. 4, top left), 50 < τ < 100 Myr (Fig. 4, top right), 100 < τ < 200 Myr (Fig. 4, bottom left), and 200 < τ < 500 Myr (Fig. 4, bottom right). The appearance of these maps changes as a function of age. In the youngest age interval (5 < τ < 50 Myr), the density distribution traces massive star-forming regions, which are visible as prominent overdense clumps in Figs. 1 and 3 (see also Zari et al. 2021; Poggio et al. 2021). The distributions of stars in the less young age intervals appear smoother and do not show dense clumps, certainly not beyond the age of one dynamical timescale or galactocentric rotation: ∼200 Myr. However, arch-like or spiral structures are visible for all age groups. The strength of spiral arms is quantified in Fig. A.1.

4.2. Radial trends of the young star density

For potential comparison with in particular external galaxies (González Delgado et al. 2016), we now look at the dependence of the young star density on galactocentric radius R_GC, averaging over the azimuthal angle that lies inside our 6 kpc  ×  6 kpc box, simply summing over all spaxels in a certain R_GC range:

$$\begin{array}{c}\hfill n_{\mathrm{R}}(\tau |R_{\mathrm{GC}})\equiv \frac{1}{S}{\displaystyle \sum _{\text{spaxels at}{R}_{\mathrm{GC}}}}b(\tau |\mathit{x},{\mathit{\alpha }}_{\mathrm{best}}),\end{array}$$(12)

where S is the total area of these spaxels.

Figure 5 shows the radial distribution of n_R(τ_i | R_GC) for three different age intervals τ_i: 5 < τ < 10 Myr, 50 < τ < 100 Myr, and 400 Myr < τ < 1 Gyr. While the older age bins only show a shallow gradient, the radial density gradient among stars τ < 10 Myr is quite steep, and shows a prominent peak at R ∼ 6 kpc and hints of two peaks at R ∼ 8 kpc and R ∼ 10 kpc.

Figure 6 shows the azimuthally averaged mean birthrate parameter n_R(τ_i|R) for different galactocentric radii. The different galactocentric annuli do not only vary significantly in their normalisation, but also their relative age distributions. While all show a distinct increase at the youngest ages, this increase is fractionally largest well outside the solar radius. The age distribution well inside the solar radius may have had a second peak between 100 and 300 Myr ago (which may be dominated by the contributions along the galactocentric angle ϕ = −0.115; see Fig. A.1).

4.3. Effective overall star formation rate

As the final step in the presentation of our results, we show the age distribution when summed over the entire 6 × 6 kpc² area. We define the total ‘effective’ birthrate b_tot as

$$\begin{array}{c}\hfill b_{\mathrm{tot}}(\tau )={\displaystyle \sum _{\text{all spaxels}\mathit{x}}}b(\tau |\mathit{x},{\mathit{\alpha }}_{\mathrm{best}})\end{array}$$(13)

by summing over all the j spaxels in the map (and not normalising by the area). Below we discuss to what extent this age distribution reflects the overall star formation rate across the Galactic disc.

Units of b_tot(τ) are of stars/year, which we converted to the more conventional M_⊙ yr⁻¹ by multiplying by the mean stellar mass derived from our assumed IMF ⟨M⟩ = 0.22 M_⊙. This is shown in Fig. 7, along with b_tot(τ) on the left Y-axis and SFR_eff on the first right Y-axis. Both values refer to the entire 6 × 6 kpc² area, a much larger area than previous quantitative studies for this age range, but of course not the ‘entire disc’.

Figure 7 shows that the age distribution of stars is approximately uniform between 20 Myr and 1 Gyr at 0.4 yr⁻¹ (or 0.08 M_⊙ year⁻¹. Remarkably, the effective star formation rate seems to have been three times higher in the last ∼10 − 20 Myr. It appears that the Galactic disc, or at least the large patch we consider here, has just lit up in young stars.

5. Discussion and summary

In the previous sections, we estimate the recent star formation history of the extended solar neighbourhood and found that it has recently increased by a factor of ∼3. We now estimate the global star formation rate of the Galactic disc and put this result into context by comparing with previous literature studies. We then use simulations to quantify to what extent global conclusions can be extrapolated from local measurements.

5.1. Comparison to existing work

Assuming that the fraction of the disc that we are probing is representative of the entire Milky Way (see Sect. 5.2), we can derive the ‘total’ SFR of the Milky Way. To do so, we considered the youngest bin, τ ≲ 10 Myr, as the best approximation to the ‘current SFR’, and we extended SFR_{eff, 5 − 10 Myr} to the disc. Since the surface density of star formation decreases strongly as a function of radius, we presume that the star formation rate surface density is described well by an exponential radial profile (Eq. (14)) within 4.5–15 kpc, with a scale length of L = 3.5 kpc (Zari et al. 2021, see Appendix B.2):

$$\begin{array}{c}\hfill {\mathrm{SFR}}_{\mathrm{tot}}=2\pi {\displaystyle {\int }_{R_{\min }}^{R_{\max }}}{S}_0{e}^{-r/L}\mathrm{d}r.\end{array}$$(14)

We chose R_min = 4.5 kpc and R_max = 15 kpc following Reed (2005) and Davies et al. (2011) and thus excluded the central region of the Galaxy. In the radial range covered by our data, ∼5 and ∼11 kpc from the Galactic centre, the total SFR is

$$\begin{array}{c}\hfill {\mathrm{SFR}}_{5-11\mathrm{kpc}}=\mathrm{\Delta }\varphi {\displaystyle {\int }_5^{11}}{S}_0{e}^{-r/L}\mathrm{d}r,\end{array}$$(15)

where Δϕ ≈ 0.7 rad corresponds to the range of azimuth angles covered. The ratio f = SFR_5 − 11/SFR_tot ≈ 8% reflects the fraction of the MW disc considered in this study. We thus find that the total SFR of the Milky Way implied by this model is $3.3_{-0.6}^{+0.7}{M}_{\odot }{\mathrm{yr}}^{-1}$ (see Fig. 7).

Recent analogous estimates of the star formation rate of the Milky Way are between 1.65  ±  0.19 M_⊙ yr⁻¹ (Licquia & Newman 2015) and 1.9 ± 0.4 M_⊙ yr⁻¹ (Chomiuk & Povich 2011), similar to but lower than the value we find. Reasons for this discrepancy may reflect differences in the assumed IMF, the spatial coverage, and the age of the tracers used; these are discussed below.

In their Sect. 3, Chomiuk & Povich (2011) reviewed the different methods used to estimate the Milky Way SFR. Such methods rely on Lyman continuum photon rates (e.g., Mezger 1978; Smith et al. 1978; McKee & Williams 1997; Murray & Rahman 2010), infrared emission and YSO counts (Robitaille & Whitney 2010; Davies et al. 2011), or massive star counts and supernova rates (Reed 2005; Diehl et al. 2006). In their study, Chomiuk & Povich (2011; and thus also Licquia & Newman 2015) normalised previous results to the same IMF, adopting the IMF from Kroupa & Weidner (2003, KW03) as a reference: ξ(m)∝m^−Γ with Γ = 1.3 for masses between 0.1 and 0.5 M_⊙ and Γ = 2.3 for masses between 0.5 and 100 M_⊙. Since our sample is predominately composed of massive stars of which most will be in unresolved binary (or multiple) systems (see e.g., Sana et al. 2014), it seems appropriate to adopt an IMF that is explicitly corrected for unresolved binaries (Kroupa 2001, 2002) reported in Sect. 3. The choice of the IMF impacts, of course, the relation between the observed and modelled luminosity function n(M_K|τ) and the implied ‘underlying’ total stellar mass (Eqs. (4) and (6)). In our formulation of the problem, this affects the inference of the birthrate parameter, α_i, and the mean stellar mass that is used to convert the SFR to units of M_⊙ yr⁻¹, which for the KW03 IMF is ⟨M⟩ = 0.63 M_⊙. To quantify the impact of the choice of IMF on our results, we re-ran the analysis using KW03 IMF. This yields an analogous estimate for the total SFR of 2.3 ± 0.4 M_⊙ yr⁻¹, which is consistent with the estimates of Chomiuk & Povich (2011) and slightly higher than those of Licquia & Newman (2015).

There are strong assumptions behind the space distribution of different tracers used to study the global SFR or the Milky Way (including ours, see Sect. 5.2). For instance, the distance to HII regions are biased against distant sources, and studies that employ HII emission to compute the Galactic SFR usually account for this by doubling the luminosity of sources in our half of the Galaxy (Mezger 1978; Smith et al. 1978; McKee & Williams 1997; Murray & Rahman 2010). In other cases, distances are not available; for example, the space distribution of YSOs assumed by Robitaille & Whitney (2010; see their Fig. 2) is quite different than what we used in this work, which is instead more similar to what was proposed by Reed (2005) and Davies et al. (2011).

Finally, the lifetimes of the tracers used in other studies are different. The estimate (0.68 − 1.45 M_⊙ yr⁻¹) of Robitaille & Whitney (2010) assumed a YSO lifetime of ≈2 Myr, which is younger than the minimum age considered here (5 Myr). Murray & Rahman (2010) used a lifetime for HII regions of 3.7 Myr following McKee & Williams (1997), which is slightly longer than the 3 Myr used by Mezger (1978), and again younger than our minimum age. Reed (2005) assigned lifetimes from ≈3.5 Myr to ≈20 Myr based on the luminosities of his sample of O- and B-type stars. We note that our estimate of the global effective star formation rate is already 3× lower for stars that are 20–50 Myr old.

The recent increase in effective SFR that we find is qualitatively similar to that found by Ruiz-Lara et al. (2020). Ruiz-Lara et al. (2020) modelled Gaia DR2 observed colour–magnitude diagrams to infer the star formation history within 2 kpc around the Sun and found three star formation enhancements that occurred ∼5.7, 1.9, and 1.0 Gyr ago and a hint of a fourth possible star formation burst spanning the last 70 Myr. Ruiz-Lara et al. (2020) did not express their results in M_⊙ yr⁻¹, complicating a direct comparison. They proposed that such enhancements are due to recurrent interactions between the Milky Way and the Sagittarius dwarf galaxy (Ibata et al. 1994; Laporte et al. 2018). We confirm the recent enhancement in star formation activity of the Galactic disc with a higher level of confidence and we constrain it to the last 10 Myr. Given that we find distinct temporal variations in the Galaxy’s effective star formation rate over the last 100 Myr (Fig. 7), matching the exact age distribution of the tracers would be crucial for a more quantitative comparison with Cignoni et al. (2006) and Ruiz-Lara et al. (2020).

Figure 4 shows that the location of the density enhancements broadly coincides between the youngest and older stellar age intervals. This suggests that the spiral pattern is associated not only with young star-forming regions but also with the underlying mass density (Eskridge et al. 2002; Binney & Tremaine 2008), at least for τ < 500 Myr. Consequently, Fig. 4 represents evidence that those spiral arms that show up in all our mono-age maps are indeed stellar-mass overdensities. Spiral arms that are seen in stellar populations of different ages are also observed in external galaxies. For instance, Meidt et al. (2021) compared the distribution of molecular gas (tracing star formation) and 3.6 μm emission (tracing old disc populations) across 67 star-forming galaxies with different morphological properties, often finding coincidence. However, not all spiral arms in the Milky Way disc, inferred from masers, for example, are reflected in stellar-mass overdensities in our maps among stars with τ ≲ 0.5 Gyr (see top left panel of Fig. 4).

5.2. How closely does the best-fit age distribution reflect the overall disc SFR?

The area sampled by our OBA sample is only ≈8% of the entire Galactic disc. Large-scale spatial variations in the star formation rate density or stellar ages, such as large-scale spiral arms, could bias the generalisation from this sample to the entire Galactic disc. To test how much variance could be expected, we examine the final snapshot of cosmological simulations of Milky Way-like galaxies. We used the TNG50 cosmological simulation (Pillepich et al. 2019; Nelson et al. 2019), whose numerical resolution reaches that of zoomed-in simulations (reaching a baryonic mass of 8.5 × 10⁴ M_⊙). IN TNG50, we selected Milky Way-like galaxies based on their stellar mass, discyness, and morphology. Specifically, we require the following:

• 10^10.5 ≤ M_⋆/M_⊙ ≤ 10^11.2

• over 40% of stars on near-circular orbits (with circularities ≥0.7, Genel et al. 2015)

• a central bar

• a non-zero star formation rate at z = 0,

leading to N = 78 galaxies. We centred and align a Cartesian coordinate system with the vertical axis aligned with the angular momentum of the disc. To mimic our observed sample, we placed a solar neighbourhood at (R, z, φ) = (8, 0, φ) and selected stars in a cylinder centred on the Sun with a 3 kpc radius. To investigate how the location φ of the Sun influences the results, we varied the azimuth of the Sun and measured the local star formation history (from the local age distribution), thereby estimating the amplitude biases in local datasets.

Sampling all galaxies and all azimuthal angles for the position of the hypothetical sun, we derived the distribution of

$$\begin{array}{c}\hfill {log}_{10}\left(\frac{{\dot{M}}_{\star }/\mathrm{d}S(\mathrm{patch})}{{\dot{M}}_{\star }/\mathrm{d}S(\mathrm{annulus})}\right),\end{array}$$

where the denominator is the number of stars per age interval in the patch, and the denominator is the same quantity in the entire annulus. We find that azimuthal variations of the recent ‘effective’ SFR in a patch of R = 3 kpc reflect the global average typically within factor of ∼1.5. These variations are shown in Fig. 8 as orange bands, reflecting the conceptual uncertainties of extrapolating the age distribution measured within our volume to a global star formation rate of the Galactic disc. This figure shows that these uncertainties are considerable. Nonetheless, the increase in the frequency of very young stars we find reflects a global recent increase in the Galaxy SFR. While it is unlikely that the Sun has just moved into a patch of higher star formation activity, it is not obvious that there are plausible physical mechanisms that would boost the Milky Way’s global star formation rate on timescales shorter than ∼20 Myr. Our comparison with simulations implies that an SFR increase in our observed volume likely reflects a global SFR increase; however, we need to be open to the possibility that our SFR estimates for the youngest age bin reflect a localised SFR fluctuation with much more modest global SFR changes across the whole Milky Way.

Fig. 8. Logarithm of the effective star formation history of the sample shown in Fig. 7. The orange bands show the expected azimuthal variations of the effective star formation history as derived by analysing Milky Way analogues in TNG50.

5.3. Summary

In this study, we modelled the observed K-band absolute magnitude distribution M_K of the sample of massive young stars presented in Zari et al. (2021) to infer the age distribution of stars within ∼3 kpc from the Sun and the effective star formation rate (SFR) of the Milky Way. We divided our volume into spaxels of different sizes and assumed that in each spaxel the M_K distribution depended only on the birthrate. We parametrised the birthrate with a step function with five age intervals to avoid imposing a particularly restrictive functional form, and we derived the parameters characterising the birthrate by comparing our predicted magnitude distribution with the data magnitude distribution by quantifying the likelihood of the data.

We derived the density distribution of stars of different ages in the Galactic plane. We found that stars with ages of < 50 Myr are clustered in compact groups that trace Galactic spiral arms. Over-densities corresponding to spiral arms are, however, also visible in the density distribution of older populations. We derived the star formation history of our sample over the last Gyr by estimating the effective birthrates (i.e., the birthrate that would lead to the observed magnitude distribution) for all the spaxels in our volume, and adding them over all the spaxels.

We found that the age distribution (or total birthrate) is almost uniform between 20 Myr and 1 Gyr, but has increased by almost a factor or 3 in the last ∼10–20 Myr. By using TNG50 simulations of Milky Way analogue galaxies, we derived that azimuthal SFR variations are considerable but smaller than the recent SFR increase. It is not clear which physical mechanisms could be responsible for this increase; hence, we cannot exclude the possibility that it might reflect a local fluctuation, due to the large number of massive star-forming regions within our volume. Based on the assumption that the local value could be extended to the entire Milky Way disc, we estimated the current (last ∼10 Myr) SFR of the entire Milky Way disc to be $3.3_{-0.6}^{+0.7}{M}_{\odot }{\mathrm{yr}}^{-1}$.

Acknowledgments

We thank the referee for their comments that improved the clarity of the manuscript. We are also grateful to D. Weinberg for valuable input. NF was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference number CITA 490888-16] through the CITA postdoctoral fellowship and acknowledges partial support from a Arts & Sciences Postdoctoral Fellowship at the University of Toronto. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research made use of Astropy, (Astropy Collaboration 2013, 2018), matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020). This work would not have been possible without the countless hours put in by members of the open-source community all over the world.

References

Appendix A: Azimuthal birthrate variations

To quantify the birthrate spatial variations, we considered five azimuths and studied the variations of b_best(τ|x, α) as a function of galactocentric radius for the same age intervals as in Fig. 5. Figure A.1 shows that the sharp peaks in birthrate for 5 < τ < 10 Myr coincide with the locations of spiral arm segments, while for the older age ranges the birthrate is lower and does not correlate strongly with the position of spiral arms.

Fig. A.1. Birthrates as a function of galactocentric radius and azimuths. Top left: Same as Fig. 3, the dashed lines correspond to the different azimuths considered. Top centre, right, and bottom: Birthrate as a function of galactocentric radius for three different age intervals at five different azimuths (as indicated in the panels). The colour scheme is the same as in Fig. 5. The thick lines correspond to the average birthrate, and the shaded areas correspond to the 16th and 84th percentiles.

All Tables

All Figures

Fig. 1. Density distribution of the filtered sample of massive hot stars from Zari et al. (2021). In the left panel, the coloured lines represent the spiral arm locations derived by Reid et al. (2019). The right panel highlights the 6 × 6 kpc region considered in this study. Such region is divided in ‘spaxels’ of different sizes, each containing the same number of sources. In both panels, the Sun is in X, Y = (0, 0).

In the text

Fig. 2. Model fits to mock datasets used as tests. Top row: black lines show the true value of the birthrate αi parameters. The orange lines show the estimated αi parameters for three different realisations of the mock datasets generated from the same initial parameters. The orange shaded areas correspond to the 1- and 2-σ levels. Bottom row: black histograms represent the magnitude distribution of three different mock datasets. The thick orange histograms represent the magnitude distribution obtained by using the best fit αi. The dashed histograms represent the magnitude distribution obtained by using the 16th and 84th percentiles for αi.

In the text

	Fig. 3. Density distribution of stars younger than 500 Myr in the Galacic plane, as derived from Eq. (11), i.e. by using the best estimates for the birthrate coefficients, αi, in each spaxel, j.
In the text

Fig. 4. Density distribution of stars with ages of 5 < t < 50 Myr (top left), 50 < t < 100 Myr (top right), 100 < t < 200 Myr (bottom left), and 200 < t < 500 Myr (bottom right). The density distribution changes as a function of time, becoming gradually smoother. Over-arching structures corresponding to spiral arms are visible for all age intervals. The solid black lines indicate the location of the spiral arms from Reid et al. (2019).

In the text

Fig. 5. Distribution of star formation rate surface density as a function of galactocentric radius. The thick lines correspond to the average star formation rate surface density for different age intervals (red: 5 < τ < 10 Myr; orange: 50 < τ < 100 Myr; yellow: 400 < Myrτ < 1 Gyr). The shaded areas correspond to the 16th and 84th percentiles. The left Y-axis shows the star formation rate expressed as the total number of stars born; the right Y-axis shows the more conventional amount of stellar mass formed, which requires the adoption of an IMF (here, Kroupa 2001).

In the text

Fig. 6. Star formation rate density in different galactocentric radii, as a function of time. The left and right Y-axes show the SFR expressed as the total number of stars born and as the stellar mass formed, respectively. The thick lines correspond to the average star formation rate surface density for three intervals in galactocentric radius. The shades areas correspond to the 16th and 84th percentiles. Ages below 5 Myr are not included in our analysis, and thus the corresponding area in the plot is hatched.

In the text

Fig. 7. Effective star formation history of the sample, in units of yr−1 (left y-axis) and M⊙ yr−1 (first right y-axis). The second right y-axis corresponds to the total SFR obtained by multiplying by the ratio f = SFR5 − 11/SFRtot (see Sect. 5). Vertical bars correspond to errors on the parameters. Ages below 5 Myr are not included in our analysis, and thus the corresponding area in the plot is hatched.

In the text

	Fig. 8. Logarithm of the effective star formation history of the sample shown in Fig. 7. The orange bands show the expected azimuthal variations of the effective star formation history as derived by analysing Milky Way analogues in TNG50.
In the text

Fig. A.1. Birthrates as a function of galactocentric radius and azimuths. Top left: Same as Fig. 3, the dashed lines correspond to the different azimuths considered. Top centre, right, and bottom: Birthrate as a function of galactocentric radius for three different age intervals at five different azimuths (as indicated in the panels). The colour scheme is the same as in Fig. 5. The thick lines correspond to the average birthrate, and the shaded areas correspond to the 16th and 84th percentiles.

In the text