11institutetext: Instituto de Astrofísica de Canarias, E-38200 La Laguna,
Tenerife, Spain 22institutetext: Departamento de Astrofísica, Universidad de
La Laguna, E-38205 La Laguna, Tenerife, Spain 33institutetext: INAF –
Astronomical Observatory of Abruzzo, via M. Maggini, sn, 64100 Teramo, Italy
44institutetext: INFN, Sezione di Pisa, Largo Pontecorvo 3, 56127 Pisa, Italy
55institutetext: Universidad de Granada, Departamento de Física Teórica y del
Cosmos, Campus Fuente Nueva, Edificio Mecenas, E-18071, Granada, Spain.
66institutetext: Instituto Carlos I de Física Teórica y computacional,
Universidad de Granada, E-18071 Granada, Spain 77institutetext: Isaac Newton
Group of Telescopes 88institutetext: Institute of Astronomy, Madingley road,
Cambridge, CB3 0HA, UK 99institutetext: Leibniz-Institut für Astrophysik
Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany 1010institutetext:
Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700
AV, Groningen, The Netherlands 1111institutetext: Université Côte d’Azur,
Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Bd de
l’Observatoire, CS 34229, 06304 Nice cedex 4, France 1212institutetext: INAF -
Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti
93/3, I-40129 Bologna, Italy

# Chronology of our Galaxy from Gaia Colour-Magnitude Diagram-fitting
(ChronoGal). I. The formation and evolution of the thin disk from the Gaia
Catalogue of Nearby Stars

Carme Gallart 11 2 2 Francisco Surot 11 Santi Cassisi†, 3344 Emma Fernández-
Alvar†, 11 2 2 David Mirabal†, 2 2 Alicia Rivero †, 22 Tomás Ruiz-Lara†, 55 6
6 Judith Santos-Torres†, 22 7 7 Guillem Aznar-Menargues 22 Giuseppina
Battaglia 11 2 2 Anna B. Queiroz 11 2 2 Matteo Monelli 1122 EugeneVasiliev 88
Cristina Chiappini 99 Amina Helmi 1010 Vanessa Hill 1111 Davide Massari 1212
Guillaume F. Thomas 1122

###### Abstract

Context. The study of the Milky Way is living a golden era thanks to the
enormous high-quality datasets delivered by Gaia, and space asteroseismic and
ground-based spectroscopic surveys. However, the current major challenge to
reconstruct the chronology of the Milky Way is the difficulty to derive
precise stellar ages for large samples of stars. The CMD-fitting technique
offers an alternative to individual age determinations to derive the star
formation history (SFH) of complex stellar populations.

Aims. We aim to obtain a detailed ”dynamically evolved” SFH (deSFH) of the
solar neighbourhood, and the age and metallicity distributions that result
from it. We define deSFH as the amount of mass transformed into stars, as a
function of time and metallicity, in order to account for the population of
stars contained in a particular volume of the MW.

Methods. We present a new package to derive SFHs from CMD-fitting tailored to
work with Gaia data, called CMDft.Gaia and we use it to analyse the CMD of the
Gaia Catalogue of Nearby Stars (GCNS), which contains a complete census of the
(mostly thin disk) stars currently within 100 pc of the Sun.

Results. We present an unprecedented detailed view of the evolution of the
Milky Way disk at the solar radius. The bulk of star formation started between
11-10.5 Gyr ago at metallicity around solar and continued with a slightly
decreasing metallicity trend until 6 Gyr ago. Between 6 and 4 Gyr ago, a
notable break in the age-metallicity distribution is observed, with three
stellar populations with distinct metallicities (sub-solar, solar, and super-
solar), possibly indicating some dramatic event in the life of our Galaxy.
Star formation then resumed 4 Gyr ago with a somewhat bursty behaviour,
metallicity near solar and average star formation rate higher than in the
period before 6 Gyr ago. The derived metallicity distribution closely matches
precise spectroscopic data, which also show stellar populations deviating from
solar metallicity. Interestingly, our results reveal the presence of
intermediate-age populations exhibiting both a metallicity typical of the
thick disk, approximately ${\rm[M/H]}\simeq-0.5$, and supersolar metallicity.

Conclusions. The many tests performed indicate that, with high precision
photometric and distance data such as that provided by Gaia, CMDft.Gaia is
able to achieve a precision of $\lesssim$ 10 % and an accuracy better than 6%
in the dating of stellar populations, even at old ages. The comparison with
independent spectroscopic metallicity information shows that metallicity
distributions are also determined with high precision, without imposing any
a-priory metallicity information in the fitting process. This opens the door
to obtaining detailed and robust information on the evolution of the stellar
populations of the Milky Way over cosmic time. As an example we provide in
this paper an unprecedented detailed view of the age and metallicity
distributions of the stars within 100 pc of the Sun.

###### Key Words.:

Galaxy: solar neighbourhood, Galaxy: stellar content; Galaxy:disk, Galaxy:
evolution, Stars: Hertzprung-Russell and C-M diagrams

$\dagger$$\dagger$footnotetext: These authors contributed in an equally
crucial way to this work

##  1 Introduction

The Milky Way (MW) is the galaxy that we can study with the greatest detail,
over its whole history, by characterising its stellar content on a star-by-
star basis. The low mass stars that formed since the first star formation
events in the Universe inform us of the rate of star formation as a function
of time, and how metals build up in successive stellar generations.

The study of the MW is living a golden era. On the one hand, the impressive
data-sets delivered by the ESA mission Gaia (Gaia Collaboration et al., 2016,
2018, 2023) are revolutionising our current view on our Galaxy (Brown, 2021).
On the other hand, several ground based spectroscopic surveys, ongoing or
planned, (e.g. LAMOST, Wang et al. 2020; Liu et al. 2019; RAVE, Steinmetz et
al. 2020b, a; GALAH, Buder et al. 2021; APOGEE, Ahumada et al. 2020; WEAVE,
Jin et al. 2023; 4MOST, Bensby et al. 2019; Chiappini et al. 2019; or DESI
Cooper et al. 2023), are complementing the Gaia mission by obtaining
spectroscopy with higher resolution and/or down to a fainter limiting
magnitude. Finally, asteroseismic missions such as Kepler (Borucki et al.,
2010) and K2 (Howell et al., 2014) have shown the potential of adding seismic
information to derive ages.

As stated in the ”Gaia Red Book”111Gaia report on the Concept and Technology
Study (sometimes referred to as the ”Gaia Red Book”), see
https://www.cosmos.esa.int/documents/29201/297049/report-science.pdf/ ”A
primary scientific goal of the Gaia mission is the determination of the star
formation histories, as described by the temporal evolution of the star
formation rate, and the cumulative numbers of stars formed, of the bulge,
inner disk, Solar neighbourhood, outer disk and halo of the Milky Way”. The
most difficult part to accomplish this goal, once the necessary data is
available, is the determination of stellar ages, since they cannot be directly
measured: they need to be inferred by comparing the observed properties with
the predictions of stellar evolution models.

The most suitable method for determining the age of a given star depends on
its mass and/or evolutionary stage, and thus deriving homogeneous ages for a
broad range of stellar types or for the full age range is virtually impossible
(Soderblom, 2010; Salaris & Cassisi, 2005, for detailed reviews). In the most
widely used method in Galactic archaeology, a set of physical parameters
derived from spectra (and/or photometry), such as effective temperature,
surface gravity and metallicity (or colours and luminosities), are compared to
a set of stellar evolution models, which predict age as a function of these
parameters (Sahlholdt et al., 2019). Isochrone fitting is in practice prone to
large uncertainties, both due to the difficulty of accurately deriving the
needed stellar parameters, and due to the biases introduced by the isochrone
interpolation. Even in the favourable case of stars with well measured
distances and accurate spectroscopic parameters, typical age errors of
individual stars are around 25% (Sanders & Das, 2018; Mints & Hekker, 2018;
Queiroz et al., 2018; Kordopatis et al., 2023), and only estimated lower in
particularly exquisite instances (e.g. Haywood et al., 2013). In spite of
this, the wealth of data from large spectroscopic surveys, and the increasing
availability of distances (initially from Hipparcos and currently from Gaia),
has led to works exploiting advanced statistical methods, in particular based
on Bayesian statistics (Pont & Eyer, 2004; Jørgensen & Lindegren, 2005), to
derive ages for large stellar samples (e.g. Holmberg et al., 2009; Feuillet et
al., 2016; Sanders & Das, 2018; Mints et al., 2019; Frankel et al., 2019;
Sahlholdt et al., 2022; Xiang & Rix, 2022; Queiroz et al., 2023).

Asteroseismology has turned into the big hope to derive precise stellar ages
(Miglio et al., 2017). When combined with spectroscopy, it provides solid
constraints on stellar mass, radius and evolutionary state (see Mathur et al.,
2012; Chaplin et al., 2020), particularly for bright red giants, which allows
to obtain ages for distant stellar samples. A number of stellar catalogues
have already exploited this combination (Martig et al., 2016; Ness et al.,
2016; Anders et al., 2017; Rendle et al., 2019). However, age errors are still
large (see Fig 22 in Pinsonneault et al., 2018). A best case scenario expected
to provide an age precision of 10% is discussed by Miglio et al. (2017) for
long duration observations like those planned with the Plato satellite.

These individual stellar age determinations require detailed and costly
observations, only possible for a tiny fraction of MW stars. This results in
very complicated selection functions. These samples allow inferring
information such as age-metallicity or age-velocity trends, but it is
basically impossible to retrieve from them the much sought holy grail of
galaxy evolution, that is the star formation history (SFH), or to produce
unbiased age and metallicity distributions directly comparable with
predictions from galaxy models.

The Colour-Magnitude Diagram fitting (CMD-fitting) technique offers a highly
complementary way to approach the problem of deriving the SFH of a composite
stellar population. In the case of the MW, because only Gaia CMDs reaching the
old main sequence turnoffs (oMSTO) are required, SFH derivation is possible
for unbiased and huge stellar samples. With current and forthcoming Gaia data
releases it will be possible, using this methodology, to obtain SFHs and
precise age distributions out to distances of several kpc from the Sun, thus
exploring all Galactic components and addressing major questions of Galactic
astronomy. In fact, the Gaia Red Book proposes the CMD-fitting methodology as
the best method to derive the SFH of the MW222(Gaia Red Book, p. 34), after
the early successes of this methodology in providing SFHs of Local Group
galaxies from deep CMDs obtained from ground based or Hubble Space Telescope
(HST) imaging.

Indeed, in extra-galactic archaeology, deep CMDs reaching the oMSTO with good
photometric accuracy and precision, analysed with the robust technique of CMD-
fitting (Tosi et al., 1991; Bertelli et al., 1992; Tolstoy & Saha, 1996;
Aparicio et al., 1997; Dolphin, 1997; Gallart et al., 1999; Dolphin, 2002;
Aparicio & Gallart, 2004; Gallart et al., 2005; Aparicio & Hidalgo, 2009;
Cignoni & Tosi, 2010), are regarded as the most direct and reliable observable
to obtain a detailed SFH and stellar age distributions of a galaxy, from its
earliest epochs to the present time. CMD-fitting has been, over 25 years, the
standard to determine detailed SFHs for Local Group galaxies, from nearby MW
satellites using ground based data, to more distant members using large
allocations of HST time. This has provided insight on a number of important
topics of near-field cosmology, such as the role of reionization on the early
SFH of dwarf galaxies (Cole et al., 2007; Monelli et al., 2010b, a; Weisz et
al., 2014; Ruiz-Lara et al., 2018), the origin of their different
morphological types (Gallart et al., 2015), the differences between the MW and
M31 satellite systems (Monelli et al., 2016; Skillman et al., 2017), the
spatial gradients (Noël et al., 2009; Cignoni et al., 2013; Meschin et al.,
2014; Rubele et al., 2018) and synchronised SFHs of the Magellanic Clouds
(Ruiz-Lara et al., 2020b; Massana et al., 2022) and the SFHs of the M31 halo,
spheroid, and outer disk (Brown et al., 2008; Bernard et al., 2012). CMD-
fitting has flourished in the context of the study of Local Group galaxies,
because they are sufficiently close to resolve their individual stars, yet far
enough that all their stars can be considered to be at the same distance,
which can be obtained accurately using various well calibrated distance
indicators (Benedict et al., 2007; Beaton et al., 2016). This is fundamental
to transform the measured apparent luminosities to absolute magnitudes that
can be compared with the predictions of stellar evolution models.

In the case of the MW, precise and accurate distances for each individual star
are necessary, and thus the first examples of CMD-fitting to derive the SFH of
stellar samples in the very close solar vicinity (within approximately 50 pc)
used Hipparcos data (Hernandez et al., 2000; Bertelli & Nasi, 2001; Cignoni et
al., 2006). However, the real breakthrough came with the availability of Gaia
data, which allowed to extend the study of the solar neighbourhood to 100-250
pc, first with Gaia DR1 (Bernard, 2018) and then with Gaia DR2 (Alzate et al.,
2021; Dal Tio et al., 2021). Additionally, samples within 2 kpc allowed to
date the accretion time of Gaia-Enceladus and the early SFH of MW thick disk
and halo (Gallart et al., 2019a), and the possible repeated influence of the
Sagittarius dSph on the SFH of the MW disk, since its accretion some 6 Gyr ago
(Ruiz-Lara et al., 2020a), as well as SFH gradients as a function of distance
from the MW plane (Gallart et al., 2019b; Mazzi et al., 2023).

In this paper, we describe in detail our current implementation of the CMD-
fitting technique to derive detailed SFHs from Gaia CMDs, which has been
upgraded in several aspects (e.g. the adopted stellar evolution library, and
the CMD-fitting code itself) with respect to the procedures used in Gallart et
al. (2019a) and Ruiz-Lara et al. (2020a). One salient characteristic of our
methodology is that no a priori assumptions are made on the age-metallicity
relation, on the metallicity distribution, or on the functional form of the
star formation rate as a function of time, SFR(t). We apply this methodology
to derive a first detailed SFH of the Solar neighbourhood, using the exquisite
Gaia Catalogue of Nearby Stars (GCNS) dataset (Gaia Collaboration, Smart et
al., 2021), therefore presenting an unprecedented detailed view of the
evolution of the Milky Way disk at the solar radius.

This paper is organised as follows: Section 2 summarises the contents of the
original GCNS dataset relevant for this study, which are complemented with
Gaia DR3 data on chemical abundances and radial velocities; Section 3 presents
CMDft.Gaia, a new suite of procedures for CMD-fitting specially tailored for
the analysis of Gaia CMDs; Section 4 describes the particular application of
CMDft.Gaia to the CMD of the GCNS while Section 5 presents the derived deSFH
and age-metallicity distributions and the robustness of these results, and
compared them with literature spectroscopic metallicity distributions and age-
metallicity relations; Section 6 discusses the evolutionary history of the
Milky Way (thin) disk in the light of the derived SFH; finally Section 7
summarised the main results and conclusions, both regarding the evolution of
the Milky Way disk and the performance of CMDft.Gaia. A number of Appendix
present complementary information on various aspects of this work.

##  2 The data

We base our analysis on the Gaia Catalogue of Nearby Stars (GCNS, Gaia
Collaboration, Smart et al., 2021), which comprises 331312 stars residing
within a sphere of 100 pc radius centred on the Sun, selected from the full
Gaia EDR3 catalogue. It is a volume-complete sample of all objects with
spectral type earlier than M8 down to the nominal $G=20.7$ magnitude limit of
Gaia.

For details on this catalogue we refer the reader to the original paper (Gaia
Collaboration, Smart et al., 2021). Here suffices to say that it originates
from a selection of all sources in Gaia EDR3 with measured parallaxes
$\omega>8\,\rm{mas}$ (corresponding to a maximum distance of 125 pc), from
which objects with spurious astrometric solutions were removed with a random
forest classifier (Breiman, 2001). Subsequently, posterior probability
densities for the true distance of each source were inferred with a simple
prior not dependent on the sky position or type of star, based on the distance
distribution of stars in GeDR3mock (Rybizki et al., 2020) with observed
parallax greater that 8 mas333For this sample of very nearby, well measured
stars, the dist50 distances obtained in this way are practically identical to
those that would be obtained by simply inverting the parallax. However, for
consistency we adopted dist50 from the GCNS as the distance, and
($\rm{dist_{84}-dist_{16})/2}$ as its error.. Finally, all the stars with a
non-zero probability of being within 100 pc were retained in the catalogue.

In terms of photometric information, apart from Gaia eDR3 data, the original
GCNS also included magnitudes from external optical and infrared catalogues.
We checked the resulting CMDs in a number of combinations of the available
filters and concluded that, for the purposes of SFH derivation, the CMD in the
Gaia bands [$G_{BP}-G_{RP},G$] was definitely superior. This is therefore the
combination we will use in the rest of the article.

Gaia Collaboration, Smart et al. (2021) do not discuss the amount of
extinction affecting the stars in the sample. We have used the Green et al.
(2019) 3D dust map to derive extinctions at the 3D position of each star and
have verified that they are negligible.

In this section we provide a bird’s-eye view of the GCNS stellar content in
terms of the CMD, kinematic and global chemical information, and how it is
globally split between the Milky Way thin and thick disc and stellar halo
components. Given that Gaia DR3 data have become available since the
publication of the GCNS, here we complement the original data-set with Gaia
DR3 line-of-sight velocities for 174221 stars, as well as metallicity ([M/H])
and [$\alpha$/Fe] abundance measurements for 23629 stars from the $Gaia$
Radial Velocity Spectrometer (RVS) as measured by the General Stellar
Parametriser-spectroscopy module, GSP-Spec (Recio-Blanco et al., 2023). We
refer the reader to Appendix A for details on how these quantities were
assembled and on how we assign probabilities of membership for a star to
belong to the thin disc, thick disc or stellar halo.

Figure 1: Summary of the content of the updated GCNS. Upper left panel: Toomre
diagram of the stars with line-of-sight velocities. The thin disk stellar
distribution is represented with violet contours and the thick disk and halo
stars as orange and black dots, respectively. The grey points in the
background represent the whole GCNS sample with kinematic data from DR3.
Middle upper panel: [Fe/H] distribution for the global (grey), thin disk
(purple), thick disk (orange) and halo (black; multiplied by 100). Upper right
panel: [$\alpha$/Fe] distribution of halo stars (black stars), thick disk
stars (orange dots and contours) and thin disk stars (purple dots and
contours). In these two panels, only the stars with abundance information have
been represented. Lower left panel: CMD of the three kinematically selected
components within 100 pc of the Sun, superimposed to the whole original GCNS
sample (in grey): thin disk (purple), thick disk (orange) and halo (black).
Middle panel: CMD of the kinematic thin disk with superimposed isochrones of
solar metallicity with a range of ages (0.03, 0.2, 0.5, 1, 2, 5, 10 and 13
Gyr) from the BaSTI-IAC (solid lines) and PARSEC (Bressan et al., 2012, dashed
lines) stellar evolution libraries. A [M/H]$=-0.5$, 5 Gyr old metallicity
BaSTI isochrone is also plot. Lower right panel: CMD of the kinematic thick
disk with superimposed isochrones from the BaSTI-IAC library: [Fe/H]=0.06 and
5, 10 and 13 Gyr; [Fe/H]=$-0.5,-1.0,-2.0$ and 10 Gyr.

Figure 1 summarises the content of this updated GCNS. In the upper panels and
in the bottom left, the thin disc stars are indicated by purple symbols or
contours, the thick disc in orange and the halo in black. The whole kinematic,
chemical and photometric samples are shown in grey in the upper left, upper
middle and lower left panels, respectively. The upper left panel represents
the Toomre diagram; note that the different components overlap slightly in
this space. The stars classified as halo mostly have retrograde orbits and
follow the distribution that has been associated to the remnants of Gaia-
Enceladus-Sausage (Belokurov et al., 2018; Helmi et al., 2018). The middle
panel presents the metallicity distribution [Fe/H] of each component and that
of the whole sample (that of the halo has been multiplied by 100 to make it
visible). Note that the thick disk distribution is only slightly shifted to
lower metallicity values compared to the thin disk distribution (mean
[M/H]${}_{thin}=-0.11$ dex, $\sigma_{thin}=0.2$ dex; mean
[M/H]${}_{thick}=-0.29$ dex, $\sigma_{thick}=0.4$ dex). Some contamination
from the thin disk may be responsible of this slightly higher metallicity
compared to the ’canonical’ metallicity of the thick disk. Finally, the upper
right panel displays the [$\alpha$/Fe] versus [Fe/H] distribution. Both disk
components, kinematically selected, have a similar range of [$\alpha$/Fe] with
thick disk stars somewhat more extended to high [$\alpha$/Fe] values and lower
[Fe/H]. The reason of the relatively similar distribution of thin and thick
disk stars in this plane is twofold: first, $\alpha$ abundances from Gaia GSP-
Spec come primarily from Calcium measurements, and this element is produced
both by SNIa and SNII, thus not resulting in such a clear-cut separation for
different populations (this is confirmed in the CNN analysis by Guiglion et
al., 2024, where the break in $\alpha$/Fe becomes more evident after combining
Gaia data with APOGEE); second, a separation of thin and thick disk stars with
kinematic (or geometric) criteria, does not necessarily reproduce a chemical
separation (Kawata & Chiappini, 2016). The kinematically selected halo stars
do have a distinct distribution, with all of them having [$\alpha$/Fe]
$\gtrsim$ 0.2 and [Fe/H]¡-0.5.

The lower left panel shows the CMD of the global GCNS catalogue and the three
kinematically selected components. One can appreciate that the dominant thin
disk component (146108 stars) reaches very bright absolute magnitudes and blue
colours in the main sequence, indicative of a very young population; the thick
disk (13153 stars) is clearly much older, and its low main sequence is located
in the blue part of that of the thin disk, indicating a lower metallicity.
Finally, the few halo stars (415 stars, in black) are all located in the blue
ridge of the other two populations, reflecting their even lower metallicity.

The lower middle and right panels of Fig.1 display the CMD of the kinematic
thin and thick disk, respectively, with isochrones superimposed. In the case
of the thin disk, solar metallicity isochrones ([Fe/H]=0.06)444Note that the
chosen solar metallicity [Fe/H]=+0.06 is the initial metallicity of the Sun so
that at its present age of $\simeq 4.5$ Gyr and as a consequence of atomic
diffusion, its measured photospheric metallicity is [Fe/H]=0. with a range of
ages (0.03, 0.2, 0.5, 1, 2, 5, 10 and 13 Gyr) from the BaSTI-IAC555The BaSTI-
IAC isochrones in this figure differ from the official release in the colours
of the faint main sequence portion of the isochrone, below $M_{G}\simeq 8$.
The reason for this corrections and how it is computed is discussed in
Appendix B. (Hidalgo et al., 2018) and PARSEC (Bressan et al., 2012) stellar
evolution libraries have been selected. The stars in this CMD are well matched
by solar metallicity isochrones within the selected age range, from the
youngest to the oldest. The youngest 30 Myr isochrone, while not obviously
needed to match very bright young stars in the main sequence, seems to match
very well the less populated red part of the broad low main sequence (fainter
than MG=8). The blue side of this low main sequence is well matched by a lower
metallicity population (see 5 Gyr isochrone with [Fe/H]= -0.5, in red).
Finally note that there are basically no stars in the subgiant branch between
the 10 and the 13 Gyr old isochrones, hinting to an scarcity of a very old
population in this sample.

The comparison between the BaSTI-IAC and PARSEC isochrones in this panel
reveals some slight differences between the two stellar evolution libraries:
the PARSEC isochrones are systematically redder in the main sequence (see
inset) and red-giant branch, while they are brighter in the sub-giant branch
compared to the BaSTI-IAC isochrones (we refer to Hidalgo et al., 2018, for a
more detailed comparison between the two independent isochrone libraries).
These differences between isochrones from different stellar evolution
libraries indicate that also systematic differences may be expected between
models and data. In the derivation of the SFH, we will quantify this effect by
allowing small shifts between the whole observed CMD and the synthetic CMD to
which it will be compared (see section 3.3.3). The differences between the
predictions of stellar evolution libraries lead to somewhat different SFHs and
age and metallicity distributions (see Section 5.2).

In the case of the thick disk, for clarity only isochrones from the BaSTI-IAC
library, in a range of metallicities, are displayed ([Fe/H]=0.06 and 5, 10 and
13 Gyr; [Fe/H]$=-0.5$, $-1.0$, $-2.0$ and 13 Gyr). The thick disk population
is basically matched by the 10 Gyr old solar metallicity isochrone, with a
minority of stars scattered around the 5 Gyr old isochrone, and very few stars
fainter (thus older) than 10 Gyr. This is similar to what was found by Miglio
et al. (2021) using isochrone ages for a Kepler sample observed by APOGEE.
These authors also found an almost coeval (age scatter of around 1 Gyr) and
old thick disk. Old ages for the thick disk were also found by Queiroz et al.
(2023) using data from different spectroscopic surveys together with Gaia DR3.
The old, lower metallicity isochrones show that there is room for a minority
lower metallicity population which, at faint magnitudes in the main sequence
(M${}_{G}<8$) becomes bluer than the bulk of the population.

This was a first qualitative assessment of the age and metallicity ranges of
the Milky Way components present in the volume covered by the GCNS. In Section
5, we will present a quantitative description of these stellar populations.

##  3 Derivation of SFHs: CMD fitting methodology

In this series of papers, we will define the SFH as the amount of mass that
has transformed into stars, as a function of time (t) and metallicity (Z), in
order to account for the population of stars contained in a particular volume
of the MW. Because the defined volumes will typically be small compared to the
MW size, and stars are expected to move away from their birth position (owing
to diffusive or dynamical processes induced by the spiral arms or the bar, see
e.g. Minchev & Famaey, 2010; Halle et al., 2015; Hayden et al., 2018; Feltzing
et al., 2020), these SFHs may be more appropriately referred to as dynamically
evolved SFHs (deSFH), and the existence of stellar migration will need to be
taken into account in the interpretation of the results. In any case, from
these deSFHs, local age and metallicity distributions can be derived and this
is a fundamentally important information in Galactic Archaeology.

The SFH can be expressed as a combination of simple stellar populations (SSPs)
with a small range of age and metallicity. A convenient way to obtain these
SSPs is from a synthetic CMD computed assuming that stars are born with a
constant probability for all ages and metallicities within a given age and Z
range and adopting other stellar population characteristics such as an initial
mass function (IMF) and a parameterisation of the binary star population. The
errors affecting the absolute colours and magnitudes of the stars, as well as
the completeness function across the CMD need to be simulated in the synthetic
CMD to make it comparable to the observed CMD. We will call this specific type
of synthetic CMD from which we derive the SSPs, mother CMD.

Once these SSPs have been defined, any arbitrary SFH (and its associated model
CMD) can be obtained as a linear combination of SSPs:

$\Psi(t,Z)=\sum_{i}\alpha_{i}\Psi_{i}$

where $\Psi_{i}$ refers to SSP $i$ and $\alpha_{i}$ is the strength attributed
to that SSP for that arbitrary SFH.

The best fit SFH (for a given mother CMD) is then obtained by comparing the
distribution of stars across the observed CMD and in an arbitrary number of
model CMDs, constructed from different combinations of SSPs defined through
sets of $\alpha$, until the best possible match is found.

Figure 2: Summary of the steps involved in the SFH derivation using CMDft.Gaia. A detailed discussion of ChronoSynth and $Dir$SFH is provided in this paper, while Dispar-$Gaia$ was presented in Ruiz-Lara et al. (2022) and will be discussed in more detail in Fernández-Alvar et al. 2024, in prep Table 1: Glossary of acronyms and specific terms Term |  Definition/description
---|---
Age (Z) bins |  Difference in age (Z) between consecutive age (Z) seed points
Age (Z) seed points |  Array of age (Z) values that define a typical ”size” of the SSPs reflecting the varying age (Z) resolution across the whole age (Z) interval
Bundle |  Regions of the CMD containing the stars that are used for the CMD-fitting. The bundle(s) are divided in colour-magnitude boxes (or ’pixels’) where the stars (real and synthetic) are counted
deSFH |  Mass, per unit time and metallicity, that has been transformed into stars somewhere in the galaxy to account for the stars that are today in the studied volume
deSFR(t) |  Marginalisation over metallicity of the deSFH. It gives the mass per unit time that has been transformed into stars somewhere in the galaxy to account for the stars that are today in the studied volume
MDFM |  Metallicity distribution function of the mass transformed into stars
MDFS |  Metallicity distribution function of the stars currently populating the studied volume
Model CMD |  CMD resulting from an arbitrary linear combination of SSPs
Mother CMD |  Specific type of synthetic CMD computed assuming that stars are born with a constant probability for all ages and metallicities within a given age and metallicity range and adopting other stellar population characteristics such as an initial mass function (IMF) and a parameterisation of the binary star population, and with observational errors and completeness simulated in it. The SSPs are derived from a mother CMD
Solution CMD |  CMD obtained with the best fit SFH, and matching the number of stars of the observed CMD inside the bundle
SSP |  Simple Stellar Population, population of synthetic stars in a small range of age and metallicity
Synthetic CMD |  Any CMD obtained from the predictions of stellar evolution models

In the following we will discuss in detail all the steps involved in our
implementation of the SFH derivation procedure, which we call CMDft.Gaia
(first introduced in Ruiz-Lara et al. 2022). It has been considerably updated
compared to previous works such as Gallart et al. (2019a) and Ruiz-Lara et al.
(2020a).

CMDft.Gaia is a suite of procedures specifically designed to deal with Gaia
data that includes i) the computation of synthetic CMDs in the Gaia bands (Sec
3.1; ii) the simulation in the synthetic CMDs of the observational errors and
completeness affecting the observed CMD after quality and reddening cuts
(Section 3.2; and iii) the derivation of the SFH itself (Section 3.3). We
explain these steps in detail below, while a summary is presented in Figure 2.

###  3.1 Synthetic CMD computation

All mother CMDs adopted in the present investigation have been computed with
our own synthetic CMD code, that results from a deep evolution/update of the
code presented in Pietrinferni et al. (2004) and Cordier et al. (2007) and
that we will call ChronoSynth. Due to the need of computing mother CMDs
hosting a huge number of synthetic stars, the current version of the code has
been parallelised in order to speed up the whole computational procedure. The
code provides magnitudes and colours of stars belonging to a synthetic stellar
population with an arbitrary SFH, as well as the total mass that has been
transformed into stars associated to that population, which includes that of
the stars already dead and thus not present in the synthetic CMD. For this
purpose, the code relies on a grid of isochrones in a wide age and metallicity
range, which depends on the adopted stellar model library.

The code has some flexibility in the types of SFH that can be adopted.
However, for the purpose of SFH derivation, the synthetic mother CMDs are
computed adopting a flat probability distribution for star to be born within
the whole defined age and metallicity range (the latter can be defined flat in
Z or log(Z)). The lower and upper limits both in age and metallicity have to
be specified in the input file, as well as the adopted IMF666The current
version of the code allows to choose between a power law (with exponent chosen
by the user), and the Kroupa et al. (1993) IMF. and the characteristics of the
binary population, which are parameterised as a function of the binary
frequency $\beta$ and minimum mass ratio qmin.

To create the synthetic stellar population, with a desired number of stars
(Nstars) down to a given limiting magnitude (Mmax), according to the adopted
SFH, a random value of the stellar age and metallicity are drawn from the
whole age/metallicity range. Then a stellar initial mass, M, is randomly
selected following the adopted IMF. These values of age, stellar mass, and
metallicity are used to interpolate in the isochrones of the selected grid to
determine the bolometric luminosity, effective temperature, and current value
of the mass777The current mass of the star, that can be different than the
initial one as a consequence of mass loss, is needed in order to properly
evaluate the actual surface gravity that is – together with the effective
temperature and chemical composition – the quantity needed to evaluate the
bolometric correction for any selected photometric passband. The occurrence of
mass loss is taken into account by adopting a stellar model grid computed by
accounting for a selected value of the mass loss efficiency (we refer to
Hidalgo et al., 2018, for a more detailed discussion on this topic). of the
synthetic star. These properties are then used to predict its absolute
magnitudes in the various Gaia DR3 photometric passbands 888Many other
photometric systems are also available in the code. on the basis of the
bolometric correction tabulations999These bolometric correction tables are
used for all the libraries of stellar models than can be selected as input in
the synthetic CMD code. This means that they are adopted also when selecting
the PARSEC library; however it is worth noting that both the BaSTI/BaSTI-IAC
library and the PARSEC one adopt very similar spectral libraries in the regime
appropriate for main sequence and sub-giant branch stars as discussed by
Hidalgo et al. (2018) and Chen et al. (2019). provided by Hidalgo et al.
(2018).

To include a given population as unresolved binary systems, the binary
fraction $\beta$ and the minimum mass ratio $q_{\rm min}$ between secondary
and primary stars have to be specified. Then, for each generated synthetic
star, an additional random number is used to determine if it is a component of
an unresolved binary. If this is the case, the mass of the secondary star is
randomly selected from the distribution given by Woo et al. (2003), that is,
the mass of the secondary of a primary star with mass Mp is randomly chosen -
according to a flat distribution - between $q_{\rm min}\times\rm{M_{p}}$ and
$\rm{M_{p}}$. The predicted properties (luminosity, effective temperature,
magnitudes) of this unresolved binary system are calculated properly adding
the fluxes of the two unresolved components. No evolution of the binary system
itself is considered.

Finally, in order to investigate the impact of the choice of different stellar
model libraries on the properties derived for the studied stellar populations,
we have implemented different model grids in ChronoSynth. The current version
of the code allows to use the following libraries: 101010In addition to the
two libraries quoted in the main text, the first release of the BaSTI
isochrone grid (hereinafter BaSTI) for both the solar scaled (Pietrinferni et
al., 2004) and $\alpha-$enhanced (Pietrinferni et al., 2006) mixture is also
kept in the code in order to allow – if needed – some comparison with previous
works. Indeed, this library is the one used in Gallart et al. (2019a) and
Ruiz-Lara et al. (2020a). The selected model set from the former library is
the one accounting for convective core overshooting during the central
H-burning stage, and mass loss according to the Reimers law with the free
parameter $\eta$ fixed at a value equal to 0.4.

-the BaSTI-IAC grid both for solar-scaled (Hidalgo et al., 2018) and $\alpha-$enhanced heavy element mixture (Pietrinferni et al., 2021). The selected grid is that accounting for diffusive processes, core convective overshooting and mass loss (with efficiency fixed by selecting a value for the free parameter $\eta$ equal to 0.3);

-the PARSEC stellar model library for a solar-scaled mixture as provided by Bressan et al. (2012).

###  3.2 Simulation of the observational errors and completeness

Prior to the comparison between the observed and synthetic star distribution
across the CMD, it is necessary to simulate the observational errors and the
completeness in the synthetic CMD. In the case of the GCNS, since it is
basically complete and the photometric and distance errors are really small,
we adopted a simplified error simulation procedure (see Section 4.2). A
comprehensive description of the general method we will adopt in this series
of papers to simulate the observational errors and completeness in the
synthetic CMD, called DisPar-$Gaia$ will be provided in Fernández-Alvar et al.
(2024, in prep), while a summary is provided in Ruiz-Lara et al. (2022) (see
also Figure 2).

###  3.3 Determination of the best fit SFH

To determine the best fit SFH through comparison of the observed CMD with
model CMDs computed from combinations of SSPs, we use a new CMD-fitting
package that we call $Dir$SFH111111Because an important feature of the code is
the use of a $Dir$ichtlet (or Voronoi) tessellation algorithm to define the
SSPs. This package is a sophisticated evolution of previous CMD-fitting
software such as IAC-pop (Aparicio & Hidalgo, 2009) and TheStorm (Bernard et
al., 2018; Rusakov et al., 2021). We will discuss the novel approach followed
by $Dir$SFH to the different steps involved in the CMD-fitting process, which
lead to extremely robust solutions providing a large amount of details in the
age-metallicity plane. We will mention the main differences with respect to
IAC-pop and TheStorm, the latter being the code used in our two first works
delivering SFHs of the MW with Gaia CMDs, namely Gallart et al. (2019a) and
Ruiz-Lara et al. (2020a).

####  3.3.1 Definition of the simple stellar populations (SSPs)

An innovative aspect of $Dir$SFH lays in the way SSPs are defined. The mother
CMD is dissected as a function of age and metallicity in a semi-random
Dirichlet-Voronoi tessellation, based on a grid of seed points in both age and
metallicity, which allows the user to define a typical ”size” of the SSPs
reflecting the varying age and metallicity resolution across the whole age and
metallicity interval. A minimum number of stars is required in each SSP, and
if this number is not reached, a SSP may be merged with neighbour ones until
it has the minimum required number of stars. Therefore, a large enough mother
diagram is necessary to avoid degrading the resolution in age and metallicity
of the derived SFH.

A large number N of ’individual’ SFHs are derived with different sets of SSPs
constructed with the above method, and the final SFH is the average of all of
them. Figure 3 shows three individual SFH solutions in which the Dirichlet-
Voronoi tessellation of the mother CMD can be appreciated. Note that each
solution is somewhat different from another, even though they share general
trends, including a low metallicity old population with very low signal, as
well as some high metallicity population.

Figure 3: Three individual solutions of the GCNS SFH are shown. The Dirichlet-
Voronoi tessellations of the mother CMD, different for each solution, can be
appreciated, while the overall shape and characteristics of the derived SFH is
preserved.

####  3.3.2 Sampling the star’s distribution across the CMD and weights of the
fit

As mentioned in the introduction of this section, the distribution of stars in
the observed CMD and in the model CMDs resulting from each combination of
SSPs, need to be compared. For this, observed and model CMDs are binned as a
function of colour and magnitude in the same way.

The amount of information on the SFH that a particular region of the CMD can
provide depends mainly on a) how separated in colour and magnitude the
populations of different ages and metallicities are in that region, b) the
accuracy with which stellar evolution models are able to predict the positions
and lifetimes of the stars (Gallart et al., 2005) in the corresponding
evolutionary phase, and c) the number of stars populating it. The main
sequence (and particularly the region around the turnoffs) and the sub-giant
branch are thus the CMD regions that provide most information on the SFH. In
contrast, the (shorter lived) red giant branch, the red clump and the
horizontal branch phases are affected both by larger uncertainties in the
stellar evolution model predictions (including more uncertain bolometric
corrections) and by the superposition in a small colour-magnitude region of
stars that encompass almost the whole range of possible ages and
metallicities. The latter leads to a poor age and metallicity resolution which
is exacerbated by a larger age-metallicity degeneracy compared to the rest of
the CMD. To take this into account, in IAC-Pop or TheStorm, several bundles in
the CMD were typically defined, each in turn divided in colour-magnitude boxes
(or ’pixels’) of different sizes where the number of stars are counted for the
comparison between the observed and the model CMD. These bundles could
exclude, or sample more coarsely, CMD regions corresponding to certain stellar
evolutionary phases, in order to modify their overall weight in the fit. This
approach has the disadvantage of a certain subjectivity in the bundle
definition, which however, was shown to have little effect on the final SFH
(see, for example Ruiz-Lara et al., 2021). $Dir$SFH uses a single bundle
defined by the user to tightly include both the observed and the mother CMD
down to a given limiting magnitude, which is used only as a delimiter of the
fitting space. Within this bundle, a weight matrix is calculated from the
mother CMD based on how precisely a given ’pixel’ of that CMD is unique in
terms of age (see next paragraph for a discussion of how these ’pixels’ are
defined). In particular, in the current implementation, the weight of each
’pixel’ in the CMD is calculated as the inverse of the variance of the stellar
ages in that ’pixel’, such that a ’pixel’ populated by a small range of age
will be given more weight (see Figure 4 for an example of the weights applied
across the CMD for the q01b03_120M_M5121212From now on, this will be the
typical naming convention for mother CMDs. In this particular case,
q01b03_120M_M5 means a CMD computed adopting qmin=0.1, $\beta$=0.3, and with
120 million stars down to $M_{G}$=5.0. The Basti-IAC stellar evolution library
is used unless stated otherwise mother CMD). This occurs, for example, in the
bright, blue parts of the CMD, which are populated exclusively by young stars.
Giving more weight to those areas defines the populations therein (mostly
young SSPs) with higher accuracy, imposing a strong constrain on the presence
of young stars in areas where the CMD has severe overlap of ages. The option
of a uniform weight across the whole CMD is also possible, and we have
verified that the results are not affected in a significant way by the
weighting scheme applied (see Figure 27).

Figure 4: Example of the weights applied across the CMD for the q01b03_120M_M5
mother CMD. The weight of each ’pixel’ in the CMD is calculated as the inverse
of the variance of the stellar ages in that ’pixel’. The bundle defining the
area where observed and model CMDs are compared is also shown. Note that the
bundle tightly delimits the mother CMD; a few ’pixels’ containing synthetic
stars outside the bundle correspond to stars that have been scattered in the
simulation of the observational errors.

The way the CMD ’pixels’ are defined in $Dir$SFH takes into account the fact
that the number of stars in each pixel is subject to statistical fluctuations.
In TheStorm this was taken into account by shifting the limits of the boxes by
a fraction of their size and calculating new SFHs, which would at the end be
averaged together and used to calculate statistical errors. In $Dir$SFH, a
’coarse’ grid of pixels is defined (typically of size 0.2 magnitudes in MG and
0.1 magnitudes in $(G_{BP}-G_{RP})$) and used to calculate a colour-magnitude
histogram. Then, the grid is shifted a number of times (25 times, for 5x5
steps in colour and magnitude in the current implementation) and new
histograms are calculated for each shift. This results in an over-sampled
representation of the CMDs, maintaining the number statistics of the coarse
grid, but adding the information of the variation of the number of stars
across the CMD in much finer steps (effectively increasing the ’resolution’ by
25 times). This preserves fine details in the stellar distribution across the
CMD that would have otherwise been destroyed by the simple coarse grid.

####  3.3.3 Search for systematic shifts between theory and observations

Uncertainties in the effective temperature scale, and/or in the bolometric
corrections adopted to transfer stellar evolution predictions from the H-R
diagram to the observational plane (in this case the Gaia photometric system),
as well as residual uncertainties in the photometric calibration, may lead to
slight overall systematic shifts between the observed and the synthetic CMDs.
An evidence of the existence of such shifts can be observed in the data shown
in the lower middle panel of Figure 1, where PARSEC and BaSTI-IAC isochrones
have been superimposed to the Gaia CMD in the absolute plane: in most
evolutionary phases, small but noticeable systematic shifts do exist between
isochrones of identical age and metallicity belonging to the two model sets
(as also discussed in Hidalgo et al. (2018)). Similar systematic shifts may be
expected between data and models.

In order to derive the size and direction of these systematic offsets, several
SFHs are calculated with the mother CMD shifted in colour-magnitude space
within a maximum specified range. Then, the residuals of these fits are
analysed, and an appropriate weighted average of the colour and magnitude
shifts leading to the smallest residuals is adopted as the best shift for the
final SFH calculation. This is a similar, but slightly more sophisticated
procedure compared to that adopted in TheStorm or IACpop (see, for example,
Figure 5 in Rusakov et al., 2021).

####  3.3.4 The minimisation algorithm

In $Dir$SFH, the goodness of the coefficients in the linear combination that
defines each model CMD is measured through a Skellam distribution (as opposed
to simple Poisson statistics in TheStorm). This statistic considers both the
observed and model CMD histograms to be stochastic in nature. For any given
colour-magnitude pixel, the difference in the number-counts between observed
and model is evaluated, and the probability of this result, considering both
inputs to be Poisson distributed, is used to calculate the goodness-of-fit of
the ensemble. This implies a critical difference with TheStorm approach: in
$Dir$SFH, cases where there are observed stars but no model CMD stars are
treated completely equivalent to the cases where there are no observed CMD
stars but model stars are present. In particular, $Dir$SFH minimises:

$S=\sum_{i}{\frac{(O_{i}-M_{i})^{2}}{O_{i}+M_{i}}}$

Where $O_{i}$ and $M_{i}$ are the number counts for each pixel in the ensemble
for the observed CMD and the model CMD, respectively.

####  3.3.5 Calculation of the final SFH and associated errors

The CMD-fitting procedure is repeated for a large number ($N\simeq 100$) of
different realisations of the SSPs and the final SFH is a weighted average of
the resulting $N$ SFHs, with the error being the dispersion of the
distribution.

####  3.3.6 The outputs of $Dir$SFH

$Dir$SFH produces two main outputs:

i) The SFH, that is, the mass transformed into stars as a function of lookback
time (age) and metallicity [M/H] in units of M☉ Gyr-1 dex-1, in the form of a
800$\times$800 grid of star formation rate values and their uncertainties as a
function of age and metallicity. From this information, the SFR(t) (in units
of M☉ Gyr-1) and the metallicity distribution function of the astrated mass
(MDFM) (in units of M☉ dex-1) are calculated by marginalising over metallicity
and age, respectively. From the SFR(t), the cumulative mass function is
calculated.

ii) A solution CMD, which is obtained by sampling the mother CMD according to
the derived SFH, until the same number of stars in the observed CMD is
obtained. In the solution CMD, each star has information of its age and
metallicity, in addition to its colour and magnitude. The solution CMD,
therefore, allows analysing the current stellar content of a given stellar
population, in terms of the number of stars currently present as a function of
their age and metallicity. In particular, the metallicity distribution
function of the stars in the sample (MDFS hereafter) can be obtained, which is
an important observable that can be compared with that derived
spectroscopically. Note that no information is obtained about the age of
individual stars in the observed CMD.

##  4 Deriving the SFH of the GCNS

In this section we will discuss the particular application of CMDft.Gaia to
the GCNS (Gaia Collaboration, Smart et al., 2021).

###  4.1 Synthetic CMDs used to derive the deSFH within 100pc

As discussed in 3.1, in addition to the adopted stellar evolution library, a
number of choices have to be made to calculate the synthetic CMDs that will be
used to derive the SFH. Along the paper we will explore the impact on the SFH
of the most relevant of these choices. In particular, the extraordinary depth
of the GCNS CMD will allow us to check different assumptions on the binary
star population using the distribution of stars in colour and magnitude in the
low main sequence. For a summary of the synthetic CMDs used in this work, see
Table 13. The choices made to compute this set of synthetic CMDs are the
following:

i) Stellar evolution library. The reference stellar evolution library used is
BaSTI-IAC with the solar scaled mixture (Hidalgo et al., 2018). Several
synthetic CMDs with a number of choices on the other parameters have been
calculated with this library. Two additional synthetic CMD has been computed
using the PARSEC solar-scaled library (Bressan et al., 2012).

ii) Age and metallicity distribution, and number of stars in the synthetic
CMDs: all synthetic CMDs have been computed with a flat age and metallicity
distribution (flat in Z) within the age and [M/H] range indicated in Table 13,
and with a specified total number of stars down to a given limiting magnitude
MG,max.

For each synthetic CMD, the age and metallicity limits as well as the number
of stars brighter than $M_{G}$=5 (which allows a homogeneous comparison of the
size of each CMD in the approximate portion used to calculate the SFH, above
the oMSTO) are given in Table 13.

iii) Binary star population: we have explored different populations of
unresolved binary systems, parameterized as a function of the fraction of
unresolved binaries $\beta$ and minimum mass ratio $q_{\rm min}$, as discussed
in Section 3.1. The $\beta$ and $q_{\rm min}$ adopted combinations are
specified in Table 13.

iv) IMF: the Kroupa et al. (1993) IMF has been used, except in one case, for
which the Salpeter IMF was adopted.

Table 2: Characteristics of the synthetic mother CMDs used for the analysis of
the SFH of the GCNS.

Synthetic Diagrams
---
Name | Age range | [M/H] | qmin | $\beta$ | Library | MG,max | N* ($\times 10^{6}$) | shift
| (Gyr) | range |  |  |  |  | (M${}_{G}\leq 5$) | (col,mag)
q01b03_60M_MG10 | 13.5-0.02 | -2.2–0.45 | 0.1 | 0.3 | BaSTI-IAC | 10 | 7.0 | (-0.036, 0.041)
q01b03_112M_MG6 | 13.5-0.08 | -2.2–0.45 | 0.1 | 0.3 | BaSTI-IAC | 6 | 68.2 | (-0.036, 0.040)
q01b03_120M_MG5 | 13.5-0.08 | -2.2–0.45 | 0.1 | 0.3 | BaSTI-IAC | 5 | 120 | (-0.035, 0.037)
q01b03_30M_MG10_parsec | 13.5-0.02 | -2.2–0.27 | 0.1 | 0.3 | PARSEC | 10 | 4 | (-0.047,0.045)
q01b03_43M_MG5_parsec | 13.5-0.02 | -2.2–0.27 | 0.1 | 0.3 | PARSEC | 5 | 43.2 | (-0.035,0.030)
q01b05_60M_MG10 | 13.5-0.02 | -2.2–0.45 | 0.1 | 0.5 | BaSTI-IAC | 10 | 6.9 | (-0.030, 0.030)
q01b07_30M_MG11 | 13.5-0.02 | -2.2–0.45 | 0.1 | 0.7 | BaSTI-IAC | 11 | 2.2 | (-0.027, 0.030)
q01b07_81M_MG5 | 13.5-0.02 | -2.2–0.45 | 0.1 | 0.7 | BaSTI-IAC | 5 | 81.3 | (-0.030, 0.030)
q04b03_30M_MG11 | 13.5-0.02 | -2.2–0.45 | 0.4 | 0.3 | BaSTI-IAC | 11 | 2.2 | (-0.034, 0.038)
q06b01_30M_MG10 | 13.5-0.02 | -2.2–0.45 | 0.6 | 0.1 | BaSTI-IAC | 10 | 3.5 | (-0.039, 0.048)

131313Columns, from left to right: identification name, which includes
information on the millions (M) of stars down to a given limiting magnitude
(MG), age and metallicity ranges, the minimum binary mass ratio (qmin), the
binary fraction ($\beta$), the stellar evolution library, the limiting
magnitude MG,max, the number of stars with a M${}_{G}\leq 5$ mag and, the
shift in colour and magnitude subtracted to the synthetic CMD.

###  4.2 Simulation of the GCNS observational errors

In the case of the GCNS, since both the photometric errors and those derived
from the distance calculation are really small and other sources of error
(such as reddening) are negligible, we adopted a simplified error simulation
procedure compared to DisPar-$Gaia$, which will be used in future papers of
this series.

We considered that the sources of uncertainty affecting the position of an
observed star in the CMD were solely the photometric errors and the error in
the determination of the distance. In order to implement these observational
effects on the mother diagrams, we first assigned to each synthetic star a
distance following the global distribution of stellar distances in the GCNS
(dist_50). This preliminary step allowed us to, applying the relation between
apparent and absolute magnitudes, move our mother CMD to the apparent plane
(we verified that extinction is totally negligible in this sample). As shown
in Riello et al. (2021) there is a clear trend of photometric errors with
magnitude, with the brighter tail (G$\lesssim 6$) presenting larger
uncertainties, and a smooth trend to larger uncertainties for fainter stars as
well. We found a similar trend for the distance errors, with more distant
stars displaying larger uncertainties, which we defined as
(dist_84-dist_16)/2.0, being ’dist_84’ and ’dist_16’ the 84th and 16th
percentiles of the distance PDF (1$\sigma$ upper and lower bounds, assuming a
gaussian distribution) provided by Gaia Collaboration, Smart et al. (2021). In
both cases (photometry and distance information), we fitted a 5th degree
polynomial ($P_{i}$) to the run of the error in the parameter $i$ as a
function of the value of such parameter, with $i$ corresponding to apparent
GBP, GRP, G or distance. Also, we characterised the running standard deviation
of each parameter ($\sigma_{i}$). Thus, to each star $s$ in the synthetic CMD,
with a given set of parameters $i$ (G, GBP, GRP, distance), we assigned random
errors ($\epsilon_{i}$) following a Gaussian centred at $P_{i}(i_{s})$ and
with sigma $\sigma_{i}(i_{s})$, with $i_{s}$ referring to the value of the
parameter $i$ for the star $s$.

Once we have the photometric and distance errors for each synthetic star, we
performed a quadratic propagation of uncertainties to translate these
observational errors into an error in the absolute magnitude and colour. These
errors were simulated by adding to the colour and magnitude of each star a
correction following a Gaussian distribution centred at zero and with sigma
equal to the propagated error in colour and magnitude. Finally, since the GCNS
is considered to be complete down to spectral type M8 (Gaia Collaboration,
Smart et al., 2021, much fainter than the stars to be used for the SFH
calculation), we have not applied any completeness simulation in the synthetic
CMD.

###  4.3 Parameterising the mother CMD and configuring $Dir$SFH

As discussed in Section 3.3, $Dir$SFH allows the user to define a number of
input parameters that will determine some details of the fit. In this section
we will discuss how these parameters have been chosen:

####  4.3.1 The arrays of age and metallicity seed points used to define the
SSPs

The decrease in the isochrone separation toward older ages results into a
decrease of the age resolution. Nevertheless, the actual age resolution as a
function of age may also depend on the characteristics of the observed CMD,
and in particular, on the photometric errors across it. The age and
metallicity seed points need to be carefully chosen to optimise the recovery
of the age and metallicity information present in the data, while avoiding
over-fitting the CMD and keeping manageable computing times.

For metallicity, and after some testing, we have adopted a typical separation
of the seed points of 0.1 dex in [M/H], which is of the order of the typical
error in spectroscopic metallicity measurements.

In order to assess the age precision and accuracy that can be achieved with a
high quality Gaia CMD such as that of the GCNS, we have designed recovery
tests based on: i) a composite CMD of four MW open clusters observed by Gaia,
and ii) seven synthetic clusters with ages 0.2, 2, 4, 6, 8, 10 and 12 Gyr, a
small age range (20 Myr) and a small metallicity range, close to solar
metallicity ([M/H]=-0.1 to -0.05), which is the metallicity of most stars in
the observed GCNS CMD. These tests consist on testing $Dir$SFH recovery using
several arrays of age seed points, which result in corresponding arrays of age
bins, as we will refer, for simplicity, to the difference in age between
consecutive age seed points. These tests will allow us to check how the age
precision and accuracy vary with different age bins. In Appendix C we describe
how these datasets have been prepared (membership selection, distance and
reddening determination in the case of the open clusters, and ChronoSynth
input parameters for the calculation of the synthetic clusters).

From a set of age seed points similar to the one used in previous
works141414The set of age seeds that we used as starting point is the
following: Ages(Gyr)= [0.08, 0.195, 0.340, 0.477, 0.606, 0.732, 0.877, 1.057,
1.293, 1.664, 2.190, 2.736, 3.295, 3.852, 4.4585, 5.231, 6.197, 7.235, 8.279,
9.324, 10.368, 11.418, 12.456, 13.5]. It is similar to the set used by Ruiz-
Lara et al. (2020a) but the age bins have been optimised by taking into
account the typical separation of the isochrones in the populated areas near
the turnoff point, as a function of age. (e.g. Ruiz-Lara et al., 2021, 2020a),
which results in what we will call XL age bins, we created three new sets of
age seed points that result in progressively smaller age bins, following the
same functional relation between bin size as a function of age as the original
set. We will call the four sets XL, L, M and S bins. We have derived the SFH
of the composite CMD of the four open clusters and that of the seven synthetic
clusters mentioned above, with the four sets of age bins, keeping unchanged
all the other parameters involved in the fit. The derived SFHs are presented
in Appendix C, while we discuss here the main conclusions regarding the
precision and accuracy of the age determination as a function of age.

For each real or synthetic cluster we have fitted a 2D gaussian151515Using a
program based on the Scikit-learn Python package and the Gaussian Mixture
Model code to the corresponding distribution of ages and metallicities of the
stars of the solution CMD in the age-metallicity plane. The square root of
each diagonal term of the corresponding covariance matrix (that is, the
standard deviations projected in the age and metallicity axis, $\sigma_{age}$
and $\sigma_{Z}$, respectively) provide a measure of the age or metallicity
precision, while the comparison of the 2D fitted gaussian centre with the
input mean age or metallicity of the corresponding synthetic cluster gives
information on the accuracy of the derived ages or metallicities.

The metallicity is recovered accurately at all ages, and with a $\sigma_{Z}$
between 0.05 and 0.10 dex (and thus, of the order of the size of metallicity
bins, 0.1 dex), with little dependency of the number of stars in the cluster
or the size of the age bins (see discussion in Appendix C and Figures 22 and
23).

The left panel of Figure 5 displays $\sigma_{age}$ as a function of the age
bin size for the four open clusters and the seven synthetic clusters. For each
cluster, four points indicating the measured $\sigma_{age}$ for the XL, L, M
and S bins are connected with a line. The age bin sizes in the $x$ axis
correspond to those at the age of each cluster. In the case of the synthetic
clusters, three sequences of $\sigma_{age}$ are displayed, corresponding to
simulations with different number of stars brighter than the faint limit of
the bundle (MG=4.2): $\simeq$ 350 per cluster (similar to the open clusters),
$\simeq$2000, such that the total number of stars in the seven clusters is
similar to that in the GCNS CMD, and $\simeq$ 14000, in order to check whether
a much larger number of stars leads to a substantially greater precision. In
the case of the open clusters, two lines indicate the results of solutions
with two different mother CMDs with different unresolved binary star
characteristics: q01b03_120M_MG5 (dashed line and open symbols) and
q01b07_81M_MG5 (solid line and filled symbols).

Figure 5: Left panel: precision analysis. $\sigma_{age}$ as a function of the
age bin size is displayed for the four open clusters (black symbols) and the
seven synthetic clusters (coloured symbols). For each cluster, four points
indicating the measured $\sigma_{age}$ for the XL, L, M and S bins are
connected with a line. The three sequences of $\sigma_{age}$ for the synthetic
clusters correspond to clusters simulated with different number of stars with
MG¡4.2: $\simeq$ 350, $\simeq$2000 and $\simeq$ 14000\. The two sequences of
$\sigma_{age}$ for the open clusters refer to solutions obtained with two
different mother CMDs with different unresolved binary star characteristics:
q01b03_120M_MG5 (dashed line and open symbols) and q01b07_81M_MG5 (solid line
and filled symbols). Right panel: accuracy analysis. The relative error of the
age determination, as a function of age, is represented. Only the synthetic
clusters are displayed. Clusters with different numbers of stars are depicted
in different colours as indicated in the labels. The size of the bins used for
each measurement is represented with the corresponding letters.

The conclusions that can be extracted from the left panel of Figure 5 are:

  * •

The $\sigma_{age}$ decreases with the size of the age bins. While this could
be expected, the fact that even for the oldest clusters (12-10 Gyr old) the
precision is better for smaller age bins is showing that, within this range of
bin sizes, we are not hitting a physical limit imposed by the decreasing
separation of the isochrones with increasing age, at least for this CMD with
very small photometric and distance errors.

  * •

The precision depends also on age, with a ’break’ between 6-4 Gyr: note that
for a very similar age bin size, the age of the 4-2 Gyr old synthetic clusters
is recovered with greater precision than that of the 6 Gyr old cluster. In
Figure 19 it can be seen that the separation of the synthetic clusters starts
to decrease faster for the synthetic clusters older than 6 Gyr. Also, the
shape of the isochrone changes around this age reflecting the transition from
radiative to convective H-burning cores. For ages older than 6 Gyr, the bin
size as a function of age is basically kept constant for each XL, L, M and S
set. In this range, it can be observed that the age of the 6 Gyr cluster
(green) is recovered with slightly better precision than the other older
clusters, for which the sequences are quite mixed, indicating little
dependence of age precision on age for ages older than 8 Gyr.

  * •

The age precision shows little dependency with the number of stars in the
population. It is remarkable that even with only a few hundred stars, $Dir$SFH
is able to determine the age as precisely as with 40x the number of stars.
This is an important finding as it shows that this methodology can be used
confidently to determine age distributions even for minority populations in
the MW.

  * •

The age of the open clusters is determined as precisely as that of synthetic
clusters of similar age, showing that the BaSTI-IAC models are able to match
the data remarkably well, and that the possible mismatch between the observed
populations and those simulated in the mother CMD (as, for example, the
characteristics of the binary star population or the IMF) do not affect
substantially the recovered age precision. In fact, in the case of the open
clusters, we have recovered the SFH with two mother CMDs with different binary
fraction (30% and 70%). For M67, the age precision is basically identical with
the two binary fractions, while for the other three clusters, a better
precision is reached with 70% binaries. This may reflect different actual
binary fractions in each individual cluster. However, except in the case of
Pleiades with 30% binary fraction, the age precision achieved is similar to
that reached for the synthetic clusters for which there is a perfect match of
the binary fraction between the mock and the mother.

  * •

The gray horizontal lines indicate the locus of 10% precision in age for ages
0.2, 2, 4 and 6 Gyr. The comparison with the sequences of the clusters
indicates that, for ages 6-12 Gyr, even for the XL age bins, the age of the
populations is recovered with precision better than 10%, a goal that is within
the best expectations of Galactic Archaeology, even when including
asteroseismology information (Miglio et al., 2017). 161616Note that, in this
case, we are referring to the age of a population of stars rather than a
single star. For the S bins, the ages of the 2-12 Gyr old clusters are
recovered with a precision of the order or better than 5%, (with better
relative precision for older clusters in each group 2-4 and 6-12 Gyr old),
while for the youngest clusters (both synthetic and open) the best resolution
achieved is 10%.

The right panel of Figure 5 displays the relative error of the age
determination (a measure of the accuracy), as a function of age. In this case,
only the synthetic clusters are displayed, as for them we can compare the
recovered age with the mean input age. Clusters with different numbers of
stars are depicted in different colours as indicated in the labels, while the
age bins used for a particular measurement are represented by the
corresponding letters. The $x$ position of the symbols has been shifted
slightly for clarity. From this figure, we can conclude that ages are
systematically overestimated by a maximum of 6%. The only instances with a
lower accuracy, up to 8%, are in the case of the XL bins for the 10 Gyr
cluster or, in the case with fewer stars, for the youngest cluster. The latter
can be easily understood as an effect of the very small number of stars, which
may result in an undersampled main sequence turnoff mimicking an older age.
For the intermediate-age range (6-2 Gyr), the accuracy is better than 4%. The
figure also shows that neither the size of the bins, nor the number of stars
in the population (except for very young ages) are systematically related to
the accuracy of the age recovery.

In Appendix C we present this study of the age accuracy and precision in more
detail.

####  4.3.2 The area in the CMD included in the fit and whether weights are
provided within it

Ruiz-Lara et al. (2021) showed that different bundle strategies had little
effect on the resulting SFH of the dwarf galaxy Leo I. With $Dir$SFH, we
consider a single bundle including the whole observed and mother CMD and we
have verified that the resulting SFH has little dependence on its exact shape.
As mentioned, we consider that the use of a single bundle removes
subjectivity, improving the repeatability of the results, and maximizes the
information used to compute the SFHs. We have paid special attention to test
whether the faint magnitude limit of the bundle would affect the precision of
the derived SFH. Inspection of the isochrones in Figure 1 (lower, middle
panel) indicates that the best age sensitivity can be expected along the main
sequence down to the oldest main sequence turnoff and on the subgiant branch.
Thus, in principle, it would be enough to sample the observed and mother CMD
down to the magnitude of the oldest turnoff of the more metal rich population
included in the mother CMD (as this is the faintest population), that is,
MG=4.2. However, it is reasonable to ask whether including a larger portion of
the main sequence below the oMSTO could provide useful information and
increase the accuracy or precision of the SFH derivation. To check this, we
derived the SFH of the synthetic and observed clusters described in (i) using
three bundles with different faint MG limit: 4.2, 4.7 and 5.2 and fitted 2D
gaussians to the age-metallicity distribution of each cluster to determine the
recovered age and metallicity and the corresponding $\sigma_{age}$ and
$\sigma_{Z}$. No significant difference in the precision or accuracy of the
derived ages is observed by changing the faint magnitude of the bundle. In
Appendix D, we show a compilation of tests carried out with the goal of
assessing the robustness of our SFH recovery. In particular, in Figure 26 we
show the SFH of the GCNS derived with the three faint bundle limits. It is
clearly seen that the solutions are basically identical. This is an important
finding as it implies that there is no gain, as far as the SFH is concerned,
in sampling the CMD deeper than the oldest and more metal-rich subgiant branch
involved, and thus, it will allow us to reach larger distances in the galaxy
than if a fainter magnitude would be necessary. As discussed above, this
result can be somewhat expected, as below the oMSTO, stars of different ages
are mixed.

We have also tested if using different weights according to the variations in
the range of stellar ages across the CMD would produce differences in the
resulting SFH (see Sec. 3.3, par. ii) and Figure 4). Figure 27 compares two
solutions calculated with and without weights across the CMD. It can be seen
that also in this case the results are basically identical.

####  4.3.3 Systematic differences in the Gaia DR3 magnitude scale and that of
stellar evolution models

For each model in Table 13, the corresponding best shift $(\delta c,\delta m)$
is calculated, and subtracted to each mother CMD before calculating the final
SFH. These shifts are also specified in Table 13. In principle, this shift
should be a systematic difference for each stellar evolution library, given a
set of bolometric corrections. However, different binary populations in the
mother CMD can also lead to small differences in the shifts as they change the
overall distribution of stars in colour. For this reason, we have computed the
shifts for each mother CMD and listed them in the last column of Table 13. It
can be seen that, for a given binary population and library, the shifts are
basically identical, particularly in the case of the BaSTI-IAC library, and
they change slightly (for a maximum of 0.01 in magnitude and/or colour) for
different binary populations. In general, it appears that models are bluer and
fainter than the GCNS by $\simeq 0.03-0.04$ mag.

####  4.3.4 Parametrisation of the unresolved binary population

Unresolved binaries in which the component stars are both in the main sequence
appear offset to brighter magnitudes and redder colours compared to the single
stars main sequence ridge line. Equal-mass binary stars appear offset by -0.75
mag from the locus of a single star of the same mass, while extreme mass ratio
binaries locate somewhat to the red and at a similar magnitude in the CMD
compared to the most massive star in the pair. Binary stars with intermediate
mass ratios populate a continuum of stars between these two extremes. The
unresolved binary population can be observed in the GCNS CMD below the oMSTO
as a parallel, less populated sequence, above and to the red of the main
sequence of single stars.

Fig. 6 shows the main sequence of the observational data (blue dots) and the
solution from the q01b03_60M_MG10 mother CMD (orange dots). The sequence of
unresolved binaries can be seen both in the data and in the solution. The
magnitude distribution of stars in the main sequence provides valuable
information on the characteristics of the binary star population. Comparing
the distribution in the observed CMD and that resulting in the best fit CMDs
for different $\beta$, and qmin, (see Table 13) will allow us to constrain the
parameters that result in a good fit of the binary sequence. We will then
adopt them to derive the final SFH of the GCNS.

Figure 6: Zoomed region of the main sequence below the oMSTO. Blue dots: CMD
of the GCNS. Orange dots: solution CMD obtained from the q01b03_60M_deep
model. Blue and yellow lines: fit to the main sequence of the GCNS and
solution CMD, respectively. Red line: fit to the solution CMD with no shift
applied.

In order to trace the position of the main sequence main locus (which should
approximately correspond to locus of the single stars) as a function of
colour, we fitted a Gaussian mixture model (GMM) to the distribution of
magnitude of the stars in colour bins of 0.01 mag. and, for each colour bin,
we calculated the maximum of the MG distribution as the peak MG value of the
dominant Gaussian component. We also calculated the colour average for each
bin. Thereafter, we performed a polynomial fit to these magnitudes and
colours. The blue and yellow line in Fig. 6 represent these fits for the GCNS
and of the solution CMD, respectively. The red line represents the fit to the
solution CMD with no shift applied. It is reassuring that the shift inferred
from the best fit to the bright part of the CMD (above the oMSTO) provides
also a good match between the observed and solution main sequence below the
oMSTO, which otherwise would be offset from each other (as the blue and red
lines are).

We then parametrise the distribution of stars in the main sequence compared to
its main locus through $\Delta$G=MG,pol-MG,∗ (see Gaia Collaboration, Smart et
al., 2021), such that MG,∗ refers to the absolute G magnitude of each star and
MG,pol is the absolute G magnitude interpolated in the polynomial fit for the
colour of the star. We then compute the histogram of the $\Delta$G values
(marginalising over the whole colour range).

Figure 7 depicts, in different colours, the $\Delta$G normalized histograms
derived from the solution CMDs for the SFHs derived from the ’deep’ mother
CMDs (those with MG,max=10-11) listed in Table 13, compared to that of the
GCNS (in black, dashed line). Note that all distributions are centred in zero,
by definition, and have a main component corresponding to the main sequence of
single stars. They also have a secondary bump on the left side of the main
component, centred at approximately $\Delta$G = -0.75, which corresponds to
the unresolved binary population. Its height depends on the binary fraction
$\beta$ and on qmin, with more prominent bumps for larger $\beta$ in which the
components have similar mass (larger qmin).

The observed distribution is closely matched by that corresponding to the
q01b03_60M_MG10 solution while the distribution that differs the most is that
of the q01b07_30M_MG11 model. The latter contains the largest fraction of
binaries among our tests, resulting in a stronger binary bump. The solutions
from q04b03_30M_MG11 and q01b05_60M_MG10, with intermediate fractions of
binaries lie in between171717q04b03_30M_MG11 has the same binary fraction as
q01b03_60_MG10, but the larger qmin implies a larger fraction of binaries with
similar mass, which are the ones that effectively contribute to the bump.. The
PARSEC model, q01b03_30Mparsec_MG10 has a slightly larger bump compared to the
equivalent model from BaSTI-IAC and to the observed distribution. It also
differs from the observed distribution and from that of the BaSTI-IAC models
by the shape of the faint part of the main sequence, with a more abrupt fall
to the zero value. Finally, all solutions show a few stars towards positive
values of $\Delta$G, which are associated to the presence of synthetic metal-
poor stars lying under the main sequence (see Fig. 6). The fact that this
feature is missing in the observed GCNS CMD may indicate an underlying subtle
mismatch between the empirical data and the theoretical models. In any case,
this only affects a minority of stars, so it won’t have a substantial effect
on the conclusions. Some observed stars are in this region, but they are much
more scattered in colour.

Taking into account these results, we will adopt $\beta$=0.3 with qmin=0.1 for
the final SFH of the GCNS. This value is basically compatible with the results
by Belokurov et al. (2020a), who studied the problem of unresolved binaries in
Gaia DR2 data based on the renormalised unit weight error (RUWE) parameter
from the Gaia catalogue. They exploit the fact that in the case of stars
belonging to an unresolved system, the motions of the centre of light and mass
are decoupled, and thus, assuming a single-source for the astrometric model
fails. They found that, using this ruwe parameter, they can identify such
unresolved systems. As part of their analysis, they study how the binary
fraction evolves across the CMD. From their Figure 9, we can see how the
average unresolved binary fraction, in the region of the CMD that we analyse,
ranges from $\sim$ 12$\%$ (faint main sequence) to $\sim$ 50$\%$ (bright main
sequence), with average values compatible with the 30$\%$ we are finding. A
similar result is found by Penoyre et al. (2022) for the GCNS.

Figure 7: Black dashed line: $\Delta$G histogram of the GCNS. coloured lines:
$\Delta$G histograms corresponding to solution CMDs with different $\beta$ and
qmin. All histograms have been normalised to the number of stars in the colour
range considered.

##  5 Results: the SFH derived from the GCNS.

In this section, we will discuss the deSFH and current age and metallicity
distributions of the stars within $100\,\rm{pc}$ of the Sun. We will focus on
the solution obtained with the q01b03_120M_MG5 mother CMD, which is the
largest synthetic CMD we have computed with the solar-scaled BaSTI-IAC stellar
evolution models (120 million stars with M${}_{G}\leq 5$), adopting
$\beta$=0.3, qmin=0.1 and a Kroupa IMF (Kroupa et al., 1993). We will use a
bundle with faint limit MG=4.2, weight of each CMD ’pixel’ calculated as the
inverse of the variance of the stellar ages in that ’pixel’ (see Section 3.3),
S age bins and 0.1 dex metallicity bins. A shift of ($\delta c,\delta m$) =
(-0.035, 0.04) has been subtracted to the colour and magnitude of the stars in
the mother CMD (see Table 13).

Figure 8 displays (left and middle panels) the observed and the solution CMD
of the GCNS, with the bundle including the stars that have been used for the
fit superimposed. The bundle is significantly larger than the area covered by
the observed CMD since it has to include the whole mother CMD. Since the
latter has a much larger metallicity range than the observed and solution CMD,
it covers a larger range in colour (see Figure 4). The right-hand panel shows
the residuals of the fit. Note the high quality of the fit, with no
significant trends or structures in the residuals, which in most ’pixels’ in
the CMD are within $\pm 1\sigma$ with deviations up to $\pm$ 3$\sigma$ in only
a few pixels.

Figure 8: Observed (left) and solution (middle) CMDs of the GCNS. The right
panel displays the residuals of the CMD fitting (in $\sigma$ units assuming
Poisson errors). The bundle encompassing the stars included in the fit is also
plotted in all panels.

###  5.1 The deSFH of the solar neighbourhood within 100 pc of the Sun.

Figure 9: deSFH of the stars within $100\,\rm{pc}$ of the Sun, derived with
the BaSTI-IAC stellar evolution library. Bottom left: Age-metallicity
distribution of the mass transformed into stars (astrated mass) as a function
of time and [M/H]. The iso-contours indicate the presence of a small amount of
star formation at old age and low metallicity. Upper left: star formation rate
as a function of time, SFR(t), integrated over metallicity. In the two left
panels, the inset shows and expanded view of the last 2 Gyr. Lower right:
MDFM. Upper right: cumulative distribution of the astrated mass as a function
of time. Two lines indicate the 50 and 90 percentiles. In the two lower
panels, the red line indicates solar metallicity, [M/H]=0. In all cases the
direction of the $x$ axis is such that old age/low metallicity is on the left.

Figure 9 displays what we consider our best calculation of the deSFH of the
stars currently within 100 pc of the Sun. The bottom left panel shows the age-
metallicity distribution of the mass transformed into stars as a function of
lookback time (age) and metallicity, [M/H]. Old ages are on the left. The
colour bar indicates the star formation rate in units of M☉ Gyr-1 dex-1. The
upper left panel displays the deSFR(t), which is the marginalization over
metallicity of the deSFH, in units of M☉ Gyr-1. The lower right panel displays
the metallicity distribution function of the mass transformed into stars
(MDFM), in units of M☉ dex-1. Finally, the upper right panel shows the
cumulative distribution of the mass transformed into stars as a function of
lookback time. The lines indicate the 50 and 90 percentiles.

Figure 10: Number of stars with M${}_{G}<4.2$, as derived from the BaSTI-IAC
stellar evolution library. Bottom left: Age-metallicity distribution of alive
stars within 100 pc of the Sun. Upper left: number of stars as a function of
their age. Lower right panel: MDFS. Upper right: fraction of stars as a
function of age. Two lines indicate the 50 and 90 percentiles. In the two
lower panels, the red line indicates solar metallicity [M/H]=0. The number are
calculated for stars with M${}_{G}<4.2$.

Figure 10 displays an alternative view of the characteristics of the stellar
populations in the solar neighbourhood, in terms of the number of stars
currently present as a function of their age and metallicity181818In this
figure, the number of stars is that inside the bundle, that is, stars with
M${}_{G}<4.2$. This number could be extrapolated to include the stars with
M${}_{G}>4.2$ using the IMF.. The bottom left panel shows the distribution in
age and metallicity of the stars currently alive. The upper left panel
displays the number of stars as a function of their age while the lower right
panel presents the stellar metallicity distribution function (MDFS). Finally,
the upper right panel shows the fraction of stars as a function of their age.

Note that the features in Figures 9 and 10 are very similar, the main
difference being their relative strength at young and old ages, since a
fraction of the mass that has been transformed into stars at any age is not
anymore in the form of currently alive stars, and this fraction varies with
time. For example, the peak of star formation that can be observed 3 Gyr ago
in the SFR(t) plot (Figure 9, upper left panel) has approximately twice the
intensity of the peak occurred 10 Gyr ago, while the ratio of number of stars
is approximately 4:1 (see Figure 10, upper left panel).

The most detailed and rich view of the history of the stellar mass in the
solar neighbourhood is provided by the panels showing the age-metallicity
distribution in Figures 9 and 10 (bottom left panels). We remind the reader
that, in the case of the deSFH, this is the mass (per unit time and
metallicity) that has been transformed into stars somewhere in the Galaxy to
account for the stars that are today in the studied volume. The corresponding
panel in Figure 10 describes the main features of the stellar populations
currently located in the solar neighbourhood. These age-metallicity
distributions can be described as follows:

-We measure an age of $\simeq$ 11 Gyr and a metallicity of ${\rm[M/H]}=-0.16$ for the oldest stars. Taking into account the systematic difference between input and recovered ages found in the tests with synthetic clusters (see Figure 5), this age could be reduced to about 10.4 Gyr. Thereafter, a sequence of progressively younger and more metal rich stars is observed, up to age $\simeq$ 9.7 Gyr (or 9.2 Gyr considering the possible systematics) and ${\rm[M/H]}=0.25$.

-Between $\simeq 9.5$ and 6 Gyr ago, two main metallicity sequences, somewhat disjoint in time, exist: at an earlier time, stars have metallicity slightly above solar, while after $\simeq$ 8 Gyr ago, the main population has solar metallicity, and a less prominent population with ${\rm[M/H]}=-0.25$) can be observed.

-Six Gyr ago a dip in the SFR(t) can be observed, followed by a large scatter in metallicity, with three main populations: one with [M/H]$=-0.4$, another at solar metallicity, and a third one with a very narrow age range (slightly older than 4 Gyr) and super-solar metallicity. The latter one coincides with a star formation gap at the expected solar metallicity, and is followed by a conspicuous break of star formation (see the upper left panel displaying the deSFR(t)).

-Since 4 Gyr ago, the star formation proceeds at an average higher rate until the present time, even though still showing a bursty behaviour. Between 4 and 2 Gyr ago, populations at different metallicities, solar and below solar (${\rm[M/H]}\simeq-0.2$) coexist. Finally, since 2 Gyr ago, the majority of the stars have solar metallicity, with a small fraction of super-solar stars.

-Almost no star formation is detected with [M/H]$\lesssim-0.5$, even though the whole metallicity range between [M/H]$=-2.2$ and [M/H]=0.45 was present in the mother CMD used to determine the SFH, and no a priori chemical evolution law was imposed. Only a very small amount of star formation is detected with ages between 12 and 10 Gyr and metallicity around [M/H]$\simeq-1.5$ (it is barely visible in Figures 9, see contours, and 10). It translates in around one hundred stars in this range of age and metallicity which can be only hinted in the left lower panel of Figure 10. We will discuss this population in Section 6. This paucity of low metallicity stars is not surprising since the GCNS sample is primarily composed of thin disk stars: more than 90% of the stars that we can classify based on kinematics are found to belong to the thin disk; see Figure 1 and Appendix A.

The right lower panels of Figures 9 and 10 show that the metallicity
distribution function peaks very close to solar metallicity with the great
majority of the stars having [M/H] in the range 0.2 to -0.3, a minority
population down to metallicity [M/H]$\simeq$ -0.5, and an even less relevant
population up to [M/H]$\simeq$0.4. The very small old metal poor population
that we mentioned earlier can be noticed as a small peak around
${\rm[M/H]}\simeq-1.3$

Finally the upper right panels of Figures 9 and 10 show the fraction of the
total mass transformed into stars as a function of age, and the fraction of
stars older than a given age, respectively. The lines indicate the 50 and 90
percentiles: it can be seen that 50% and 90% of the stellar mass was reached
4.5 and 0.9 Gyr ago, respectively, while 50% and 90% of the stars are older
than 3 and 0.6 Gyr. These numbers picture an average young solar
neighbourhood.

###  5.2 Robustness of the solution

In Section 4 we have discussed a number of choices that have been tested,
regarding both parameters adopted to calculate the mother CMDs and the input
parameters of the fit, in order to determine their impact on the final SFH.
These parameters are: the number of stars in the mother CMD, the size of the
age bins, the faint magnitude limit of the bundle, the weights across the CMD,
and the unresolved binary stars parameterisation. Additionally, a mother CMD
has been computed with the PARSEC stellar evolution models. SFHs determined
with different choices in each one of these aspects are shown in a number of
figures in Appendix D. Here, for its relevance, in Figure 11 we show the
solution obtained with the PARSEC stellar evolution models.

Figure 11: deSFH of the stars within 100 pc of the Sun derived with the PARSEC
synthetic CMD. Bottom left: Age-metallicity distribution of the astrated mass
as a function of time and [M/H]. The inset shows an expanded view of the last
2 Gyr. The iso-contours indicate the presence of a small amount of star
formation at old age and low metallicity. Upper left: star formation rate as a
function of time, SFR(t), integrated over metallicity. Lower right: MDFM.
Upper right: cumulative distribution of the astrated mass as a function of
time. Two lines indicate the 50 and 90 percentiles. In the two lower panels,
the red line indicates solar metallicity [M/H]=0. In all cases the direction
of the $x$ axis is such that old age/low metallicity is on the left.

Concentrating first in the various solutions obtained with the BaSTI-IAC
stellar evolution library, we see that the main conclusions regarding the SFH
of the GCNS do not change with the different parametrisations which, from
oldest to youngest ages, can be summarised as follows (see Figure 9: i) the
oldest substantial191919Very minor older populations are also found, with low
and high metallicity, [M/H]$\simeq-1.5$, and [M/H]$\simeq$0.3. stellar
populations are measured as $\simeq$11.5 Gyr old; ii) there are several
episodes of star formation with decreasing age and increasing metallicity in
the subsequent $\simeq$ 2 Gyr, until 9-9.5 Gyr ago; iii) there are hints of
splitting of the age-metallicity distribution between 10 and 6 Gyr ago, with a
large scatter and different episodes of star formation with different
metallicities around $\simeq$5 Gyr ago; iv) there is a sudden decrease in the
star formation activity around 4 Gyr followed by an increased average star
formation rate, which lasts until the present time; v) the star formation
activity is bursty at all ages.

Most of these features hold in the solution with the PARSEC models, which we
show in Figure 11, where the BaSTI-IAC solution has been superimposed in thin
dashed lines in the panels displaying the deSFR(t), the MDFM, and the
cumulative distribution of the mass transformed into stars. The main
differences are that the overall population is retrieved more metal poor with
PARSEC (see the MDFM in the lower right panel of Figure 11), and that the
metallicity range for each given age is narrower, with less obvious
metallicity splits. However, the range of metallicity around 10 Gyr ago is
also present, as well as an split around 8 Gyr ago. A feature that is hinted
in the BaSTI-IAC solutions and is more evident in the PARSEC one is a somewhat
decreasing metallicity trend since the oldest age ($\simeq 11$ Gyr ago) until
4 Gyr ago. At this age, after the same break in the SFH, the metallicity
starts to increase to reach solar metallicity at the present time. Comparing
the deSFR(t) in the upper left panel of Figure 11, it can be noticed that the
same episodes of star formation are present in both cases, even though the
relative importance of some of them may be somewhat different, and age shifts
are observed, specially between 9 and 4 Gyr ago. It can also be noticed that
no significant star formation before 11 Gyr ago is measured with the PARSEC
models, while star formation starts to rise 12 Gyr ago in the BaSTI-IAC
solution, with metallicity slightly below solar. Finally, the shape of the
cumulative SFR after 10 Gyr ago indicates a slightly more constant average SFR
in the case of the solution based on the PARSEC library, with 50% and 90% of
the stellar mass reached 5 and 1 Gyr ago - as opposed to 4.5 and 0.9 Gyr for
the case of the BaSTI-IAC model library.

Going back to the comparison of the different SFHs obtained with the BaSTI-IAC
models, some differences in the solutions can be associated with changes in
some of the parameters: a larger number of stars in the mother CMD (Figure 25)
or smaller age bins (Figure 28) result in sharper features in the SFH. The
effect of the age bin size is most evident. For example, while there are hints
of the split in the age-metallicity distribution at intermediate-ages
regardless of the size of the bins, they are more clearly visible with the
smallest (S) bins. The reality of this latter feature, which appears in most
of the solutions with different parameters, deserves further attention. In
Section C.3 we have discussed an experiment designed to check whether this
metallicity split may be real: we derived the SFH for a composite CMD of 14
synthetic clusters of 7 different ages and two metallicities separated by a
small metallicity gap (0.15 dex). Figure 24 displays the age-metallicity
distribution recovered from the CMD of this composite population. It can be
seen that the two metallicities are best separated for ages of 8 Gyr or older,
where the metallicity is more precisely recovered (see also Figure 22), and of
course, they are best separated also for the smallest (S) age bins. For
younger ages, the existence of the two metallicities is still hinted in the
cluster’s age-metallicity distribution, but in a less obvious way. We
conclude, thus, that the different metallicity sequences seen in the SFH of
the GCNS are real, at least for stars older than 2 Gyr, age when the
metallicity distribution gets narrower.

###  5.3 Comparison with spectroscopic metallicity distribution functions.

The comparison of the MDFS resulting from the SFH calculation with that
obtained from spectroscopy can provide a totally independent, external check
of our results, as no assumptions on the age-metallicity distribution or the
metallicity distribution function have been made in the SFH derivation.

Figure 12 shows the comparison of our derived MDFS with two spectroscopic MDF:
that for the metallicities derived by $Gaia$ GSP-Spec for the same volume
(middle panel; see Section 2), and that for the Fuhrmann et al. (2017) volume
complete sample of stars located within 25 pc of the Sun (right panel) and
with T${}_{eff}\geq 5300$ K. We show both the MDFS from our preferred deSFH
obtained from the BaSTI-IAC mother CMD displayed in Figure 9 (solid coloured
lines) and that obtained from the PARSEC mother CMD, displayed in Figure 11
(dashed coloured lines).

Figure 12: Comparison of the MDFS derived from CMD-fitting (both with the
BaSTI-IAC and PARSEC stellar evolution models, with two spectroscopic MDF.
Left panel: completeness in the CMD of the metallicity measurements from
$Gaia$ GSP-Spec (blue-yellow histogram) and CMD of the stars measured by
Fuhrmann et al. (2017), in brown, superimposed on the CMD of the GCNS, in
grey. Middle panel: MDF obtained from $Gaia$ GSP-Spec metallicities (Recio-
Blanco et al., 2023) for stars in the GCNS. Right panel: MDF derived from
Fuhrmann et al. (2017) stars within 25 pc of the Sun.

For the comparison with the Gaia GSP-Spec sample, we took into account the
completeness and typical errors of the spectroscopic metallicity measurements
within the GCNS volume. The completeness is represented in the left panel of
Figure 12, indicated as the number of stars with spectroscopic [Fe/H] divided
by the total number of stars in the GCNS, in pixels of ($\delta$ col, $\delta$
mag)= (0.05, 0.1) mag across the CMD. The GCNS CMD is shown in the background
as grey dots. It can be seen that there is a sharp cut at $G_{BP}-G_{RP}\simeq
0.75$ corresponding to the T${}_{\rm{eff}}=8000$ limit of GSP-Spec (Creevey et
al., 2023; Recio-Blanco et al., 2023). This cut implies that main sequence
stars younger than $\simeq$2 Gyr are not represented in the Gaia GSP-Spec MDF
(MDFGaia thereafter). The figure also shows that RGB stars have Gaia
metallicity determinations with high completeness, while the completeness is
around 50% for low main sequence stars. Given that, according to our derived
age-metallicity distributions, the metallicity range of stars younger than
$\simeq$ 2 Gyr is within the limits of the general metallicity range, and that
stars up to $\simeq$ 1 Gyr old are present in the RGB, we do not expect that
the incompleteness of the spectroscopic metallicity measurements across the
CMD will significantly affect the shape of the MDFGaia. In any case, for
consistency, we simulated the Gaia GSP-Spec incompleteness and errors as
explained below.

In the middle panel of Figure 12 the MDFGaia (black line) is compared with our
MDFS. The curves are histograms of the number of stars per 0.05 dex bins,
normalised to the total number of stars. The green lines corresponds to the
metallicity distribution of the whole solution CMD, while the pink lines are
the same MDFS after completeness correction, MDF${}_{S}^{C}$. To perform this
correction, we randomly selected, for each CMD pixel, a fraction of stars from
the solution according to the completeness. As expected, the two lines are
basically identical, confirming that the metallicity distribution across the
CMD is basically uniform. Note that the peak of the MDFS and the MDFGaia occur
basically at the same metallicity in the case of the BaSTI-IAC solution (solid
line), indicating that the metallicity scale of our solution from these models
is the same as the GSP-Spec metallicity scale202020Slightly different $(\delta
c,\delta m)$ shifts would result in a slightly different peak of the age-
metallicity distribution, but it is reassuring that our best estimate of this
shift provides a solution that leads to a remarkable agreement with the
spectroscopic MDFGaia. This is not the case for the PARSEC models (dashed
line): in this case, the peak metallicity of the solution MDFS is more metal
poor than that of MDFGaia. The difference in shape of the distributions can be
attributed to the different dispersion of the metallicity measurements. In the
case of the SFH, the average standard deviation of the metallicity derived
from the open clusters (for the S age bins) is $\sigma_{[M/H]}^{OC}\simeq$0.06
dex. The standard deviation of the differences between the GSP-Spec
metallicity measurements, based on Table D.1. in Recio-Blanco et al. (2023),
is around 0.1 - 0.12 dex. In order to simulate this different dispersion in
the MDFS, we used a more conservative value of 0.15 dex, that we call
spectroscopic sigma. We then add a normal distribution of noise to the
metallicity values of the solution, with a standard deviation equal to the
square root of the quadratic difference of the spectroscopic sigma minus
$\sigma_{[M/H]}^{OC}$. The blue lines are the MDF${}_{S}^{C}$ after adding
this extra dispersion (MDF${}_{S,d}^{C}$). Note that the width of
MDF${}_{S,d}^{C}$ are now practically identical to that of MDFGaia.

The right panel of Figure 12 displays the comparison of our MDFS (from BaSTI-
IAC, solid green line, and PARSEC, dashed green line) with the MDF obtained by
Fuhrmann et al. (2017). It is based in high-resolution spectroscopic data
obtained with the FOCES (Pfeiffer et al., 1998) and BESO (Steiner et al.,
2006) échelle spectrographs for the northern and southern samples,
respectively. The fainter sources were observed with the FEROS échelle
spectrograph (Kaufer et al., 1999). This catalogue provides very high quality
[Fe/H] and [Mg/Fe] determinations, with typical errors of 0.08 dex in [Fe/H].
Since these errors are similar to our estimated errors in the metallicity
determination of the clusters, we do not apply any extra scatter to our MDFS.
However, since our metallicity values represent the global metallicity of the
stars, [M/H], we have transformed the [Fe/H] values into [M/H] using [Mg/Fe]
as a proxy for [$\alpha$/Fe], and following the relation

$\mathrm{[M/H]}=\mathrm{[Fe/H]}+\log\left(0.694\times
10^{\mathrm{[\alpha/Fe]}}+0.306\right)$ (1)

from Salaris et al. (1993), with coefficients slightly different from the
original reference in order to take into account the reference solar mixture
adopted in the BaSTI-IAC models, based on the heavy element distribution
provided by Caffau et al. (2011). The agreement between Fuhrmann et al. (2017)
MDF and our MDFS derived from the BaSTI-IAC solution is outstanding: the peak
metallicity and the width and shape of the distributions are basically
identical, including its asymmetry slightly biased to lower metallicities.
Additionally, both MDF show a small peak at [M/H] $\simeq-0.4$, which was
erased when degrading our metallicity errors to match those in the MDFGaia. In
Fuhrmann et al. (2017), most stars with this metallicity belong to the low
[Mg/Fe] branch, with ${\rm[Mg/Fe]}\simeq-0.15$, and are consequently
classified as Pop I stars. In our solutions, the [M/H]$\simeq-0.4$ stars have
ages around 5 Gyr 212121See also the position of the 5 Gyr old, [M/H]=-0.5
isochrone compared to the thin and thick disks CMDs in Figure 1. This striking
agreement of our MDFS with a high quality spectroscopic MDF supports that the
interesting low-metallicity feature in the age-metallicity distribution at
intermediate ages (see Figures 9 and 10) and the many structures and splits in
our age-metallicity distributions are real, including possibly the bump at the
high metallicity end, [M/H]$\geq 0.2$, which is present also in the Fuhrmann
et al. (2017) MDF. Even though the shape is different from the high
metallicity feature in our MDFS, the fraction of stars encompassed in this
feature is similar, and thus, this also supports that the high metallicity
features that can be observed in our solution may be real. The most striking
of those features is the one at age $\simeq$ 4 Gyr and [M/H]=0.4.

The agreement between MDFS obtained from the PARSEC solution and the Fuhrmann
et al. (2017) MDF is clearly worse: the peak metallicity does not coincide
and, more importantly, the width and the fine structure of the MDF are in
worse agreement. This reinforces our choice of the BaSTI-IAC library as the
reference stellar evolution library for this work.

###  5.4 Comparison with literature age-metalliicty distributions

Figure 13: Age-metallicity distributions from the literature. The red line
indicates solar metallicity. Panels from top to bottom show results from
different sources: $Dir$SFH results within 100 pc obtained in section 5, with
error bars indicating $\sigma_{age}$ (see Figure 5); StarHorse results for
LAMOST MRS DR7, APOGEE DR17 and GALAH DR3; Xiang & Rix (2022) results for
LAMOST LRS DR7; and Kordopatis et al. (2023) results for $Gaia$ GSP-Spec. Each
panel shows the number of stars used to derive the age-metallicity
distribution, constrained to similar ranges of Rg, Rapo and Rperi as the
sample presented in this paper. Error bars indicate the mean uncertainty in
the age (as provided for each catalogue) within specific ranges (0-3, 3-6,
6-9, 9-12, 12-14 Gyr).

Figure 13 displays a number of age-metallicity distributions adapted from
recent literature data (Xiang & Rix, 2022; Queiroz et al., 2023; Kordopatis et
al., 2023), as indicated in the labels and on the figure caption. The samples
analysed in these works are very heterogeneous not only in terms of the volume
covered, but also in the data used (from different spectrocopic survey
samples) and the resulting selection functions. They all have in common the
use of Bayesian isochrone fitting techniques to derive ages. The precision in
the ages determined for individual stars depends on the accuracy of the
spectroscopic measurements and the stellar evolutionary state, as highlighted
in previous studies (Mints & Hekker, 2017; Sanders & Das, 2018). Isochrone
fitting methods are inherently linked to the choice of stellar evolution
models. In Figure 13, ages are derived by comparison with the PARSEC models
(used in Queiroz et al. 2023 and Kordopatis et al. 2023) and the Yonsei-Yale
(Demarque et al., 2004) models (used in Xiang & Rix, 2022), which introduce
differences in age scales compared to the BaSTI-IAC models used in the
$Dir$SFH determination. In addition, each Bayesian method incorporates
different prior assumptions. Furthermore, APOGEE (Majewski et al., 2017),
GALAH (De Silva et al., 2015), LAMOST (Cui et al., 2012), and Gaia RVS (Recio-
Blanco et al., 2023) each incorporate specific assumptions in the derivation
of spectroscopic [Fe/H] and alpha abundances, which we present in terms of
global metallicity using eq.1.

In order to mitigate the effect of the different selection functions and make
the Galactic stellar populations represented in the different panels more
comparable to those present in the volume of 100 pc radius analysed in this
paper (top panel), we have selected, from Kordopatis et al. (2023) only those
stars with guiding Radius $7\leq R_{g}\leq 9$ and Z${{}_{\rm max}}\leq 500$
pc, similar to the majority of the stars in our sample (see Figure
14)222222The orbital parameters represented in Figure 14 are taken from
Kordopatis et al. (2023). Then, we have selected the same stars in the Xiang &
Rix (2022) and Queiroz et al. (2023) samples.

The features in the four central panels are very similar and characterised by
a basically flat age metallicity relation with metallicity slightly below
solar and prominent maximum in the stellar distribution around 3 Gyr ago. The
metallicity distribution seems to get slightly narrower around 4 Gyr ago with
hints of widening between 6 and 4 Gyr ago. The absence of stars with ages
younger than $\simeq$ 2 Gyr is entirely due to selection effects, mainly
caused by the Teff and logg limits of the stars actually present in these
catalogues, and the further limits imposed in the papers for the stars
actually used to derive the ages. The bottom panel, corresponding to
Kordopatis et al. (2023) ages determined from Gaia RVS data, displays the most
numerous stellar sample, and it also shows a high concentration of stars
around 3 Gyr ago. However, in this case, a more prominent maximum is seen
around 1.5 Gyr ago. This difference with the other panels can be attributed to
the fact that (Xiang & Rix, 2022) and Queiroz et al. (2023) restrict their
analysis to subgiant branch stars or main sequence turnoff and subgiant branch
stars, respectively, while Kordopatis et al. (2023) don’t have this
restriction and thus, include in their sample brighter, and thus younger
stars, as well as RGB stars. Note however the large errors of this sample,
specially at young ages.

These age-metalliticy distributions provide much less detail, owing to larger
age errors, than the one derived in this paper from CMD-fitting, and displayed
in the top panel. A detailed comparison of the age-metallicity distributions
obtained with both methodologies is out of the scope of this work. However, it
is possible to identify the overdensity of 3 Gyr old stars with below solar
metallicity with the various star formation events that we resolve between 4
and 2 Gyr ago. The younger overdensity in the Kordopatis et al. (2023) age-
metallicity distribution coincides with the star formation episodes that we
identify in the last 1 Gyr. Finally, the slightly increased metallicity
dispersion between 6 and 4 Gyr ago could be the same that we observe, resolved
in three star formation episodes at different metallicities.

This comparison shows the superior ability of CMD fitting to derive age
distributions compared to those obtained from the combination of individual
ages. This is thanks to the fact that CMD-fitting uses, in addition to the
colours and magnitudes of the stars, the extra information on the relative
lifetimes of stars of different masses and in different evolutionary stages,
provided by stellar evolution theory, and on the relative number of stars of
different masses, according to externally derived IMF. It also has the
advantage that the required data are only photometry and parallaxes from Gaia,
and thus the selection function is easier to control (and in this particular
case, the sample can be considered complete). Of course, the important caveat
is that ages are not obtained for each individual star, and this precludes the
direct association of ages to other important properties, such as chemical
abundances or kinematics, which help disentangling different stellar
populations in the MW. Both sets of techniques are therefore highly
complementary and necessary to advance in our understanding of the
evolutionary processes that have shaped the MW.

##  6 Discussion: the evolutionary history of the Milky Way thin disk

Figure 14: Orbital parameters for stars in the GCNS. Upper left panel: guiding
radius (Rg); upper right panel: Zmax; lower left panel: apocentric radius
(Rapo); lower right panel: pericentric radius (Rperi) radius.

Figure 14 displays the guiding radius (Rg, upper left panel), Zmax (upper
right panel), apocentric (Rapo) and pericentric (Rperi) radius (lower left and
right panels, respectively), for stars in the GCNS. We have adopted for these
parameters the values provided by Kordopatis et al. (2023). It can be seen
that the majority of the stars currently within 100 pc of the Sun have guiding
radius between 7 and 9 kpc and Zmax up to 0.5 kpc. In their apocenters they
can travel as far as 11 kpc while in their pericenters they can approach the
Galactic centre as close as 5 kpc. We live, therefore, in a quite cosmopolitan
place, with a good fraction of the inhabitants possibly born quite far away.
As a result, the age-metallicity distribution derived in this paper from the
exquisitely observed bubble of 100 pc radius from the Sun contains, in fact,
information of a much larger volume. This local stellar population, however,
is heavily dominated by thin disk stars, as discussed in Section 2 and
graphically shown in Figure 1. Indeed, our kinematic selection classifies only
$\simeq$ 10% of the GCNS stars as thick disk. Most of the stars have low
[$\alpha$/Fe], and so they would be chemically classified as thin disk as
well. The results in this paper, therefore, provide a sharp view of the
formation and evolution of the MW thin disk, while still providing hints of
the nature of the stellar populations in the thick disk and the halo.

We shall start our discussion in chronological order, by noting the small
amount of stars older than 10 Gyr, with [M/H]$\lesssim-1$, which are
metallicity values typical of the stellar halo. This population is shown with
contours in Figure 9. Based on previous derivations of the SFH of
kinematically selected halo stellar populations (Gallart et al., 2019a;
Massari et al., 2023), the stars with this age-metallicity location would
correspond to the GES remnant. Interestingly, we don’t see in our solution
signs of a population older than $\simeq 11$ Gyr and [M/H] between -0.5 and
-0.3 which according to Gallart et al. (2019a) and Massari et al. (2023) would
correspond to the red sequence of the halo CMD (see also Haywood et al.,
2018b) also called the Splash (Belokurov et al., 2020b). In our solution CMD,
for M${}_{G}\leq 4.2$, 96 out of 21357 stars232323When quoting number of stars
in a particular population, we refer to stars inside the analysed area of the
CMD. have [M/H]$<-0.6$ and ages older than 9 Gyr. However, when classifying
the stars in the GCNS based on their kinematics, only 9 stars end up being
classified as halo within this magnitude range. Therefore, most of the stars
that produce the signal that we observe at old age and low metallicity must
have disc kinematics. This would be in qualitative agreement with the numerous
claims in the literature about metal-poor stars in disc-like orbits, possibly
tracing an early rotating disk component (Morrison, 1990; Haywood et al.,
2013; Chiba & Beers, 2000; Carollo et al., 2010; Sestito et al., 2020; Di
Matteo et al., 2020; Fernández-Alvar et al., 2021; Cordoni et al., 2021;
Mardini et al., 2022; Xiang & Rix, 2022; Belokurov & Kravtsov, 2022; Feltzing
& Feuillet, 2023; Chandra et al., 2023; Nepal et al., 2024b; Fernández-Alvar
et al., 2024) and opens interesting prospects about the possibility of
studying the ages of stars in the metal poor disk by analysing larger samples
of stars pertaining to different Milky Way volumes. However, owing to the
small number of stars in this particular sample, our results on this matter
should be considered as preliminary.

The thin disk is expected to form after the early epoch of frequent mergers
and to grow from the smooth accretion of gas (Brook et al., 2012; Agertz et
al., 2021; Yu et al., 2023; Grand et al., 2017). Since most stars in our
sample belong to the thin disk, our results should be able to provide tight
constraints regarding the oldest age of the bulk of the Milky Way thin disk
stars. Our deSFH indicates that these stars have near solar metallicity and
are $\simeq$ 11 Gyr old (after taking into account our age precision at old
ages, and $\simeq$ 10.4 Gyr old after additionally correcting by age
systematics).242424This is actually an upper limit, as the 11 Gyr old
population in Figures 9 and 10 could actually originate in the stars that we
have kinematically classified as thick disk in Section 2. However, deriving
SFHs for kinematically distinct populations is beyond the scope of this paper,
and will be the subject of Fernández-Alvar et al. (2024, in prep). This age
indeed coincides very well with the time of the last major merger experienced
by the Milky Way, that of Gaia-Enceladus-Sausage (GES, Belokurov et al., 2018;
Helmi et al., 2018) dated to have occurred around 11-10 Gyr ago (Gallart et
al. 2019a, Bonaca et al. 2020, Montalbán et al. 2021, although see also Donlon
et al. 2023). The high metallicity of the oldest stars in the thin disk must
be the result of a rapid initial chemical enrichment in our Galaxy in an epoch
of intense star formation activity expected to take place in its inner parts
(in an inside-out scenario in which star formation is more efficient in the
innermost regions, e.g. Larson, 1976; Chiappini et al., 1997, 2001; Haywood et
al., 2018a; Matteucci, 2021), and it is followed by a trend of slightly
decreasing metallicity with time. Inverted age-metallicity relations have been
inferred by other works (Sahlholdt et al., 2022; Xiang & Rix, 2022; Imig et
al., 2023, see also Section 5.4), and different explanations have been
proposed. The most natural one could be that this trend is the result of more
steady star formation fuelled by gas accretion while the radius of the MW disk
progressively grows (Chiappini et al., 2001; Minchev et al., 2014; Ciucă et
al., 2021). In this scenario, the combined effects of gas dilution and metal
enrichment by ongoing star formation may be the origin of the observed trend
(Weinberg et al., 2017) for the main populations formed between $\simeq$ 11
and 6 Gyr ago. The effects of radial migration (Grand et al., 2015; Minchev et
al., 2018; Okalidis et al., 2022) and satellite infall (Grand et al., 2018;
Agertz et al., 2021; Lu et al., 2022, see below) are expected to also play a
role.

A dramatic event in the history of the Milky Way must have happened around 6
Gyr ago, when a dip in the star formation activity is observed, followed by a
clear disruption of the previous relatively steady age-metallicity trend. This
event may be the first infall of the Sagittarius dwarf galaxy (Sgr), thought
to have occurred around this epoch (Law & Majewski, 2010; Laporte et al.,
2018), and seen to have possibly caused repeated bursts of star formation
coinciding with its pericentric passages (Ruiz-Lara et al., 2020a).
Cosmological simulations find that accretion of large satellites is usually
accompanied by the accretion of gas (both from the gas in the satellite itself
252525Sagittarius must have been gas-rich at the time of its accretion, as
star formation close to the centre of the remnant galaxy has been seen to
proceed for many Gyr after its infall, possibly coinciding with the subsequent
pericentric passages (see Siegel et al., 2007). and from cosmological cold
flows, Ruiz-Lara et al. 2016, Grand et al. 2020, Agertz et al. 2021), and
result in a drastic chemical and kinematic reorganisation of the disk,
accompanied by disk growth. For example, Grand et al. (2018) found that the
accretion of a gas-rich satellite about 6 Gyr ago supplied fresh, metal-poor
gas which reinvigorated star formation in the outer disc. This resulted in a
resumption of disc growth and distinct features in the age-metallicity
distribution. Another example is the work of Agertz et al. (2021), who find
that a satellite accretion occurring at $z\simeq 1.5$ ($\simeq$ 9 Gyr ago)
results in the rapid formation of an extended disk which is more metal poor
than the existing disk when the satellite was accreted, leading to a feature
in the age-metallicity distribution of their simulated disk that is highly
reminiscent of the low and intermediate metallicity population that we infer
around 5 Gyr ago. What is missing in their predicted age-metallicity
distribution is the high-metallicity population that we infer in ours. In our
case, this may indicate that the star formation triggered by this accretion
event may have reached the inner galaxy as well, where star formation would
take place in a more metal enriched environment. Stellar migration to the
solar radius could possibly be enhanced at this time by a growing bar (Nepal
et al., 2024a) which formation could have been triggered by the interaction
with Sgr itself (e.g. Purcell et al. 2011; however see also Merrow et al. 2023
for claims that the GES could be linked to the bar formation). Thus, both
populations with metallicity clearly different (higher and lower) than solar,
observed between 6 and 4 Gyr ago, could be present in the analysed volume
owing to radial migration (from the inner and outer disk, respectively). It is
interesting to note that the three populations follow an age sequence, with
the most metal poor population being the oldest, and the more metal rich, the
youngest. Additionally, the number of stars in them decreases from the more
metal-poor to the most metal-rich population (with $\simeq 1000,700$ and 400
stars respectively). This could indicate a propagation of the star formation,
possibly with decreasing intensity, from the outer to the inner disk over a
period of 1.5-2 Gyr.

Figure 15: CMD of the two open clusters contained within 100 pc: Hyades and
Coma Berenices. The superimposed isochrones correspond to the age of the two
main star formation episodes in the last 1 Gyr, occurred $\simeq$ 900 and 650
Myr ago

After the striking dip observed 4 Gyr ago in the deSFH (Figure 9) and the age-
metallicity distribution (Figure 10), star formation in the extended solar
neighbourhood starts again and proceeds at an average higher rate than ever
before, and in an episodic manner. The stars formed in the first episode after
the 6-4 Gyr old break have metallicity very close to (or just slightly below)
solar, but a further metallicity decrease is observed between 3 and 1.5 Gyr
ago, with most of the stars having [M/H] $\simeq-0.2$ and only a minor
population with solar metallicity. This age can be linked to another
pericentric passage of Sgr (Purcell et al., 2011; Laporte et al., 2018, with
exact timing depending on the models), and even to a moment in which the MW
starts to experience the combined influence of both Sgr and the LMC (Besla et
al., 2007; Laporte et al., 2018; Vasiliev, 2023). Therefore, the same
mechanisms discussed above, with possibly less intensity owing to a stripped,
less massive Sgr may be at play here. Finally in the most recent 2 Gyr, the
deSFH exhibits a notable bursty pattern (which is observable due to enhanced
age resolution in this time range) and a discernible upward trend in
metallicity, with the majority of stars having metallicity very close to
solar. This pattern can be attributed to both the gradual increase in local
metallicity and the limited radial migration that likely occurred within this
relatively short time span, with the majority of stars in the analysed volume
born in close proximity to their current locations (Minchev et al., 2013). An
exception to this is the significant $\simeq$ 600 Myr old, metal-rich
population, containing around 250 stars (see inset in the lower left panel of
Figure 9), which presence in the solar neighbourhood must again be attributed
to radial migration. Finally, the ages of the two strongest episodes of star
formation that we derive in the last 1 Gyr, occurring $\simeq$ 650 and 900 Myr
ago, coincide very approximately with the ages of the two open clusters
contained within 100 pc, Hyades and Coma Berenices, respectively (see Figure
15 and Gaia Collaboration, Smart et al., 2021). However, the presence of the
stars belonging to these clusters in our sample cannot account for the whole
strength of the bursts, which contain $\simeq$ 800 stars each, while the
number of stars in the cluster CMD within the same magnitude range
(M${}_{G}\leq 5$) is an order of magnitude smaller.

##  7 Summary and conclusions

This paper presents for the first time a detailed dynamically evolved SFH
(deSFH) within 100 pc of the Sun, as well as the distribution of age and
metallicity of the stars in the same Milky Way volume. They have been derived
with a new set of codes to perform CMD fitting, called CMDft.Gaia. In this
section we will present our main conclusions both regarding the stellar
content of the solar neighbourhood and the performance of CMDft.Gaia.

###  7.1 The stellar populations within 100pc of the Sun

The main conclusions regarding the evolutionary history of the disk near the
solar radius, obtained from the deSFH and age and metallicity distributions of
the GCNS, which contains mostly thin disk stars inhabiting a range of
galactocentric radius between at least 7 and 9 kpc, and Zmax¡500pc, can be
summarised as follows:

-The bulk of the stellar population contained in the analysed volume form an age-metallicity sequence, with metallicities around solar. The oldest stars in this main population have an age of about 11 Gyr (or 10.4 Gyr considering systematic effects) and a metallicity of ${\rm[M/H]}=-0.16$, increasing to ${\rm[M/H]}=0.25$ within the next Gyr.

-Our solution requires around 100 stars older than 10 Gyr and with metallicities around ${\rm[M/H]}=-1.5$. This number is one order of magnitude larger than the stars that we classify as halo based on their kinematics, and thus we could be dating the population of the metal poor disk stars found in many previous works.

-Between $\simeq 9.5$ and 6 Gyr ago, two distinct metallicity sequences emerged. In the earlier phase, stars exhibited metallicity slightly above solar. Subsequently, around 8 Gyr years ago, the primary population has solar metallicity, with a less prominent group with ${\rm[M/H]}\simeq-0.25$).

-Stars between 6 and 4 Gyr old are clumped in three well separated populations with little difference in age but large metallicity differences: the oldest with [M/H]$=-0.4$, the next with solar metallicity, and the youngest with a very narrow age range $\simeq 4$ Gyr old and supersolar metallicity. The latter one coincides with a star formation gap at the expected solar metallicity, and is followed by a conspicuous break of star formation.

-Stars younger than 4 Gyr are more numerous (per age range) compared to older ages, and form clumps as a function of age and metallicity: between 4 and 2 Gyr ago, populations with a range of metallicities, including solar and sub-solar (${\rm[M/H]}\simeq-0.2$), coexisted. However, since 2 Gyr ago, the majority of stars have exhibit solar metallicity.

-The stellar metallicity distribution derived with CMDft.Gaia is in excellent agreement with spectroscopically derived ones. In particular, it compares extremely well (even in its details) with that obtained by Fuhrmann & Chini (2021) with their high-resolution, high signal-to-noise, complete sample of stars within 25 pc of the Sun. The peak metallicity and width of the two distributions are virtually identical, and the same two lower metallicity components, with ${\rm[M/H]}\simeq-0.2$ and $\simeq-0.5$, are observed. Interestingly, our age-metallicity distributions disclose that these are intermediate-age populations, including the ${\rm[M/H]}\simeq-0.5$ component, with metallicity typical of the thick disk, which here we infer to be $\simeq$ 5 Gyr old.

The picture that emerges from the results in this paper is thus that of a thin
disk that begins to form the bulk of its stars $\simeq$ 11-10 Gyr ago, from
gas chemically enriched in an earlier epoch of high star formation rate,
possibly triggered by mergers (Gallart et al. 2019a, see also Orkney et al.
2022, Belokurov & Kravtsov 2022). Several Gyr of star formation with slightly
decreasing metallicity follow, possibly the result of steady star formation
fuelled by gas accretion while the radius of the Milky Way disk progressively
grows. Then, around 6 Gyr ago, some dramatic event seems to be required to
cause the break in the star formation activity after which three stellar
populations with very different metallicities are observed. This event is
likely the first infall of the Sgr dwarf galaxy, and we hypothesise, in
agreement with results of cosmological simulations, that it was accompanied by
a substantial amount of gas accretion which led to an important chemical and
kinematic reorganisation of the disk, accompanied by further disk growth.
Stellar migration to the solar radius could possibly be enhanced at this time