11institutetext: Max-Planck-Institut für Astronomie, Königstuhl 17 D-69117
Heidelberg, Germany
22institutetext: Canadian Institute for Theoretical Astrophysics, University
of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada

# Did the Milky Way just light up?
The recent star formation history of the Galactic disc

Eleonora Zari  Neige Frankel 11 22 Hans-Walter Rix 11

(Received -; accepted -)

We map the stellar age distribution ($\lesssim 1$ Gyr) across a
6kpc$\,\times\,$6kpc area of the Galactic disc, in order to constrain our
Galaxy’s recent star-formation history. Our modelling draws on the sample of
Zari et al. (2021) that encompasses all presumed disc OBA stars ($\sim
500,000$ sources) with $G<16$. To be less sensitive to reddening, we do not
forward model the detailed CMD distribution of these stars, but instead the
K-band absolute magnitude distribution, $n(M_{K})$ among stars with $M_{K}<0$
and $T_{\mathrm{eff}}>7000$ K at a certain positions $\@vec{x}$ in the disc as
a step function with five age bins, $b(\tau\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha})$, logarithmically-spaced in age
from $\tau=5$ Myr to $\tau\sim 1$ Gyr. Given a set of isochrones and a Kroupa
(2001) initial mass function, we sample $b(\tau\,|\,\@vec{x},\@vec{\alpha})$
to maximise the likelihood of the data $n(M_{K}\,|\,\@vec{x},\@vec{\alpha})$,
accounting for the selection function. This results in a set of mono-age
stellar density maps across a sizeable portion of the Galactic disc. These
maps show that some, but not all, spiral arms are reflected in overdensities
of stars younger than 500 Myr. The maps of the youngest stars ($<10$ Myr)
trace major star forming regions. The maps at all ages exhibit an outward
density gradient and distinct spiral-like spatial structure, which is
qualitatively similar on large scales among the five age bins. When summing
over the maps’ area and extrapolating to the whole disc, we find an effective
star-formation rate over the last 10 Myr of $\approx
3.3\mathrm{M_{\odot}/yr}$, higher than previously published estimates that had
not accounted for unresolved binaries. Remarkably, our stellar age
distribution implies that the star-formation rate has been three times lower
throughout most of the last Gyr, having risen distinctly only in the very
recent past. Finally we use TNG50 simulations to explore how justified the
common identification of local _age distribution_ with global _star formation
history_ is: we find that the global star-formation rate at a given radius in
Milky-Way-like galaxies is approximated within a factor of $\sim 1.5$ by the
young age distribution within a 6kpc$\,\times\,$6kpc area near $R_{\odot}$.

###### Key Words.:

stars: early-type; Galaxy: structure; Galaxy: disk;

##  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 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 CMD position, or
simply by the fact that they are short-lived. But 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\lesssim 500$ pc) molecular clouds in the Milky Way, nearly
complete lists of young stellar objects (YSOs) are available, whose ages ($t$)
can be derived by measuring their CMD positions or infrared excesses. The SFR
is therefore given by $\langle\dot{M}_{*}\rangle=N_{YSO}\leavevmode\nobreak\
\langle M_{YSO}\rangle/t$, where $<M_{*}>$ is the mean mass of YSOs and is
usually assumed to be $\langle M_{*}\rangle=0.5\,\mathrm{M_{\odot}}$ (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, e.g. Robitaille & Whitney (2010), Chomiuk & Povich (2011),
and Licquia & Newman (2015). Specifically, Robitaille & Whitney (2010) used
the census of young stellar objects (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 $\mathrm{M_{\odot}\,yr^{-1}}$. 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
$<\dot{M_{*}}>=1.9\pm 0.4\leavevmode\nobreak\ \mathrm{M_{\odot}}/yr$. Licquia
& Newman (2015) revisited this work and estimated $<\dot{M_{*}}>=1.65\pm
0.19\leavevmode\nobreak\ \mathrm{M_{\odot}}/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)
corotating with the Sun were not necessarily born within that volume. Second,
we know from external disk galaxies that the instantaneous SRF 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 that any estimate of
$\dot{\mathrm{M}}_{*}\sim\mathrm{M}_{*}(<\Delta t)/\Delta t$ will depend on
the age range $\Delta t$ of the young stars considered.

In this study, we set out to model and determine the age distribution over a
substantive fraction ($6\times 6$ kpc around the Sun) of the Galactic disc,
using upper main sequence, i.e. OB(A), stars. To do this we devise and apply a
non-parametric model111a model for which we do not assume a specific
functional form, but have parameters for 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 explore to which 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) have 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 Section 2. We lay out our methodology in Section
3. Our results, presented in Section 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 Section 5.

Figure 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\times 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)$.

##  2 Data: Luminous Young Disc Stars

The sample that we used in this study consists of $\sim 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 stars. Then, they
applied several cuts in colour-colour space to exclude bright red giant- 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 ($|z|<300$ pc)
and with small vertical velocities $v_{z}$.

Z21’s sample extends to distances of $\sim$5 kpc from the Sun (Fig. 1, left).
However, in this study we focus on a 6kpc$\times$6kpc region of the Galactic
disc, centred on the Sun, shown in Fig. 1 (right), where we can reasonably
assume that our sample, as selected in Eq. 3, is complete. To derive our
spatially resolved SFR map of the Galactic disc, we divide our volume into 256
”spaxels” of different sizes, each containing the same number of stars ($\sim
1000$). We model the magnitude distribution of the sources in each spaxel, as
described in the next Section.

##  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 star dereddening, which
would be difficult if not impossible (but see e.g. Ruiz-Lara et al., 2020).
Therefore, we resort to modelling a more restricted set of observables, 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, as a statistic $n\left(M_{K}\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha}\right)$, the absolute K-band
magnitude distribution of sample members, which we take to be a function of
position in the Galaxy $\@vec{x}$. And we assume that at each $\@vec{x}$ there
is an underlying age or birthrate distribution,
$b\bigl{(}\tau\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha}\bigr{)}$, whose temporal dependence is described by a
set of parameters $\@vec{\alpha}$, 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:

$\displaystyle n\left(M_{K}\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha}\right)=\int d\tau\leavevmode\nobreak\ \int
dm_{0}\leavevmode\nobreak\ \int dT_{\mathrm{eff}}\ \ \times$ $\displaystyle\ \
\ \ \ S_{c}(M_{K},T_{\mathrm{eff}})\ \
p(M_{K},T_{\mathrm{eff}}\leavevmode\nobreak\ |\leavevmode\nobreak\
\tau,m_{0})\ \ \xi(m_{0})\ \ b\bigl{(}\tau\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha}\bigr{)}.$ (1)

In Eq.3 the term $S_{c}(M_{K},T_{\mathrm{eff}})$ reflects the observational
selection function of the sample as specified in Eq. 3;
$p(M_{K},T_{\mathrm{eff}}|\tau,m_{0})$ reflects the probability density of a
star to have absolute magnitude $M_{K}$ and effective temperature
$T_{\mathrm{eff}}$ at a given birth mass $m_{0}$ and age $\tau$; and
$\xi(m_{0})$ is the initial mass function in units of
$\mathrm{mass}\leavevmode\nobreak\ [M_{\odot}]^{-1}$. The integral over age
$\tau$ should in principle extend from 0 to $t_{Hubble}$. In practice, we
consider minimum and maximum ages from $\tau_{min}$ to $\tau_{max}$, as the
very youngest stars are too enshrouded to be in the sample, and stars
$\tau\gtrsim 1$ Gyr will be too faint in their main-sequence phase. In the
next section we describe the single terms of Eq. 3.

###  3.1 Functional forms of the model components

Birthrate:  We specify the birthrate $b\bigl{(}\tau\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha}\bigr{)}$ as a function of time
$\tau$ at a given position in the Galaxy $\@vec{x}$, given parameters
$\@vec{\alpha}$. 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 parameterize it by
step functions in time that can be different in each spatial pixel $\@vec{x}$:

$b\bigl{(}\tau\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha}\bigr{)}=\sum_{\mathrm{age\leavevmode\nobreak\
bin}\leavevmode\nobreak\
i}\alpha_{i}\bigl{(}\@vec{x}\bigr{)}\leavevmode\nobreak\ \chi_{A_{i}}(\tau),$
(2)

where $A_{i}$ are age intervals, the rate coefficients $\alpha_{i}$ are the
elements of $\@vec{\alpha}$ and must satisfy $\alpha_{i}>0$, and the indicator
function $\chi_{A_{i}}=1$ if $\tau\in 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:

$\displaystyle A_{1}$ $\displaystyle=5-10\,\mathrm{Myr}$ $\displaystyle A_{2}$
$\displaystyle=10-30\,\mathrm{Myr}$ $\displaystyle A_{3}$
$\displaystyle=30-100\,\mathrm{Myr}$ $\displaystyle A_{4}$
$\displaystyle=100-300\,\mathrm{Myr}$ $\displaystyle A_{5}$
$\displaystyle=300\,\mathrm{Myr}-1\,\mathrm{Gyr},$

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

Initial mass function:  We use a broken power law Kroupa IMF Kroupa (2001,
2002) whose piecewise components $\xi(m)\approx m^{-\Gamma}$ are:

$\displaystyle\Gamma$ $\displaystyle=0.3\hskip 14.22636pt\mathrm{for}\,\,\
0.01<\mathrm{M/M_{\odot}}<0.08$ $\displaystyle\Gamma$ $\displaystyle=1.8\hskip
14.22636pt\mathrm{for}\,\,\ 0.08<\mathrm{M/M_{\odot}}<0.5$
$\displaystyle\Gamma$ $\displaystyle=2.7\hskip 14.22636pt\mathrm{for}\,\,\
0.5<\mathrm{M/M_{\odot}}<1$ $\displaystyle\Gamma$ $\displaystyle=2.3\hskip
14.22636pt\mathrm{for}\,\,\ 1<\mathrm{M/M_{\odot}}<120.$

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
$\xi(m)$. We discuss in Section 5 how our results would have been affected by
using an alternate, simpler IMF model (e.g., Chomiuk & Povich, 2011;
Robitaille & Whitney, 2010).

Selection function:  The term $S_{c}$ in Eq. 3 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 Section 2 and
Z21 as a function of $M_{K}$ and $T_{\mathrm{eff}}$. We presume that
$S_{c}(M_{K},T_{\mathrm{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_{\mathrm{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_{\mathrm{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. 3 has non-zero
contributions only from parts of the isochrones where both $M_{K}$ and
$T_{\mathrm{eff}}$ are within the “selected” range. In practise, we
approximate the sample selection function as:

$S_{c}(M_{K},T_{\mathrm{eff}})=\cases{1}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
{M_{K}<\min\\{0,M_{K}^{max}(\@vec{x})\\}}\leavevmode\nobreak\
\wedge\leavevmode\nobreak\ T_{\mathrm{eff}}>7000\mathrm{K}\\\
0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{otherwise.}{}$ (3)

The condition $M_{K}<\min\\{0,M_{K}^{max}(\@vec{x})\\}$ arises from the fact
that for the most distant patches $\@vec{x}$ the apparent magnitudes $G$ of
stars with absolute magnitudes $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
derive $M_{K}^{max}(\@vec{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).

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 $\tau$:
$p(M_{K}|\tau)$ (in units of $\mathrm{mag}^{-1}$). To compute $p(M_{K}|\tau)$
we first apply our selection function $S_{c}(M_{K},T_{\mathrm{eff}})$, and
marginalise (integrate) over all $T_{\mathrm{eff}}$ and $m_{0}$. We use 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
sample the Kroupa IMF and assign to each star the appropriate absolute
magnitude $M_{K}$ by finding the nearest isochrone point according to the mass
of the star. Finally, we count the number of stars in bins of width $\Delta
M_{K}=0.25$ mag, between $M_{K,min}$ and $M_{K,max}$ and we normalise such
number for the total number of stars and the bin width $\Delta M_{K}$.
Therefore, we obtain the quantity $n(M_{K},\tau)$ that corresponds to the
integral over the initial mass and effective temperature
$m_{0},T_{\mathrm{eff}}$ of Eq. 3:

$\displaystyle n(M_{K},\tau)$ $\displaystyle=\int dm_{0}\int
dT_{\mathrm{eff}}$ $\displaystyle
S_{c}(M_{K},T_{\mathrm{eff}})\leavevmode\nobreak\
p(M_{K},T_{\mathrm{eff}}\leavevmode\nobreak\ |\leavevmode\nobreak\
\tau,m_{0})\leavevmode\nobreak\ \xi(m_{0}).$ (4)

Since we do not know the individual ages of the stars in bin $\@vec{x}$, we
also need to integrate over isochrone age $\tau$ before comparing with the
data:

$n(M_{K}\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha})=\int_{\tau_{min}}^{\tau_{max}}n(M_{K},\tau)b\bigl{(}\tau\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha}\bigr{)}d\tau,$ (5)

which we do numerically. As the isochrones are logarithmically spaced, it
proves useful to define the logarithmic age $\omega=\log\tau$ and transforming
integration variables to get

$n\bigl{(}M_{K}\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha})=\int_{i,\omega_{min}}^{i,\omega_{max}}n(M_{K}\leavevmode\nobreak\
|\leavevmode\nobreak\ 10^{\omega_{i}})\leavevmode\nobreak\
b(10^{\omega_{i}})\leavevmode\nobreak\
\ln(10)e^{\ln(10)\omega_{i}}\leavevmode\nobreak\ d\omega$ (6)

where we dropped the arguments $\@vec{x}$ and $\@vec{\alpha}$ from $b$ for
conciseness. Now, Eq. 6 represents the predicted $M_{K}$-distribution at a
given position $\@vec{x}$, given the model assumptions about
$b\bigl{(}\tau\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha}\bigr{)}$.

###  3.2 The data likelihood function

We compare the predictions from Eq. 6 to the observed magnitude distribution
$\\{M_{K,i}\\}$ detected in the patch $\@vec{x}$ by quantifying and optimising
the likelihood of the data. We assume that the stars in our data set are
sampled by an inhomogeneous Poisson point process with a rate function
$\lambda=n(M_{K}\leavevmode\nobreak\ |\@vec{x},\@vec{\alpha})$. This leads the
joint probability of the data, given the model parameters $\@vec{\alpha}$ of
$p(N=N_{\star},\\{M_{K,i}\\}|\@vec{\alpha})$ (the likelihood) to be the
product of the two terms specified below. The first term sets the
normalisation, and informs on the total size of the sample,

$\frac{\left[\int n(M_{K}\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha})dM_{K}\right]^{N}}{N!}\exp{\left[-\int
n(M_{K}\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha})dM_{K}\right]}$ (7)

states that the probability of observing N stars in a spaxel follows a
Poissonian distribution. The second term

$\prod_{i}\frac{n(M_{K,i}\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha})}{\int n(M_{K}\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha})dM_{K}},$ (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:

$\mathcal{L}\left(\leavevmode\nobreak\ \\{M_{K,i}\\}\leavevmode\nobreak\
|\leavevmode\nobreak\
\@vec{x},\@vec{\alpha}\right)=\Biggl{(}\prod_{i=1}^{N_{*}}n(M_{K,i}\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha})\Biggr{)}\leavevmode\nobreak\
e^{-N_{*,tot}},$ (9)

where the normalization factor is:

$N_{*}^{tot}(\@vec{x})=\int_{M_{K,min}}^{M_{K,max}=0}n(M_{K}\,|\,\@vec{x},\@vec{\alpha})\leavevmode\nobreak\
dM_{K}.$ (10)

Using uniform priors between 0 and 0.1, we sample the $\alpha_{i}(\@vec{x})$
parameters characterising the birthrate $b\bigl{(}\tau\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha}\bigr{)}$ 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$\sigma$ levels. We call $b\bigl{(}\tau\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha}_{best}\bigr{)}$ the birthrate
obtained by using the best estimates of the $\alpha_{i}$ parameters.

###  3.3 Tests on mock data

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

To test our model, we generate a mock magnitude distribution
$n(M_{K}|\@vec{x},\@vec{\alpha})$ (see Eq. 6) by using a known birth rate
$b\bigl{(}\tau\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha}_{true}\bigr{)}=\sum_{i}\alpha_{i,true}\chi_{A_{i}}(\tau)$
(where we chose here $\@vec{\alpha}_{true}=\@vec{\alpha}_{best}$), and the
same selection function, set of isochrones, and IMF as used in our model. We
then use our fitting procedure described in the previous section to derive the
parameters $\alpha_{i,\mathrm{pred}}$ and compare them with the true values
$\alpha_{i,\mathrm{true}}$. To use this test in a physically meaningful range
of the parameter space and to anticipate the results (Section 4,) the
$\alpha_{i,\mathrm{true}}$ are chosen to reproduce the observed number of
stars in a spaxel at a given position $\@vec{x}$.

Figure 2 (top row) shows the true (black lines) and predicted values (orange
lines) for $\alpha_{i}$, for three realisation of the same mock magnitude
distribution. The shaded areas correspond to $1\sigma$ and $2\sigma$ levels.
Figure 2 (bottom row) shows the corresponding true mock magnitude distribution
(black histogram) and the predicted magnitude distribution for the best fit
$\alpha_{i}$ value (thick solid orange histogram) and the $1\sigma$ values
(dashed thin histograms). The predicted parameters are compatible within
1-2$\sigma$ with the true values and they are always consistent with each
other. The increase in birth rate 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.

##  4 Results

In this Section we present the results obtained by fitting our model to the
data presented in Section 2. We compute the values and uncertainties for the
coefficients $\alpha_{i,j}$ that imply $b(\tau_{i}\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x}_{j},\@vec{\alpha}_{j})$ for all the spaxels
$\@vec{x}_{j}$ in the map of Fig. 1 (right) by sampling the likelihood of Eq.
9. Example values of the $\alpha_{i}$ parameters and their uncertainties are
given in Table 1 and will be available on CDS.

Table 1: Values and uncertainties for three spaxels, and corresponding coordinates. The full result table will be available upon request and on CDS. $x_{\mathrm{min}}$ – $x_{\mathrm{max}}[kpc]$ |  $y_{\mathrm{min}}$ – $y_{\mathrm{max}}$ [kpc] | $\alpha_{\alpha_{16}}^{\alpha_{84}}$ | $\alpha_{\alpha_{16}}^{\alpha_{84}}$ | $\alpha_{\alpha_{16}}^{\alpha_{84}}$ | $\alpha_{\alpha_{16}}^{\alpha_{84}}$ | $\alpha_{\alpha_{16}}^{\alpha_{84}}$
---|---|---|---|---|---|---
|  | 5-10 Myr | 10-30 Myr | 30-100 Myr | 100-300 Myr | 300 Myr - 1 Gyr
2.23–2.99 | 2.57–2.99 | $0.009_{0.008}^{0.01}$ | $0.0006_{0.0002}^{0.0016}$ | $0.00025_{0.00007}^{0.0006}$ | $0.0026_{0.0024}^{0.0029}$ | $0.00009_{0.00002}^{0.0002}$
2.23–2.99 | 1.96–2.57 | $0.011_{0.009}^{0.013}$ | $0.0022_{0.0006}^{0.0046}$ | $0.00033_{0.00009}^{0.0008}$ | $0.0010_{0.0006}^{0.0014}$ | $0.0015_{0.001}^{0.002}$
1.86–2.23 | 2.49–2.99 | $0.008_{0.007}^{0.009}$ | $0.0018_{0.0007}^{0.0034}$ | $0.0007_{0.0002}^{0.0014}$ | $0.0019_{0.0016}^{0.0022}$ | $0.00026_{0.00009}^{0.00052}$

These results are essentially number- or mass-density maps of mono-age stellar
populations (in five age bins) across our surrounding 6kpc$\,\times\,$6kpc
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.

Figure 3: Density distribution of the 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 $\alpha_{i}$ in each spaxel $j$.

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

Figure 3 shows the number surface density of stars
$n(\Delta\tau\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha})$, younger than 500 Myr, defined as:

$n(\Delta\tau\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha})\equiv\frac{1}{\mathrm{Area}(\@vec{x})}\int_{\tau_{min,i}}^{\tau_{max,i}}b_{best}(\tau\leavevmode\nobreak\
|\leavevmode\nobreak\ \@vec{x},\@vec{\alpha}_{best})\leavevmode\nobreak\
d\tau,$ (11)

where the integral of $b(\tau\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha}_{best})$ over time gives the number of stars formed
between $\tau_{min}<\tau<\tau_{max}$ and $\mathrm{Area}(\@vec{x})$ is the
spaxel area in kpc2. 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 with the initial mean mass per star, $\langle
M_{*}\rangle=0.22\,\mathrm{M_{\odot}}$ for our Kroupa IMF.

Figure 4: Density distribution of stars with ages $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).

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

###  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 6kpc$\,\times\,$6kpc box, simply summing over all spaxels in a
certain $R_{GC}$ range:

$n_{R}(\tau\leavevmode\nobreak\ |\leavevmode\nobreak\ R_{GC})\equiv\
\frac{1}{S}\sum_{\mathrm{spaxels\leavevmode\nobreak\ at\leavevmode\nobreak\
R_{GC}}}b(\tau\leavevmode\nobreak\ |\leavevmode\nobreak\
\@vec{x},\@vec{\alpha}_{best}),$ (12)

where $S$ is the total area of these spaxels.

Figure 5 shows the radial distribution of $n_{R}(\tau_{i}\leavevmode\nobreak\
|\leavevmode\nobreak\ R_{GC})$, for three different age intervals $\tau_{i}$:
$5<\tau<10$ Myr, $50<\tau<100$ Myr, and $400\,\mathrm{Myr}<\tau<1$ Gyr. While
the older age bins only show a shallow gradient, the radial density gradient
among stars $\tau<10$ Myr is quite steep.

Figure 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<\tau<10$
Myr; orange: $50<\tau<100$ Myr; yellow: $400<$ Myr $\tau<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 the more conventional amount of stellar mass formed, which requires the
adoption of an IMF (here Kroupa, 2001).

Figure 6 shows the azimuthally averaged mean birth rate parameter
$n_{R}(\tau_{i}|R)$ for different Galactocentric radii. The different
Galactocentric annuli do not only vary significantly in their normalization,
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 about 300 Myr ago.

Figure 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 younger than 5 Myr are not included in our
analysis, and thus the corresponding area in the plot is hatched.

###  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 entire $6\times 6\leavevmode\nobreak\ $ kpc2
area. We define the total ”effective” birth rate $b_{tot}$ as

$b_{tot}(\tau)=\sum_{all\leavevmode\nobreak\ spaxels\leavevmode\nobreak\
\@vec{x}}b(\tau\,|\,\@vec{x},\@vec{\alpha}_{best})$ (13)

by summing over all the $j$ spaxels in the map (and _not_ normalising by the
area). We shall discuss below to what extent this age distribution reflects
the overall _star formation rate_ across the Galactic disc.

Units of $b_{tot}(\tau)$ are # of stars/year, which we convert to the more
conventional $\mathrm{M_{\odot}\,yr^{-1}}$ by multiplying by the mean stellar
mass derived from our assumed IMF $<M>=0.22\,\mathrm{M_{\odot}}$. This is
shown in Fig. 7, where we show $b_{tot}(\tau)$ on the left Y-axis, and
$\mathrm{SFR_{eff}}$ on the first right Y-axis. Both values refer to the
entire $6\times 6\leavevmode\nobreak\ $ kpc2 area, a much larger area than
previous quantitative studies for this age range, but of course not the
“entire disc”.

Fig. 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 $\sim
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.

Figure 7: Effective star formation history of the sample, in units of
$\mathrm{yr^{-1}}$ (left y-axis) and $\mathrm{M_{\odot}\,yr^{-1}}$ (first
right y-axis). The second right y-axis corresponds to the total SFR obtained
by multiplying by the ratio $f=\mathrm{SFR_{5-11}/SFR_{tot}}$ (see Section 5).
Vertical bars correspond to errors on the parameters. Ages younger than 5 Myr
are not included in our analysis, and thus the corresponding area in the plot
is hatched.

##  5 Discussion and Summary

In the previous sections we have estimated the recent star formation history
of the extended solar neighbourhood and found that it has recently increased
by a factor of $\sim$3\. We now estimate the global star formation rate of the
Galactic disc and we 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 Section 5.2), we can derive the “total” SFR of
the Milky Way. To do so, we consider the youngest bin, $\tau\lesssim$ 10 Myr,
as the best approximation to the “current SFR” and we extend
$\mathrm{SFR_{eff,5-10Myr}}$ 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 well described 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):

$\mathrm{SFR}_{tot}=2\pi\int_{R_{min}}^{R_{max}}S_{0}e^{-r/L}dr.$ (14)

We chose $R_{min}=4.5\,\mathrm{kpc}$ and $R_{max}=15\,\mathrm{kpc}$ following
Reed (2005) and Davies et al. (2011) and thus exclude the central region of
the Galaxy. In the radial range covered by our data, $\sim 5$ and $\sim 11$kpc
from the Galactic centre, the total SFR is:

$\mathrm{SFR}_{5-11\leavevmode\nobreak\
\mathrm{kpc}}=\Delta\phi\int_{5}^{11}S_{0}e^{-r/L}dr,$ (15)

where $\Delta\phi\approx 0.7$ rad corresponds to the range of azimuth angles
covered. The ratio $f=\mathrm{SFR_{5-11}}/\mathrm{SFR_{tot}}\approx 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}\,\mathrm{M_{\odot}/yr}$ (see Figure 7).

Recent analogous estimates of the star formation rate of the Milky Way are
between $1.65\pm 0.19\,\mathrm{M_{\odot}\,yr^{-1}}$ (Licquia & Newman, 2015)
and $1.9\pm 0.4\,\mathrm{M_{\odot}\,yr^{-1}}$ (Chomiuk & Povich, 2011),
silimar 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.

Assumed IMF. In their Section 3, Chomiuk & Povich (2011) review 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), infra-red emission and YSO counts
(Robitaille & Whitney, 2010; Davies et al., 2011), or massive stars 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 as reference the IMF from IMF (Kroupa &
Weidner, 2003, KW03): $\xi(m)\propto m^{-\Gamma}$ with $\Gamma=1.3$ for masses
between 0.1 and 0.5 $\mathrm{M_{\odot}}$ and $\Gamma=2.3$ for masses between
0.5 and 100 $\mathrm{M_{\odot}}$. 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 Section 3. The choice of the IMF impacts, of course, the relation
between the observed and modelled luminosity function
$n(M_{K}\leavevmode\nobreak\ |\leavevmode\nobreak\ \tau)$ and the implied
“underlying” total stellar mass (Eqs. 3.1 & 6). In our formulation of the
problem, this affects the inference of the birthrate parameter $\alpha_{i}$,
and the mean stellar mass that is used to convert the SFR to units of
$\mathrm{M_{\odot}}\,yr^{-1}$, which for the KW03 IMF is
$<M>=0.63\,\mathrm{M_{\odot}}$. To quantify the impact of the choice of IMF on
our results, we re-run the analysis using KW03 IMF. This yields an analogous
estimate for the total SFR of $\leavevmode\nobreak\ 2.3\pm
0.4\,\mathrm{M_{\odot}}\,yr^{-1}$, which is consistent with the estimates of
Chomiuk & Povich (2011) and slightly higher than Licquia & Newman (2015).

Space distribution. There are strong assumptions behind the space distribution
of different tracers used to study the global SFR or the Milky Way (including
ours, see Section 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 proposed by Reed (2005) and Davies et al. (2011).

Mean Age of SFR Tracers. Finally, the lifetimes of the tracers used in other
studies are different. Robitaille & Whitney (2010)’s estimate
($0.68-1.45\,\mathrm{M_{\odot}yr^{-1}}$) assume a YSO lifetime of $\approx 2$
Myr, which is younger than the minimum age considered here (5 Myr). Murray &
Rahman (2010) use a lifetime for HII regions of 3.7 Myr following McKee &
Williams (1997), slightly longer than the 3 Myr used by Mezger (1978), and
again younger than our minimum age. Reed (2005) assigns lifetimes from
$\approx 3.5$ Myr to $\approx 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\times$ 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 $\sim$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_{\odot}$/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 with higher confidence the recent enhancement in star
formation activity of the Galactic disc 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 $\tau<500$ Myr. Consequently, Fig. 4 represents evidence
that those spiral arms that show up in all our monoage 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) compare the distribution of molecular gas (tracing star formation)
and $3.6\,\mu 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 e.g.
from masers, are reflected in stellar mass overdensities in our maps among
stars with $\tau\lesssim 0.5$ Gyrs (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 $\approx$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 use the TNG50 cosmological simulation (Pillepich
et al., 2019; Nelson et al., 2019), whose numerical resolution reaches that of
zoom-in simulations (reaching a baryonic mass of $8.5\times 10^{4}M_{\odot}$).
IN TNG50, we select Milky Way-like galaxies based on their stellar mass,
diskyness and morphology. Specifically, we require:

  * •

$10^{10.5}\leq M_{\star}/M_{\odot}\leq 10^{11.2}$,

  * •

over 40% of stars are on near circular orbits (with circularities $\geq 0.7$,
Genel et al., 2015),

  * •

have a central bar,

  * •

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

leading to $N=78$ galaxies. We centre and align a Cartesian coordinate system
with the vertical axis aligned with the angular momentum of the disc. To mimic
our observed sample, we place a solar neighbourhood at
$(R,z,\varphi)=(8,0,\varphi)$ and select stars in a cylinder centred on the
Sun with a 3 kpc radius. To investigate how the location $\varphi$ of the Sun
influences the results, we vary the azimuth of the Sun and measure 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 derive the distribution of

$\log_{10}\left(\frac{\dot{M}_{\star}/dS(\mathrm{patch})}{\dot{M}_{\star}/dS(\mathrm{annulus})}\right),$

where the denominator is the number of stars per age interval in the pacth,
and the denominator 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 $\sim 1.5$. These
variations are shown in Fig.8 are 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 star formation rate. It is unlikely that the Sun has just moved into a
patch of higher star formation activity of late.

Figure 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 have 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 $\sim 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 birth rate. 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 younger than $<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 $\sim$ 10-20 Myr. Such an increase can be related to the large number of
massive star-forming regions within our volume. By using TNG50 simulations of
Milky Way analogue galaxies, we derived that azimuthal SFR variations are
considerable but smaller than the recent SFR increase. Based on the assumption
that the local value could be extended to the entire Milky Way disk, we
estimated the _current_ (last $\sim 10$ Myr) SFR of the entire Milky Way disk
to be $3.3_{-0.6}^{+0.7}\,\mathrm{M_{\odot}/yr}$.

###### Acknowledgements.

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 et al., 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.

##  Appendix A Azimuthal birthrate variations

To quantify the birthrate spatial variations, we consider five azimuths and we
study the variations of $b_{best}(\tau|\@vec{x},\@vec{\alpha})$ as a function
of Galactocentric radius for the same age intervals as in Fig. 5. Figure 9
shows that the sharp peaks in birthrate for $5<\tau<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.

Figure 9: Top left: same as Fig. 3, the dashed lines correspond to different
azimuths. Top center, 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.

## References

  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bouwens et al. (2011) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2011, ApJ, 737, 90
  * Bovy et al. (2012) Bovy, J., Rix, H.-W., Liu, C., et al. 2012, ApJ, 753, 148
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Chen et al. (2015) Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068
  * Chiappini et al. (1997) Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765
  * Chomiuk & Povich (2011) Chomiuk, L. & Povich, M. S. 2011, AJ, 142, 197
  * Cignoni et al. (2006) Cignoni, M., Degl’Innocenti, S., Prada Moroni, P. G., & Shore, S. N. 2006, A&A, 459, 783
  * Davies et al. (2011) Davies, B., Hoare, M. G., Lumsden, S. L., et al. 2011, MNRAS, 416, 972
  * de Rossi et al. (2009) de Rossi, M. E., Tissera, P. B., De Lucia, G., & Kauffmann, G. 2009, MNRAS, 395, 210
  * Diehl et al. (2006) Diehl, R., Halloin, H., Kretschmer, K., et al. 2006, Nature, 439, 45
  * Eskridge et al. (2002) Eskridge, P. B., Frogel, J. A., Pogge, R. W., et al. 2002, ApJS, 143, 73
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Genel et al. (2015) Genel, S., Fall, S. M., Hernquist, L., et al. 2015, ApJ, 804, L40
  * González Delgado et al. (2016) González Delgado, R. M., Cid Fernandes, R., Pérez, E., et al. 2016, A&A, 590, A44
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90
  * Ibata et al. (1994) Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Nature, 370, 194
  * Kauffmann & Haehnelt (2000) Kauffmann, G. & Haehnelt, M. 2000, MNRAS, 311, 576
  * Kennicutt & Evans (2012) Kennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531
  * Kreckel et al. (2018) Kreckel, K., Faesi, C., Kruijssen, J. M. D., et al. 2018, ApJ, 863, L21
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Kroupa (2002) Kroupa, P. 2002, Science, 295, 82
  * Kroupa & Weidner (2003) Kroupa, P. & Weidner, C. 2003, ApJ, 598, 1076
  * Lada et al. (2010) Lada, C. J., Lombardi, M., & Alves, J. F. 2010, ApJ, 724, 687
  * Laporte et al. (2018) Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286
  * Licquia & Newman (2015) Licquia, T. C. & Newman, J. A. 2015, ApJ, 806, 96
  * Madau et al. (1996) Madau, P., Ferguson, H. C., Dickinson, M. E., et al. 1996, MNRAS, 283, 1388
  * McKee & Williams (1997) McKee, C. F. & Williams, J. P. 1997, ApJ, 476, 144
  * Meidt et al. (2021) Meidt, S. E., Leroy, A. K., Querejeta, M., et al. 2021, ApJ, 913, 113
  * Mezger (1978) Mezger, P. O. 1978, A&A, 70, 565
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
  * Murray & Rahman (2010) Murray, N. & Rahman, M. 2010, ApJ, 709, 424
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2
  * Pillepich et al. (2019) Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196
  * Poggio et al. (2021) Poggio, E., Drimmel, R., Cantat-Gaudin, T., et al. 2021, A&A, 651, A104
  * Reed (2005) Reed, B. C. 2005, AJ, 130, 1652
  * Reid et al. (2019) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131
  * Rix & Bovy (2013) Rix, H.-W. & Bovy, J. 2013, A&A Rev., 21, 61
  * Rix et al. (2021) Rix, H.-W., Hogg, D. W., Boubert, D., et al. 2021, AJ, 162, 142
  * Robitaille et al. (2008) Robitaille, T. P., Meade, M. R., Babler, B. L., et al. 2008, AJ, 136, 2413
  * Robitaille & Whitney (2010) Robitaille, T. P. & Whitney, B. A. 2010, ApJ, 710, L11
  * Ruiz-Lara et al. (2020) Ruiz-Lara, T., Gallart, C., Bernard, E. J., & Cassisi, S. 2020, Nature Astronomy, 4, 965
  * Rybizki et al. (2022) Rybizki, J., Green, G. M., Rix, H.-W., et al. 2022, MNRAS, 510, 2597
  * Sana et al. (2014) Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15
  * Smith et al. (1978) Smith, L. F., Biermann, P., & Mezger, P. G. 1978, A&A, 66, 65
  * Tang et al. (2014) Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Williams et al. (2011) Williams, B. F., Dalcanton, J. J., Johnson, L. C., et al. 2011, ApJ, 734, L22
  * Zari et al. (2021) Zari, E., Rix, H. W., Frankel, N., et al. 2021, A&A, 650, A112


# ATLAS Tile calorimeter calibration and monitoring systems

Marija Marjanović
on behalf of the ATLAS collaboration
Laboratoire de Physique de Clermont
Université Clermont-Auvergne, CNRS-IN2P3
Aubière, France
<EMAIL_ADDRESS>
###### Abstract

The††footnotetext: Copyright 2018 CERN for the benefit of the ATLAS
Collaboration. Reproduction of this article or parts of it is allowed as
specified in the CC-BY-4.0 license ATLAS Tile Calorimeter (TileCal) is the
central section of the hadronic calorimeter of the ATLAS experiment at LHC.
This sampling calorimeter uses steel plates as absorber and scintillating
tiles as active medium. The light produced by the passage of charged particles
is transmitted by wavelength shifting fibers to photo-multiplier tubes (PMTs),
located in the outer part of the calorimeter. The readout is segmented into
about 5000 cells, each one being read out by two PMTs in parallel. To
calibrate and monitor the stability and performance of the full readout chain
during the data taking, a set of calibration sub-systems is used. The TileCal
calibration system comprises Cesium radioactive sources, laser, charge
injection elements, and an integrator based readout system. Combined
information from all systems allows to monitor and to equalize the calorimeter
response at each stage of the signal evolution, from scintillation light to
digitization. Calibration runs are monitored from a data quality perspective
and used as a crosscheck for physics runs. Data quality in physics runs is
monitored extensively and continuously. Any problems are reported and
immediately investigated. The data quality efficiency achieved was 99.6% in
2012, 100% in 2015, 98.9% in 2016 and 99.4% in 2017.

Based on LHC Run-I experience, all calibration systems were improved for Run-
II. TileCal performance during LHC Run-II, (2015-2017), is discussed,
including calibration, stability, absolute energy scale, uniformity and time
resolution. Results show that the TileCal performance is within the design
requirements and has given essential contribution to reconstructed objects and
physics results.

###### Index Terms:

calorimeter, calibration, monitoring

##  I Introduction

The ATLAS detector [1] is a multipurpose particle detector at the Large Hadron
Collider (LHC) at CERN. It is designed to identify and measure the properties
of elementary particles created in proton-proton ($pp$) or heavy ions
collisions. It is composed of several sub-detectors, one being the Tile
Calorimeter (TileCal) [2], [3], a hadronic calorimeter which allows for the
measurements of jets, hadronically decaying tau leptons and missing transverse
energy. TileCal is also used in the Level 1 trigger. Figure 1 shows the
electromagnetic and hadronic calorimeters of ATLAS. TileCal surrounds the
Liquid Argon (LAr) barrel electromagnetic and endcap hadronic calorimeters.

TileCal is a sampling calorimeter: it consists of an absorber layer made of
iron plates, that are used to stop the particles, and a detection layer made
of scintillating tiles. It has around 5000 cells each one read by two fibers,
making a total of 1000 readout channels. TileCal covers the central region of
ATLAS ($|\eta|<1.7$). It consists of one barrel, and two extended barrels in
$\eta$, and 64 modules in $\phi$. Each cylinder has three longitudinal layers
(A, B/C and D). One additional layer (E) covers gap and crack regions. Figure
2 shows the TileCal cell and scintillator structure.

Light is produced in the scintillating tiles when the particle passes through
it. This light is then converted into electric currents by photomultiplier
tubes (PMTs). Signal from the PMTs is then shaped and amplified using two
gains (high and low gain) with a ratio 1:64, with a 10-bit ADC by digitizers.
Signals are sampled and digitized at 40 MHz. Events are retained if they pass
Level 1 trigger. There is also an integrator system that measures the
integrated current from the PMTs.

Figure 1: ATLAS inner detector and calorimeters. Tile Calorimeter consists of
one barrel and two extended barrel sections and surrounds the LAr barrel
electromagnetic and endcap hadronic calorimeters. Figure 2: Schematic showing
the TileCal cell and scintillator structure.

##  II TileCal Calibration systems

The signal coming from each channel is converted into energy and that relation
has to be carefully calibrated and monitored. There are several factors that
affect it, like stability of the high voltage supply, the stress of the PMTs,
the failures of the readout electronics and also the ageing of the optics. In
order to monitor the effect of all these changes, several calibration systems
are used [4] as shown in Figure 3. The paths of different calibration and
monitoring systems are partially overlapping, allowing for cross-checks and an
easier identification of component failures. The calibration constants are
being calculated and applied to convert measured ADC counts into energy. The
reconstructed energy of each channel, E(GeV), is derived from the raw
response, A(ADC), as follows:

$E(GeV)=A(ADC)\times C_{ADC\to pC}\times C_{pC\to GeV}\\\ \times
C_{Cesium}\times C_{Laser}$ (1)

where $C_{ADC\to pC}$ stands for conversion of the ADC counts to charge in pC.
$C_{pC\to GeV}$ is the relation between the detector response, in pC, and the
energy deposit, in GeV. Finally, $C_{Cesium}$ and $C_{Laser}$ are used to
maintain the energy scale constant. While $C_{pC\to GeV}$ was fixed during
dedicated test beam campaigns, the remaining calibration constants are
provided by individual systems.

Figure 3: Flow diagram of the readout signal path of different TileCal
calibration systems.

###  II-A Cesium calibration

The calibration of optic components of TileCal and PMTs is performed using a
movable Cesium radioactive $\gamma$ source [5]. The source of Cesium passes
through the whole calorimeter and emits $\gamma$ rays of well known energy
(662 keV). It uses an integrator read-out system, a different readout system
than the one used for physics. The cesium system is used for checking the
quality of the optics and full system readout (scintillators, fibers, PMTs and
electronics), and to equalize the response of all read-out channels. It also
monitors the cell electromagnetic scale in time.

The Cesium system was upgraded for the Run-II of the LHC. This included
improvements of stability and safety (new water storage system, lower pressure
and precise water level metering). The precision of single channel response is
now better than 0.3%. Figure 4 shows the variation of TileCal response over a
period of 8 years measured in Cesium calibration runs averaged over all cells
in a given sample [6].

Figure 4: Variation of TileCal response over a period of 8 years measured in
Cesium calibration runs averaged over all cells in a given sample.

Figure 5 presents the change of the signal in TileCal cells between February
2015 and January 2018 as a function of $\eta$ in the three longitudinal
samplings. Signal is measured in Cesium calibration runs. Every point
represents average of 64 cells (over $\phi$). The largest drift is observed
for A-cells.

Figure 5: Change of the signal in TileCal cells between February 2015 and
January 2018 as a function of $\eta$ in three longitudinal samplings.

###  II-B Laser calibration

The calibration of PMT gains is performed twice a week using a laser
calibration system [7] during dedicated runs. The laser system measures the
drift seen in PMT response with respect to the last cesium scan. It sends
light pulses to the PMTs with a wavelength close to the one of physical
signals (532 nm) in the absence of collisions. Controlled amount of light is
sent to the PMTs through $\sim$400 fibers. The mean gain variation of the 9852
TileCal channels is computed cell by cell. For each cell, the gain variation
is defined as the mean of the gaussian function that fits the gain variation
distribution of the channels associated to this cell. A total of 64 modules in
$\phi$ are used for each cell while known pathological channels are excluded.
The observed down-drift mostly affects cells at inner radius, that are the
cells with higher current. Maximal drift is observed in E and A cells which
are the cells with the highest energy deposits.

Figure 6 shows the mean gain variation (in %) in the ATLAS TileCal PMTs that
read the signal deposited in each channel, as a function of $\eta$ and radius,
between the 24 May 2016 and the 27 October 2016 (given as an example).

Figure 6: The mean gain variation (in %) in the ATLAS TileCal PMTs that read
the signal deposited in each channel, as a function of $\eta$ and radius,
between the 24 May 2016 and the 27 October 2016 (given as an example).

Figure 7 depicts the mean gain variation (in %) of the channels of the ATLAS
Tile calorimeter as a function of the polar angle ($\phi$) and the layer. The
PMT mean gain variation is defined as the mean of the gaussian function that
fits the gain variation distribution of the channels, arranged into $\phi$
bins. One $\phi$ bin corresponds to one module per partition. The channels
with known problems are not taken into account.

Figure 7: The mean gain variation (in %) in the ATLAS TileCal PMTs that read
the signal deposited in each channel, as a function of $\phi$ and radius,
between the 24 May 2016 and the 27 October 2016 (given as an example).

The laser system underwent an upgrade for Run-II of the LHC. The electronics
and optical components were improved, which led to better control of the
emitted light.

###  II-C Charge injection system

Charge Injection System (CIS) injects a signal of known charge and measures
the electronic response. It checks the full ADC range (0-800 pC) and has two
gains for each PMT (low gain and high gain). CIS calibration is taken at least
twice a week and it measures the pC/ADC conversion factor and corrects for
non-linearities in low gain. It is also used to calibrate analog Level 1
trigger.

Figure 8 shows the detector-wide CIS calibration constant averages of all the
high gain ADC’s for each CIS calibration run from 15 May 2017 to 05 December
2017 (given as an example), plotted as black circles. The CIS constants from a
typical channel (LBC19, Channel 13) are additionally plotted as blue triangles
for comparison. The RMS values on the plot are indicative of the fluctuation
present in calibrations. In addition, there is a 0.7% systematic uncertainty
present in individual calibrations, represented by the blue error bars.
Problematic channels are not included when calculating the mean.

Figure 8: CIS calibration constant averages of all the high gain ADC’s for
each CIS calibration run from 15 May 2017 to 05 December 2017, plotted as
black circles. The CIS constants from a typical channel (LBC19, Channel 13)
are additionally plotted as blue triangles for comparison.

The CIS measurements have a precision of 0.7% and the average of all channels
doesn’t change by more than 0.03%.

###  II-D Integrator system

The integrator system [8] integrates PMT signals over a large time window
($\sim 10\mu s$). It monitors and measures the response to the source during
Cesium scans, and during physics runs, it measures the detector response to
the minimum bias events. The minimum bias events are soft parton interactions
that dominate the high energy $pp$ collisions. The integrator system also
provides an additional way to monitor the instantaneous luminosity in ATLAS as
it allows for an independent measurement given an initial calibration
(luminosity coefficient). This is possible as the measured currents are
linearly dependent on the instantaneous luminosity. The response stability is
used to produce calibration constants in absence of Cesium calibration and to
calibrate E-cells and MBTS (Minimum Bias Trigger Scintillator cells), which
are not instrumented by Cesium system. The integrator system monitors the full
optical chain.

Figure 9 presents the average current vs. instantaneous luminosity measured by
the integrator system during the $pp$ collision data taking in 2015. The
currents were measured by PMT 17 of module 22 in the C-side Extended Barrel.
This channel belongs to cell D5, which is stable within 1-2% over the 2015
data taking period. The considered runs are spread over the whole data taking
period and fulfil the quality requirements used for physics analyses. A linear
fit has been performed. The fit parameters are given in the upper panel of the
plot. The lower panel shows the ratio of the current over the fitted function.
A good linear description of the current vs instantaneous luminosity can be
observed. The stability of the integrator gains for each channel is better
than 0.05% and the average stability is better than 0.01%.

Figure 9: Average current (measured by PMT 17 of module 22 in the C-side
Extended Barrel) vs. instantaneous luminosity measured by the integrator
system during the $pp$ collision data taking in 2015. The dashed region
presents the integrated luminosity delivered to the ATLAS detector.

The laser and integrator systems are used in parallel to monitor TileCal. The
integrator system monitors PMT response drift and scintillator ageing, while
the laser system monitors only PMT response drifts. Figure 10 shows the
variation of the response to integrator (labeled Minimum Bias) and laser
systems for cell A13 in the inner layer of the Extended Barrel, covering the
region $1.2<|\eta|<1.3$, as a function of time. The response variation is
derived with respect to a reference cell D6 ($1.1<|\eta|<1.3$), which has
exhibited less than 1% drift throughout the 2017 collision period. Each
Minimum Bias point represents the average of the response variation of a
subset of A13 channels (one of the most irradiated cells), corresponding to
the central 80% of the total distribution. Each laser point represents the
average of all channels. Integrator and laser data cover the period from the
beginning of June to halfway November 2017. The response variation versus time
measured by the integrator system has been normalized to the response
variation measured by the laser system on June 12th, corresponding to the
first point in the plot. The integrated luminosity is the total delivered
during the proton-proton collision period of 2017\. As already observed in
previous years the down-drifts of the PMT gains (seen by laser) coincide with
the collision periods, while up-drifts are observed during machine development
periods. The same response could be measured by Cesium, but the integrator
measurement is used, as it is possible to do it with better time resolution
and throughout the year. In 2016 and 2017 a systematic difference is observed
between the laser and integrator systems and it is attributed to scintillator
ageing due to irradiations. Therefore an extra calibration from the integrator
system is applied on some channels.

Figure 10: Variation of the response to integrator (labeled Minimum Bias) and
laser systems for cell A13 (one of the most irradiated cells) in the inner
layer of the Extended Barrel, covering the region $1.2<|\eta|<1.3$, as a
function of time.

Figure 11 shows the relative light yield $I/I_{0}$ of scintillators and
wavelength- shifting fibres for the A13 cell as a function of the deposited
dose for the year 2015-17. The relative light yield is derived from the
difference in the response to minimum bias ($\Delta R_{MB}$) events and laser
pulses ($\Delta R_{Las}$) and interpreted as the light yield loss in
scintillators and fibres due to irradiation, and it is defined as
$I/I_{0}=1+(\Delta R_{MB}-\Delta R_{Las})/100\%$. The response variation is
derived with respect to a reference cell D6. Minimum bias and laser data cover
the 2015-17 $pp$ collision period, and the integrated luminosity is the total
delivered during the period. Vertical error bars represent the statistical and
systematic uncertainties of the measurement. The horizontal bars represent the
RMS of the different dose values within the cell volume and do not include
systematics. For the nominal dose value the average is taken. The black
vertical line represents the expected dose by the end of Run-III (450 fb-1).

Figure 11: Relative light yield $I/I_{0}$ of scintillators and wavelength-
shifting fibres for the A13 cell as a function of the deposited dose.

###  II-E Time calibration

A precise time calibration is crucial as it is needed to reconstruct the
energy recorded in each cell. It is used to set the phase so that a particle
that is travelling from the collision point at the light speed would give a
signal whose measured time would be equal to zero. The timing information is
also used in some physics analysis, for example those that search for
R-hadrons, using time of flight. Time calibration is calculated using jets and
is monitored during data taking using laser system. Laser pulses are sent to
TileCal during empty bunch crossing (1-2 Hz frequency) to calibrate timing.

Figure 12 shows the cell time resolution in jet events ($pp$ collision data at
$\sqrt{s}$=13 TeV with 25 ns bunch spacing) as a function of the energy
deposited in Long Barrel cells. Run-to-run differences are accounted for. All
cells belonging to reconstructed jets with $p_{T}>$20 GeV are considered,
after the usual event and jet cleaning procedures are applied. The closed
circles correspond to Gaussian $\sigma$, the open squares indicate the RMS of
the underlying time distributions. The resolution in Long Barrel is slightly
better than in Extended Barrel cells (even if E-cells are excluded), due to
physically smaller cells in Long Barrel compared to Extended Barrel.

Figure 12: Cell time resolution in jet events ($pp$ data at $\sqrt{s}$=13 TeV,
25 ns bunch spacing) as a function of the energy deposited in Long Barrel
cells.

##  III Performance of the TileCal

###  III-A Noise

The control of the noise also plays an important role in physics measurements
.The total noise in each cell of the calorimeter comes from two sources:
electronic noise and noise coming from the pile-up. Electronic noise is
controlled by measurements taken in the dedicated runs when there is no signal
in the detector. Pile up contribution to noise comes from the multiple
interactions that take place at the same bunch crossing, or from the events
from the previous/ following bunch crossings. Electronic noise stays at the
level below 20 MeV for most of the cells. Noise is measured regularly with
calibration runs. In 2014 new power supplies were installed giving better
electronics stability and more gaussian noise shape.

Figure 13 shows the noise distributions in different TileCal cells as a
function of $\eta$ for zero-bias data collected in 2016 (runs 310247, 310249,
310341, 310370) and Pythia 8 Monte Carlo simulation with A3 Minimum Bias tune
at a center-of-mass energy of 13 TeV with a bunch spacing of 25 ns selected
with an average number of interactions per bunch crossing of 30. The noise was
estimated as the standard deviation of the measured cell energy distribution.
The data (Monte Carlo) is represented with closed (open) markers. The cells of
Layer A, Layer BC, Layer D and Layer E are shown with different colors.

Figure 13: Noise distributions in different TileCal cells, represented as a
function of $\eta$.

Figure 14 presents the TileCal noise dependence with the average number of
interactions per bunch crossing, $\langle\mu\rangle$, for several zero-bias
data runs of 2016 with 25 ns bunch spacing at a center-of-mass energy of 13
TeV (runs 310247, 310249, 310341, 310370) and Pythia 8 Monte Carlo simulation.
The noise was estimated as the standard deviation of the energy distribution
per cell. The data (Monte Carlo) is represented with closed (open) markers.
The noise is shown for one TileCal tower with cells A5, BC5 and D2 in
$0.4<\eta<0.5$. Only $\langle\mu\rangle>15$ are shown in data were the
reliable statistics is observed.

Figure 14: Noise dependence with the average number of interactions per bunch
crossing, $\langle\mu\rangle$.

###  III-B Data quality

TileCal channels are monitored for hardware issues, timing offsets and
miscalibrations. Problematic channels are identified and masked. During
maintenance periods all identified issues are fixed. Data Quality efficiency
was 100% in 2015, 98.8% in 2016 and 99.4% in 2017. Figure 15 presents the
evolution of cell masking in TileCal. It shows the percentage of all cells and
channels in the detector that are masked as a function of time starting from
December 2010. The recent detector status corresponding to the 2nd June, 2018
is noted. The hatched area represents the maintenance period of the detector.
The legend includes the amount of masked cells (0.89%) and masked channels
(1.46%) at that time.

Figure 15: Evolution of Cell Masking in TileCalorimeter.

Figure 16 shows the number of cells masked in each $\eta-\phi$ bin. The bin is
chosen to be $0.1\times 0.1$, corresponding to the calorimeter tower
definition. TileCal cells can have one or two readout channels. The analog
signal from each channel is then amplified in separate high (HG) and low (LG)
branches and digitized by two ADCs. The bin color corresponding to integer
numbers is the integrated number of cells in calorimeter tower. The non-
integer numbers mean the integrated fraction of masked ADCs in tower.

Figure 16: Representation of the number of TileCal masked cells in $\eta-\phi$
plane.

###  III-C Single particle response

An important characteristic of the TileCal is the ratio of energy at
electromagnetic scale to track momentum $\langle E/p\rangle$ for isolated
charged hadrons in minimum bias events. It is used to evaluate calorimeter
uniformity and linearity during data taking.

Figure 17 shows the TileCal response to single isolated charged hadrons,
characterized by energy over momentum ($E/p$), measured using $1.6~{}{\rm
fb^{-1}}$ of proton-proton collision data at 13 TeV collected in 2015. The
fraction of events with additional simultaneous proton-proton collisions
(pile-up) was negligible. The data is compared to simulated events generated
using PYTHIA 8.186. The energy is reconstructed from topological clusters
matched to a track in a cone of $\Delta R<0.2$. Tracks are required to pass
minimum quality criteria and $p>2~{}{\rm GeV}$. To reduce contamination from
neutral hadrons and muons, the energy deposited in the electromagnetic
calorimeter is required to be less than $1~{}{\rm GeV}$, and the fraction of
energy deposited in the TileCal at least $70\%$. Black dots represent data and
the blue line represents simulation. Statistical uncertainties are shown for
data and simulation. The lower plot shows the ratio of data to simulation.
Good agreement between data and MC prediction is observed. The average of the
$E/p$ distribution is 0.67 (below one), as expected for a sampling non-
compensating calorimeter.

Figure 17: Calorimeter response to single isolated charged hadrons.

Figure 18 shows the TileCal response to single isolated charged hadrons,
characterized by the mean of the energy over momentum ($\langle E/p\rangle$)
as a function of momentum, integrated over the pseudo-rapidity and $\phi$
coverage of the calorimeter, measured using $1.6~{}{\rm fb^{-1}}$ of $pp$
collision data at 13 TeV collected in 2015. The fraction of events with
additional simultaneous proton-proton collisions (pile-up) was negligible. The
data is compared to simulated events generated using PYTHIA 8.186.

Figure 18: TileCal response to single isolated charged hadrons, characterized
by the mean of the energy over momentum ($\langle E/p\rangle$) as a function
of momentum, integrated over the pseudo-rapidity and $\phi$ coverage of the
calorimeter, measured using $1.6~{}{\rm fb^{-1}}$ of $pp$ collision data at 13
TeV collected in 2015.

###  III-D Muons

Muons from cosmic rays are used to study in situ the electromagnetic energy
scale and intercalibration of Tile cells. A good energy response uniformity
between calorimeter cells is observed, and the response non-uniformity in
$\eta$ is measured to be below 5%.

Figure 19 presents the profile of energy deposition in the LB-BC layer as a
function of the track $\phi$-coordinate impact point, obtained using 2015
cosmic data. The average response over all central region cells in each module
is shown by symbols of different colors. $\phi$ = -1.57 corresponds to
vertical track. Only statistical uncertainties are shown.

Figure 19: Profile of energy deposition in the LB-BC layer as a function of
the track $\phi$-coordinate impact point.

###  III-E Jet performance

Good description of the cell energy distribution and of the noise in the
calorimeter is crucial for the building of topoclusters which are used for jet
and missing transverse energy reconstruction

Figure 20 shows the distributions of the energy deposited in the TileCal cells
from collision data at $\sqrt{s}$ = 13 and 0.9 TeV superimposed with Pythia
minimum bias Monte Carlo and randomly triggered events from filled and empty
bunch crossings. A Level 1 minimum bias trigger scintillator (MBTS) trigger is
required for the 13 TeV data. To ensure exactly one interaction has occurred
per bunch crossing, only events having a single reconstructed primary vertex
are selected. Each distribution is normalized to unit area. Good agreement is
observed in Tile cell energy distribution.

Figure 20: Distributions of the energy deposited in the TileCal cells.

Figure 21 shows the jet response ratio of the data to the Monte Carlo
simulation as a function of $p_{T}$ for three in situ techniques combined to
determine the energy scale correction: Z+jet (squares), $\gamma$+jet (full
triangles) and multijet (empty triangles). The error bars showing the
statistical and the total uncertainties are presented. The total uncertainty
is derived by adding in quadrature statistical and systematic uncertainties.
The results are shown for anti-$k_{T}$ jets with distance parameter of R = 0.4
calibrated with the electromagnetic (EM) and jet energy scale (JES) scheme
followed by $\eta$ intercalibration and using the 2016 dataset. The result of
the combination of the in situ techniques is shown as the dark line. The outer
band indicates the total uncertainty resulting from the combination of in situ
techniques, while the inner dark band presents the statistical component of
the uncertainty. The Z and $\gamma$ analyses are described in [9], the
multijet analysis in [10], and the combination in [11]. Jet energy resolution
is around 1% for $p_{T}>100$ GeV and is consistent over all jet energy scale.
Constant term is within expected 3%.

Figure 21: Distribution of the jet energy resolution as a function of the
$p_{T}$ of the jet.

##  IV Conclusion

The ATLAS Tile Calorimeter is an important sub-detector of the ATLAS detector
at the LHC. It is the hadronic sampling calorimeter made of steel plates which
act as absorber and scintillating tiles as active medium. Control of its
energy is essential to measure the energy of jets, hadronically decaying tau
leptons and missing transverse energy. The TileCal calibration system consists
of Cesium radioactive sources, laser, charge injection components, and an
integrator based readout system. Combined information from all systems allows
for an efficient monitoring and correction of fine instabilities of TileCal
cells response. Intercalibration and uniformity are monitored with isolated
charged hadrons and cosmic muons. Data quality in physics runs is monitored
extensively and continuously. All problems are reported and addressed. The
data quality efficiency achieved was 99.6% in 2012, 100% in 2015, 98.9% in
2016 and 99.4% in 2017. The stability of the absolute energy scale at the cell
level was maintained to better than 1% during LHC data-taking. Following the
experience gained during LHC Run-I, all calibration systems were improved for
Run-II. TileCal performance during LHC Run-II, (2015-2018), including
calibration, stability, absolute energy scale, uniformity and time resolution
show that the TileCal performance is within the design requirements and has
given essential contribution to reconstructed objects and physics results.

## References

  * [1] ATLAS Collaboration, “The ATLAS Experiment at the CERN Large Hadron Collider,” JINST 3, S08003 (2008).
  * [2] ATLAS Collaboration, “Readiness of the ATLAS Tile Calorimeter for LHC collisions,” Eur. Phys. J. C 70, 1193 (2010).
  * [3] ATLAS Collaboration, “Operation and performance of the ATLAS Tile Calorimeter in Run 1,” arXiv:1806.02129 [hep-ex].
  * [4] ATLAS Tile Calorimeter system Collaboration, “Design of the front-end analog electronics for the ATLAS tile calorimeter,” Nucl. Instrum. Meth. A 551, 469 (2005).
  * [5] ATLAS Collaboration “Cesium monitoring system for ATLAS Tile Hadron Calorimeter,” Nucl. Instrum. Meth. A 494, 381 (2002).
  * [6] ATLAS Tile Calorimeter Subsystem, Approved Plots https://twiki.cern.ch/twiki/bin/view/AtlasPublic/ApprovedPlotsTile
  * [7] ATLAS Tile Calorimeter system Collaboration, “The Laser calibration of the Atlas Tile Calorimeter during the LHC run 1,” JINST 11, no. 10, T10005 (2016)
  * [8] S. Fracchia, ”Irradiation effect on the response of the scintillators in the ATLAS Tile Calorimeter”, ATL-TILECAL-PROC-2014-013
  * [9] ATLAS Collaboration. ”Determination of the jet energy scale and resolution at ATLAS using ${Z/\gamma}$–jet events in data at $\sqrt{s}=8$ TeV”, ATLAS-CONF-2015-057
  * [10] ATLAS Collaboration, ”Data-driven determination of the energy scale and resolution of jets reconstructed in the ATLAS calorimeters using dijet and multijet events at $\sqrt{s\ }=8~{}TeV$”, ATLAS-CONF-2015-017
  * [11] ATLAS Collaboration, ”Monte Carlo Calibration and Combination of In-situ Measurements of Jet Energy Scale, Jet Energy Resolution and Jet Mass in ATLAS”, ATLAS-CONF-2015-037