# Uncovering the birth of the Milky Way through accurate stellar ages with Gaia

Carme Gallart, Edouard Bernard, Chris Brook, Tomás Ruiz-Lara, Santi

Cassisi, Vanessa Hill, Matteo Monelli

► To cite this version:

Carme Gallart, Edouard Bernard, Chris Brook, Tomás Ruiz-Lara, Santi Cassisi, et al.. Uncovering the birth of the Milky Way through accurate stellar ages with Gaia. Nature Astronomy, 2019, 3 (10), pp.932-939. 10.1038/s41550-019-0829-5. hal-02380439

HAL Id: hal-02380439

https://hal.science/hal-02380439v1

Submitted on 8 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Uncovering the birth of the Milky Way through accurate stellar ages with Gaia

Carme Gallart<fnref n="1" />,²*, Edouard J. Bernard³, Chris B. Brook<fnref n="1" />,², Tomás Ruiz-Lara<fnref n="1" />,², Santi Cassisi⁴,⁵,
Vanessa Hill³ and Matteo Monelli<fnref n="1" />,²

Knowledge of the ages of the stars formed over a galaxy's lifetime is fundamental to an understanding of its formation and evolution. However, stellar ages are difficult to obtain since they cannot be measured from observations, but require comparison with stellar models<fnref n="1" />. Alternatively, age distributions can be derived by applying the robust technique of colour-magnitude diagram fitting², which until now has been used primarily to study nearby galaxies. Accurate distances to individual Milky Way stars now provided by the Gaia spacecraft mission³ have allowed us to derive ages from a thick-disk colour-magnitude diagram and from the two-sequenced colour-magnitude diagram of the kinematically hot local halo⁴, whose blue sequence has been linked to a major accretion event, Gaia-Enceladus⁵,⁶. Because accurate stellar ages were lacking, the time of the merger and its role in our Galaxy's early evolution remained unclear. Here we show that the stars in both halo sequences share identical age distributions, and are older than most of the thick-disk stars. The sharp halo age distribution cutoff at ten billion years ago can be identified with the time of accretion of Gaia-Enceladus to the Milky Way. Together with state-of-the-art cosmological simulations of galaxy formation⁷, these robust ages allow us to order the early sequence of events that shaped our Galaxy. We identify the red-sequence stars as the first stars formed within the Milky Way progenitor, and their kinematics indicate that these stars constitute the long-sought in situ halo of the Milky Way.

The accurate parallaxes and luminosities provided in the second data release (DR2) of the Gaia mission³ have allowed us to measure the intrinsic colours and magnitudes of stars located in a large volume of the Milky Way, encompassing different Galactic structural components. From these, we have produced colour-magnitude diagrams (CMDs) in units of absolute magnitudes and colours, which is what is required to derive star-formation histories² and stellar age distributions, by comparing them with theoretical CMDs derived from stellar evolution models⁸.

Fig. 1a,b show the CMD of two sub-populations of Milky Way stars taken from a parent population that lies within a sphere of about 2 kiloparsecs (kpc) around the Sun, as observed by the Gaia space observatory. Within this volume, accurate distances and absolute magnitudes can be derived directly from parallaxes. The CMD in Fig. 1a contains about 60,000 stars from this spherical region with large tangential velocities relative to the Sun (greater than 200 km s⁻<fnref n="1" />). Stars with such high velocities are classified in this study as belonging to a kinematically defined stellar halo<fnref n="1" />. The CMD in Fig. 1b contains some half a million stars from the same

spherical region but selected to be at least 1.1 kpc above or below the Galactic plane. At this distance from the plane, the majority of stars are expected to belong to the thick disk⁹, rather than to the young thin-disk component. We have excluded from this sample the stars with high velocity (greater than 200 km s⁻<fnref n="1" />) that have been included in the halo CMD. We note that our definition of the thick disk is morphological rather than chemical or kinematical: that is, some stars with thin-disk kinematics⁴ would be included in our sample. For simplicity, in this paper we will refer to the two samples we have just defined as 'halo' and 'thick disk' (see Methods for further information on the sample selection criteria). As is readily apparent, the CMD of the kinematically hot halo population contains two distinct sequences: a 'blue' sequence on the left, and a 'red' sequence on the right. This peculiar morphology indicates the presence of two distinct sub-populations within the halo. The blue sequence has been associated with a very important merging event—when the Gaia-Enceladus galaxy merged with the Milky Way—which occurred early during the formation of the Milky Way⁶, and of which the first hints date back to at least 15 years ago¹⁰⁻¹³. The nature of the red sequence has been less clear. It has been associated with the Milky Way's thick disk⁵, with support for this association coming from the characteristics of its chemical composition¹¹.

We compared the distribution of stars in these observed CMDs with that in CMDs computed with state-of-the-art stellar evolution models⁸ and found the combination of ages and metallicities (that is, the amount of chemical elements other than hydrogen and helium) of the stars that best reproduce the observed CMDs (see Methods). These best-fit CMDs are displayed in Fig. 1c,d. The model CMDs clearly succeed at reproducing the various features observed in the empirical CMDs.

Figure 2 displays the distribution of ages and metallicities of these best-fit model CMDs, divided into the blue and red sequences in the case of the halo CMD. The age distributions clearly demonstrate that the two sequences in the halo CMD are composed of stars that are coeval and formed at the earliest possible times in the life of the Universe (peak age 13.4 billion years (Gyr), with half of the stars already formed by 12.3 Gyr ago). The difference in colour is a result of the red- and blue-sequence stars having different metallicities. The models also clearly show that the thick-disk population spans a wider age range than the two halo populations: its stars closely follow the age distribution of the halo for approximately the first 2 Gyr, but the thick disk keeps forming stars after this initial period until about 6 Gyr ago, with its age distribution reaching a second peak later on (at about 9.5 Gyr ago, with half of the stars already formed by 10.5 Gyr ago).

<footnote n="1">¹Instituto de Astrofísica de Canarias, La Laguna, Spain. ²Departamento de Astrofísica, Universidad de La Laguna, La Laguna, Spain. ³Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France. ⁴INAF - Astronomical Observatory of Abruzzo, Teramo, Italy. ⁵INFN, Sezione di Pisa, Pisa, Italy. *e-mail: carme.gallart@iac.es</footnote>
<figure><img src="image_2.png" /></figure>

**Fig. 1** | Milky Way halo and thick-disk observed and modelled CMDs. a,b, The observed Gaia CMDs for the halo (a) and thick-disk (b) sub-populations. Colour is expressed on the x axis as the difference of the Gaia magnitudes in the BP and RP passbands. Absolute magnitudes in the Gaia G passband are indicated on the y axis. CMDs are represented as Hess diagrams with a logarithmic normalization of the number of stars in the colour-magnitude bins. Two sequences are visible in a: the blue sequence on the left and the red sequence on the right. c,d, The best-fit model CMDs that have been fit to the halo (c) and thick-disk (d) observed CMDs. The solutions for the fit using stellar evolution models with α-enhanced chemical composition are shown (Methods). The solutions obtained using stellar evolution models with solar-scaled chemical composition are not noticeably different. The dashed line in b and d indicates the position of the gap between the blue and red halo sequences. See Methods for more information on the sample selection and the CMD fitting technique.

These age results contrast with literature values for the ages¹²,¹⁴,¹⁵ assigned to stars associated with these components, on which previous conclusions about their nature relied⁶. These former age determinations indicated that the two halo populations had different ages, with blue sequence stars about 2 Gyr younger than stars in the red sequence. However, these ages were based either on (1) a relatively small number of stars, making use of spectroscopically determined stellar parameters that allow less-robust age estimates than absolute luminosities and colours or (2) indirect age estimates based on low-resolution spectroscopy. They contrasted with ages inferred with simple isochrone fitting taking into account average metallicities for stars in the two sequences⁴,⁶. Our robust age distributions for the thick-disk and halo components provide key insights into the first events that shaped the Milky Way's structure, providing evidence that complements previous studies, which largely dealt with chemical and kinematic information⁵,⁶,¹⁶⁻²¹.

That the two halo populations are coeval means not only that the red and blue populations formed stars at similar times, but

that they also stopped forming at similar times. Additionally, the well established relation²² between a galaxy's mass and the amount of metals contained in their stars means that the stars in the red sequence of the halo, being more metal-rich, must have formed in a galaxy that was more massive than the one where the stars in the blue sequence were formed. Both these findings are precisely what is expected if these two populations were involved in a merger event, with the red-sequence stars belonging to the main progenitor of the Milky Way, and the blue-sequence stars belonging to a smaller accreted galaxy, Gaia-Enceladus. This is shown in Fig. 3. Fig. 3a shows the age-metallicity distribution derived from the observed CMDs, with the red- and blue-sequence halo stars indicated as red and blue contours. Fig. 3b shows the signature in the age-metallicity plane left by a merger event in a simulated Milky Way analogue galaxy (see Methods), where the main progenitor stars are shown as red contours and those from the accreted galaxy as blue contours. The merger can be identified as two tracks in this plot, with the higher-metallicity stars belonging to the main
<figure><img src="image_3.png" /></figure>

**Fig. 2** | Milky Way halo and thick-disk stellar age and iron [Fe/H] distributions. **a,** The distribution of stellar ages of the thick disk (black histogram), and the blue (blue histogram) and red (red histogram) halo sequences according to their best-fit model CMDs, adopting the α-enhanced stellar evolution models. The solutions obtained using the solar-scaled models show similar features (see Supplementary Information). **b,** The distribution of [Fe/H] of the best-fitting models for the blue and red halo sequences as well as the thick-disk stars (solid lines), compared to the [Fe/H] distributions for the stars in each sample that had measured [Fe/H] in the LAMOST or GALAH spectroscopic surveys<sup>66-68</sup> (filled histograms). The fair agreement between the [Fe/H] distributions inferred with the CMD fitting technique and those measured spectroscopically provides an external check of the reliability of the results.

<figure><img src="image_4.png" /></figure>

**Fig. 3** | Age-metallicity relations for the two halo populations compared to a simulated Milky Way analogue. **a,** Age-metallicity relation for the best-fit model of the halo CMD. [M/H] on the y axis is used to represent the overall metallicity of the stars compared with that of the Sun. [M/H] = 0 corresponds to solar metallicity. The blue and red contours outline the distribution of stars in the blue and red halo sequences, respectively, overlaid on the general distribution of halo stars. **b,** As for **a** but for a simulated Milky Way analogue (see Methods). Blue and red contours depict the age-metallicity distribution of simulated stellar particles belonging to an accreted galaxy (analogous to Gaia-Enceladus) and the Galaxy progenitor, respectively. The last major merger can be identified as two tracks in the age-metallicity plot, with the higher-metallicity stars belonging to the main progenitor and the lower-metallicity stars coming from the merged satellite. The agreement, which we expected to be only qualitative, is remarkable. The time and mass of the last major merger, which can be seen to differ between the simulation and the real Milky Way, determine the details of the ages and metallicities of the red and blue halo star populations.

progenitor and the lower-metallicity stars coming from the less-massive merged galaxy.

The difference in the amount of metals in the red and blue sequences indicates that the accreted Gaia-Enceladus galaxy had about 30% (ref. ²³) of the mass of stars in the main progenitor of the Milky Way, although we stress that this ratio remains quite uncertain. This would indicate a total mass ratio of around 4:1 between the two galaxies, given the relation between stellar mass and total mass²⁴. This estimated mass ratio is similar to that derived in previous works⁶ using different arguments. Regardless of the exact mass ratio, this encounter heated some of the main progenitor stars that had been forming in a disk-like structure, to the extreme kinematics that led them to be classified as halo stars. The existence of such halo stars

has been predicted by cosmological simulations of Milky-Way-type galaxies¹⁷,²⁵ and has also been postulated from observational evidence⁵,¹¹,¹². Our age determination reveals that the stars heated to halo-like kinematics were among the first formed in the Milky Way, during approximately the first 3 Gyr of its evolution, just before the merger with Gaia-Enceladus took place. We can date the merger as occurring about 10 Gyr ago.

The similarity of the star-formation timescales for the two halo sequences also has important implications for the origin of their detailed chemical enrichment. It is known⁴⁻⁶,¹² that for stars with similar amounts of Fe, the blue-sequence stars (Gaia-Enceladus) have significantly less of the α-elements than do the red-sequence stars (formed in the main progenitor). The run of [α/Fe] versus [Fe]
<figure><img src="image_5.png" /></figure>

**Fig. 4** | Comparison of the stellar kinematics between observations of Milky Way stars and the simulated Milky Way analogue. **a,c,** The Toomre diagram for Milky Way stars (obtained by combining Gaia DR2 astrometry and LAMOST radial velocities$^{66}$), and for the simulated Milky Way analogue galaxy, respectively. $V_{Rv}$, $V_Z$ and $V_φ$ indicate radial, vertical and azimuthal velocities. **b,d,** Similar information as in **a** and **c** but presented as a $V_φ$ versus $V_Z$ diagram. In all panels, stars belonging to the observed halo blue sequence or the simulated accreted satellite are represented by blue shades, the observed halo red sequence or the simulated progenitor are represented as red contours, and the observed or simulated thick disk is shown in greyscale. We note that some stars in our halo sample turn out to have disk properties when the full kinematic information is taken into account. This may help to explain the presence of a small fraction of young stars in the halo sequences (Fig. 2a) as contamination from thin-disk stars in the sample. The same is true for the inferred thick-disk young population.

in different stellar systems is often used to trace the star-formation timescale of a galaxy, since $\alpha$ elements are produced in type II supernovae on shorter timescales than iron, which is mostly produced on a longer timescale by type Ia supernovae. The fact that we find similar star-formation timescales for the two components indicates that other mechanisms should be at work. In our simulations, the progenitor and satellite have different star-formation histories leading up to the merger$^{19}$, with the progenitor star-formation rate increasing while the satellite is small enough to have an almost constant star-formation rate, regulated by feedback.

A further interesting result concerns the relationship between the red-sequence (in situ) halo and the thick disk. The CMDs of these two components overlap, as seen in Fig. 1, as do their ages, chemical abundances and kinematics, as seen in Figs. 2 and 4 as well as in previous studies$^{4,5}$. The implication is that the in situ halo and thick disk have the same origin, linked to the main progenitor of the present-day Milky Way: in this self-consistent scenario, the old tail of the thick-disk population shown in Fig. 2 is populated by stars formed in the main progenitor, which were not heated to halo-like velocities. We may speculate that the metal-weak thick disk$^{26}$ would

be part of this population. Therefore, stars that would be currently associated with the thick disk were forming in the main progenitor before, during and after the merger with Gaia-Enceladus, whereas a fraction of main progenitor stars formed very early, and were heated to halo kinematics by the accretion event. The infalling gas associated with the merger probably contributed to maintaining and even stoking star formation in the early disk, the stellar age distribution of which shows a second peak of high star-formation rate about 9.5 Gyr ago. Such a formation scenario for the thick disk has been foreshadowed by cosmological galaxy-formation simulations$^{27}$ in which the thick disk forms at high redshift during a period characterized by gas-rich mergers.

A final implication of our findings regards the transition$^{28,29}$ between thick- and thin-disk formation. The derived ages for thick-disk stars indicate that this transition occurred around 8–6 Gyr ago, in agreement with previous studies$^{30}$, and that this transition was not associated with the earlier Gaia-Enceladus merger.

In conclusion, the reported new stellar age distributions, aided by state-of-the-art cosmological simulations of disk galaxy formation, delineate a clear picture of the formation of our Galaxy. In this
picture, a primitive Milky Way had been forming stars over a period of some 3 Gyr when a smaller galaxy, which had been forming stars on a similar timescale but was less chemically enriched owing to its lower mass, was accreted into it. This merger heated some of the existing stars in the main progenitor into a stellar halo-like configuration. A ready supply of infalling gas during the merger ensured the maintenance of a disk-like configuration, with the thick disk continuing to form stars at a substantial rate. Subsequently, around 8–6 Gyr ago, the gas settled into a thin disk that has continued to form stars until the present day.

## Methods

**Data selection.** A detailed and reliable determination of the age and metallicity distributions of a stellar system can be obtained from a CMD that is essentially complete down to the oldest main-sequence turn-offs¹. For a relatively metal-rich stellar population such as that of the Milky Way, the old turnoff occurs at an absolute magnitude $M_G = 4$ in the Gaia G band, but it is necessary to reach about one magnitude deeper to accurately sample this critical point in the CMD. At a distance of 2 kpc from the Sun, $M_G = 5$ corresponds to an apparent magnitude $G = 16.5$, down to which Gaia DR2 is essentially complete.

From Gaia DR2 (ref. ³) we selected all the stars inside a sphere of radius approximately 2 kpc centred on the Sun down to $M_{G,R} = 7$ (where $M_{G,R}$ is $M_G$ before applying any correction for interstellar extinction. $M_{G,R}$ was calculated directly from the apparent magnitude using the Gaia parallaxes, as $M_{G,R} = G + 5 + 5 \times \log_{10}(\omega/1,000.0)$, with $\omega$ being the parallax in milliarcseconds (mas). A parallax zero-point offset of 0.054 mas (ref. ³¹, in the sense that the corrected parallaxes are larger) has been applied. This value was preferred to the Gaia DR2 recommended value of 0.03 mas (ref. ³²) because several subsequent works have found a larger value³³⁻³⁸. The direct calculation of distance from the parallax value is valid because the great majority (over 95%) of stars in our sample have small relative parallax error (parallax over error > 5, ref. ³⁹). The stars excluded by this cutoff are more scattered in the CMD than the thick-disk and halo sequences. Even if they were to correspond to actual stars in these populations, their numbers are so low that neither the CMD fit or the derived age distribution would be affected. Finally, following refs. ⁴,³², it was also required that astrometric chi2 al/ (astrometric n good obs al - 5) < 1.44×max(1, exp(-0.4×(phot_g_mean_mag-19.5))) (only 2% of stars were removed by this cutoff) and that the stars had ($G_{BP} - G_{RP}$) measured colour.

From this parent population, we selected for the thick-disk sample all the stars in the northern and southern caps of the sphere, with distances from the Galactic plane $|Z| > 1.1$ kpc, and tangential velocities, calculated as $V_T = (4.74/\omega) \times [\mu_0^2 + \mu_\delta^2]$, where $\mu_0$ and $\mu_\delta$ are the proper motions in the right ascension and declination coordinates, respectively, lower than 200 km s⁻¹. For the halo sample, we selected all the relatively high Galactic latitude stars ($b > 30$) with $V_T > 20$ km s⁻¹. The latitude cut was performed to avoid substantial contamination from the thin disk at low Galactic latitude. We note that, since distances are actually obtained directly from parallax inversion, as the distance in parsec equal to 1,000 divided by the parallax in milliarcseconds, a small bias is introduced, such that true distances would tend to be smaller⁴⁰ by approximately 100 pc on average. This bias very slightly affects the exact spatial shape and size of the actual volume covered, but does not change the fact that we are considering a thick-disk sample (defined morphologically) and a halo sample (defined kinematically), and thus does not affect the conclusions of the paper.

Corrections for interstellar extinction were calculated on a star-by-star basis by interpolating the colour excess from a three-dimensional map¹¹. This map is strictly valid only up to heliocentric distances of 2.6 kpc and $|Z| < 0.6$ kpc, and thus, our thick-disk sample and part of the halo sample are outside its validity range. However, most of the dust is confined in a layer thinner than 0.6 kpc and consequently, the extinction to the thick-disk stars integrated up to this height is a very good approximation of the total extinction. We used published extinction coefficients¹², together with the effective temperature of each star (we neglected the dependence on metallicity). The temperature was estimated using a correlation between $G_{BP} - G_{RP}$ colour and temperature $T_{eff}$ determined from a large sample of Gaia DR2 stars with measured $T_{eff}$. Relatively small extinction values are inferred for most of the stars in both samples: 96% (or 93%) of the stars have $A_G < 0.2$ and 79% (or 70%) have $A_G < 0.1$ for the thick-disk (or halo) sample.

After retaining only stars brighter than $M_G = 5$, the thick-disk sample contains 518,968 stars and the halo sample contains 61,427 stars. The CMD of these stellar samples are displayed in Fig. 1a,b.

**Modelling of the CMD.** To derive the age and metallicity distributions of each population, the observed CMDs were modelled through the comparison with synthetic CMDs following well established techniques⁴³⁻⁵¹ that have been used for over 20 years to determine detailed star-formation histories for Local Group galaxies⁵²,⁵³. These galaxies are sufficiently close to observe individual stars, yet far enough away that all their stars can be considered to be at the same distance, which can be obtained accurately using various well calibrated distance indicators⁵³,⁵⁴.

However, applying this technique to populations of Milky Way stars has only been possible within an extremely small volume around the Sun⁵⁵⁻⁵⁸ owing to the difficulty of obtaining accurate distances for large stellar samples. This is one of the many aspects in which the Gaia mission is providing a breakthrough: in its DR2 (ref. ³), Gaia provided distances for 1.3 billion stars in the Milky Way. This distance information has now allowed us to apply the CMD-fitting technique to the different morphological components of our own Galaxy and therefore to robustly derive their age distributions. This provides an unprecedented opportunity to examine the early sequence of events that shaped the formation and structure of the Milky Way.

In the process of fitting an observed CMD, we seek the combination of simple stellar populations (that is, synthetic populations with small ranges of ages, $\Delta T \le 1$ Gyr, and metallicities, $\Delta[Fe/H] \approx 0.1$ dex) that provide the best fit. The minimization was performed with the code TheStorm²⁹,⁶⁰. This code uses a mother synthetic CMD from which the simple stellar populations are extracted. To check for possible systematic effects derived from the chemical properties of the stellar evolution models, these synthetic CMDs have been calculated using the BaSTI stellar evolution library in both its solar-scaled⁴¹ and $\alpha$-enhanced heavy-element distribution versions⁹. For the comparison itself, we make use of most of the area within the CMD populated by stars, excluding the Horizontal Branch, which is still difficult to reproduce quantitatively in stellar evolution models. Although this portion of the CMD is not included in the fit, we tried to obtain a fine match of this evolutionary stage in the observed CMDs of both the halo (where it is quite prominent) and the thick disk, where the blue extension of the Horizontal Branch is basically missing. Since, for any given metallicity, the Horizontal Branch morphology is mainly driven by the mass loss efficiency during the red-giant-branch stage, we used stellar model sets accounting for different efficiencies of this process as appropriate for each component: we used $\eta = 0.2$ for simulating the thick-disk CMD and $\eta = 0.4$ for the halo CMD, where $\eta$ is the free parameter used in combination with the Reimers' mass loss law in stellar modelling⁷. Although we chose different $\eta$ values for the different components under analysis, we should state that the results presented in this paper are not contingent upon these particular choices. Common to all mother synthetic CMDs (with $10^8$ stars) are a constant star-formation rate with flat age and metallicity distributions (14 Gyr to 30 Myr, and $-2.6 < [Fe/H] < -0.1$, for the halo, $-2.6 < [Fe/H] < 0.16$ for the thick disk) and the Kroupa initial mass function⁶². A binary fraction of 0.7 with mass ratio $q$ uniformly distributed between 0.1 and 1 was adopted. Finally, stellar evolutionary model predictions were transferred to the Gaia photometric system by using appropriate bolometric correction tables using the prescriptions for the DR2 catalogue⁶³.

The results displayed in Figs. 1–3 are obtained with the $\alpha$-enhanced models ($\alpha/[Fe] = 0.4$ (ref. ⁸), which are the most adequate for thick-disk and halo stellar populations. According to the analysis⁷ of the metallicity and $\alpha$ abundances¹¹ of a fraction of stars in the blue and red sequences of the halo CMD, red sequence stars have $\alpha/[Fe] \approx 0.3$ (with a mild dependence on metallicity), a value also typical for thick-disk stars⁴¹. Stars in the blue sequence have $\alpha/[Fe]$ values that range between $\alpha/[Fe] \approx 0.4$ at the lower-metallicity extreme and $\alpha/[Fe] \approx 0.15$ at the higher-metallicity extreme ($[Fe/H] \approx -0.8$ dex). However, we note that our conclusions, and in particular the age distributions, which are the most robust output of our CMD fitting process, are not affected by the choice of a scaled-solar heavy element distribution (Supplementary Fig. 1, equivalent to Fig. 2, is provided as Supplementary Information to demonstrate this point. It displays the results obtained adopting a solar-scaled heavy-element distribution). We expect, therefore, that the results would be even less affected by the small differences in the $\alpha/[Fe]$ that would be required to match exactly the chemical characteristics of the different populations.

Prior to the comparison between the observed Gaia CMDs and the mother synthetic CMDs, it is necessary to perform a simulation of the observational errors in the synthetic CMD, in such a way that it can be consistently compared with the observations. We considered two error sources, namely, the error in the absolute magnitudes introduced by the formal error in the parallax, and the error in the Gaia photometry. To model the first error source, which dominates the error in $M_G$, we have statistically reproduced the distributions of parallax and parallax over error of the observed sample in the synthetic CMD (the parallax over error was updated from the values provided in the Gaia DR2 archive taking into account the parallax offset). The theoretical $M_G$ values are then transformed to account for observational errors by adding a correction given by a Gaussian distribution centred at 0 and with a standard deviation given by $\Delta M_G = 2.17/\text{parallax over error}$. The maximum parallax errors in the samples used in this paper are of the order of 0.1 mas. Similarly, the errors in the Gaia photometry have been implemented by adding a dispersion of 0.004 magnitudes to the ($G_{BP} - G_{RP}$) colour of all synthetic stars at all magnitudes⁶³. We have not simulated uncertainties related to the extinction correction. These are expected to slightly smooth the CMD, and therefore the details in the star-formation history. These uncertainties are of the order of 0.02–0.03 mag in E(B–V) over most of the sky, and so one order or magnitude smaller than the dispersion of the CMD due to the parallax uncertainty (which vary between 0.0 and 0.3, with a typical value of 0.15).

The comparison between the observed and synthetic CMDs is performed using the number of stars in small colour–magnitude boxes. No a priori constraint on
the age-metallicity relation is adopted: The Storm code solves simultaneously for age and metallicity within the range covered by the simple stellar populations. For more details on the fitting procedure and error calculation, the reader is referred to refs.⁴⁹,⁶⁰,⁶⁵.

The results of this CMD fitting process are the star-formation rate as a function of time and the age-metallicity distribution. From these, it is possible to compute a best-fit CMD (displayed in Fig. 1c,d) which is then used to derive the age and metallicity distributions displayed in Figs. 2 and 3a. In the case of the halo CMD, the fitting process was performed in the CMD including both the blue and the red sequences. The age and metallicity distributions presented in Fig. 2 separately for each sequence were derived by drawing a polygon enclosing each sequence in the observed CMD (avoiding the red clump where both sequences naturally overlap) and using this polygon to assign synthetic stars in the best-fit CMD to either the blue or the red sequence. The boundary between the blue and red sequences corresponds to the fiducial line shown in Fig. 1. We have tested that this particular delimitation has a negligible impact on the results presented here.

**Cross-match with spectroscopic surveys.** We used a homogeneous catalogue⁶⁶ of distances, three-dimensional velocities, and ages for the cross-match between major Milky-Way spectroscopic surveys and Gaia DR2. We used this catalogue, together with the original Gaia DR2 information and the original [Fe/H] information from the LAMOST third data release A, F, G, K catalogue http://dr3.lamost.org/(fe b value)⁶⁷ and the GALAH first data release (fe h value)⁶⁸, to select stars matching exactly the same spatial and kinematical selection as our CMDs. The cross- identifications were performed requiring an exact match of the corresponding identification numbers in the different catalogues (Gaia source_id, LAMOST obsid and GALAH gaia_dr2_id). We selected these two surveys because (1) LAMOST includes a large number of stars on the main-sequence turnoff region of the CMD overlapping with our halo sample (2,131 stars after quality cutoffs; see below), and with our thick-disk sample (22,121 stars). SEGUE also has a significant overlap in the main-sequence region, but the survey target selection function somehow selected mostly stars on the blue main sequence. The other two high- resolution surveys in the catalogue (GES and APOGEE) contained too few stars in the main-sequence turnoff region to be used here. While GALAH offers only a small sample of stars overlapping with our halo sample (235 stars after the quality cutoff), it includes main-sequence stars covering both the blue and red sequences (125 and 82 stars, respectively), and it provides a larger overlap with the thick-disk sample (765 stars), made essentially of giant stars. The collected spectroscopic samples were used in Fig. 2 to build a metallicity distribution for the blue and red sequence halo, as well as the thick disk, and in Fig. 4 to explore the three- dimensional kinematics of the samples in cylindrical coordinates.

Some quality flags were used to clean the spectroscopic samples: both were required to have a quality flag = 0 in the compilation catalogue⁶⁶, as advised by the authors. Additionally, LAMOST stars were considered only when snrg≥30 and GALAH stars were considered only when flag_cannon=0.

**Galaxy simulation.** Cosmological simulations of the Milky Way's halo have generally been done by simulating halos of similar mass to the Milky Way (around $10^{12} M_{\odot}$), often imposing a relatively quiet merger history for the Milky Way since redshift $z \approx 1$. A range of increasingly detailed explorations of the origins of stars with high kinematics have been made using such simulations⁶⁹⁻⁷⁷.

For this study, we had a suite of 18 Milky Way mass galaxies from the NIHAO⁷⁸ and MaGICC⁷⁹ programmes, and selected the one which had a merger history and properties most like the Milky Way. The simulation used for Figs. 3 and 4 comes from a fully cosmological, hydrodynamical simulation of a Milky Way mass analogue galaxy from the MaGICC programme. Importantly, these simulated galaxies have been shown to match a wide variety of galaxy-scaling relations over a wide mass range, as well as reproducing a range of other galaxy properties such as enrichment of the circumgalactic medium⁷⁹. The particular realization, g15784, has a total mass of $1.4 \times 10^{12}$ solar masses and has been shown to be in good agreement with the Milky Way in terms of its detailed chemical and kinematic properties of the thick and thin disks⁸⁰. In the present study we show that the stars in the region around the Sun match the most important properties of the stars observed in our own Milky Way, and use this to interpret the red and blue sequences of the observed halo CMD as belonging to the main Milky Way progenitor and Gaia- Enceladus respectively, with the merger of these two proto-galaxies resulting in the in situ and vast majority of the inner stellar halo.

In the simulation, stars form when gas has become sufficiently cool ($T < 10,000$ K) and sufficiently dense ($n_{\text{th}} > 9.3 \text{ cm}^{-3}$). When gas particles meet these criteria, stars form according to: $dM/dt = c \cdot M_{\text{gas}}/t_{\text{dyn}}$, where M is the mass of stars formed in time dt, $M_{\text{gas}}$ is the mass of a gas particle, $t_{\text{dyn}}$ is the dynamical time of gas particles, and c is the star-formation efficiency—that is, the fraction of gas that will be converted into stars. Supernova feedback follows the blastwave model⁸¹ with thermal energy ($10^{51}$ erg) deposited into the surrounding interstellar medium from each supernova. The simulations also include radiation feedback from massive stars⁹². Cooling is disabled in the blast region (about 100 pc) for about 10 Myr. The simulation includes heating from a uniform ultraviolet ionizing background radiation field⁹³. Cooling takes into account both primordial gas and metals⁹⁴. The metal cooling grid is derived using CLOUDY⁹⁵.

Along with the fact that the same physical prescriptions used in this simulation resulted in matching scaling relations at a range of masses, meaning that low-mass accreted satellites contain the correct amounts of stars and gas for their halo masses, the correspondence with the observations of the thick disk gives us confidence in using this Milky Way analogue as an interpretive tool, as we have done here.

## Data availability

All data analysed in this paper are publicly available. They can be retrieved from the DR2 Gaia archive (http://gea.esac.esa.int/archive/) and from the compiled and supplemented GALAH and LAMOST catalogues⁶⁶, cross-matched as indicated in the Methods. The datasets containing information not available in public catalogues that is necessary to produce the figures can be retrieved from https://zenodo.org/record/3228143#.XRoQs3qH1mM (extended dataset hereafter). A explanatory README file is included. The extended dataset consists of the following files. (1) For Fig. 1a,b, two tables directly retrieved from the Gaia archive as described in the Methods, supplemented by extinction information on a star- by-star basis; the code used to interpolate the three-dimensional extinction maps¹¹ can be retrieved from https://github.com/edober/dust3d. (2) For Figs. 1c,d, 2 and 3a, two tables with the derived solution CMDs and two files with the data needed to define the boxes used to select stars in the blue and red sequences of the halo CMD. (3) For Fig. 3b, the necessary tables with the age, metallicity and velocity data for the main progenitor and accreted satellite for realization g15784 of the MaGICC program. The complete information on the final timestep of that simulation, together with scripts to read and plot the data are also included in the extended dataset.

Received: 26 April 2019; Accepted: 24 May 2019;
Published online: 22 July 2019

## References

1. Soderblom, D. R. The ages of stars. Annu. Rev. Astron. Astrophys. **48**, 581–629 (2010).

2. Gallart, C., Zoccali, M. & Aparicio, A. The adequacy of stellar evolution models for the interpretation of the color-magnitude diagrams of resolved stellar populations. Annu. Rev. Astron. Astrophys. **43**, 387–434 (2005).

3. Gaia Collaboration Gaia data release 2. Summary of the contents and survey properties. Astron. Astrophys. **616**, A1 (2018).

4. Gaia Collaboration et al. Gaia data release 2. Observational Hertzsprung-Russell diagrams. Astron. Astrophys. **616**, A10 (2018).

5. Haywood, M. et al. In disguise or out of reach: first clues about in situ and accreted stars in the stellar halo of the Milky Way from Gaia DR2. Astrophys. J. **863**, 113 (2018).

6. Helmi, A. et al. The merger that led to the formation of the Milky Way's inner stellar halo and thick disk. Nature **563**, 85–88 (2018).

7. Brook, C. B., Stinson, G., Gibson, B. K., Wadsley, J. & Quinn, T. MaGICC discs: matching observed galaxy relationships over a wide stellar mass range. Mon. Not. R. Astron. Soc. **424**, 1275–1283 (2012).

8. Pietrinferni, A., Cassisi, S., Salaris, M. & Castelli, F. A large stellar evolution database for population synthesis studies. II. Stellar models and isochrones for an α-enhanced metal distribution. Astrophys. J. **642**, 797–812 (2006).

9. Gilmore, G. & Reid, N. New light on faint stars. III. Galactic structure towards the South Pole and the Galactic thick disc. Mon. Not. R. Astron. Soc. **202**, 1025–1047 (1983).

10. Brook, C. B., Kawata, D., Gibson, B. K. & Flynn, C. Galactic halo stars in phase space: a hint of satellite accretion? Astrophys. J. **585**, L125–L129 (2003).

11. Nissen, P. E. & Schuster, W. J. Two distinct halo populations in the solar neighborhood. Evidence from stellar abundance ratios and kinematics. Astron. Astrophys. **511**, L10 (2010).

12. Schuster, W. J., Moreno, E., Nissen, P. E. & Pichardo, B. Two distinct halo populations in the solar neighborhood. III. Evidence from stellar ages and orbital parameters. Astron. Astrophys. **538**, A21 (2012).

13. Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E. & Deason, A. J. Co-formation of the disc and the stellar halo. Mon. Not. R. Astron. Soc. **478**, 611–619 (2018).

14. Hawkins, K., Jofré, P., Gilmore, G. & Masseron, T. On the relative ages of the α-rich and α-poor stellar populations in the Galactic halo. Mon. Not. R. Astron. Soc. **445**, 2575–2588 (2014).

15. Ge, Z. S. et al. Ages of 70 dwarfs of three populations in the solar neighborhood: considering O and C abundances in stellar models. Astrophys. J. **833**, 161 (2016).

16. Hawkins, K., Jofré, P., Masseron, T. & Gilmore, G. Using chemical tagging to redefine the interface of the Galactic disc and halo. Mon. Not. R. Astron. Soc. **453**, 758–774 (2015).

17. Bonaca, A., Conroy, C., Wetzel, A., Hopkins, P. F. & Kereš, D. Gaia reveals a metal-rich, in situ component of the local stellar halo. Astrophys. J. **845**, 101 (2017).

18. Hayes, C. R. et al. Disentangling the Galactic halo with APOGEE. I. Chemical and kinematical investigation of distinct metal-poor populations. Astrophys. J. **852**, 49 (2018).
19. Fernández-Alvar, E. et al. Disentangling the Galactic halo with APOGEE. II. Chemical and star formation histories for the two distinct populations. *Astrophys. J.* **852**, 50 (2018).

20. Mackereth, J. T. et al. The origin of accreted stellar halo populations in the Milky Way using APOGEE, Gaia, and the EAGLE simulations. *Mon. Not. R. Astron. Soc.* **482**, 3426–3442 (2019).

21. Di Matteo, P. et al. The Milky Way has no in-situ halo but it has a thick disc. Composition of the stellar halo and age-dating the last significant merger with Gaia DR2 and APOGEE. Preprint at https://arxiv.org/abs/1812.08232 (2018).

22. Erb, D. K. et al. The mass–metallicity relation at $z > 2$. *Astrophys. J.* **644**, 813–828 (2006).

23. Ma, X. et al. The origin and evolution of the galaxy mass-metallicity relation. *Mon. Not. R. Astron. Soc.* **456**, 2140–2156 (2016).

24. Behroozi, P. S., Wechsler, R. H. & Conroy, C. The average star formation histories of galaxies in dark matter halos from $z = 0–8$. *Astrophys. J.* **770**, 57 (2013).

25. Zolotov, A. et al. The dual origin of stellar halos. *Astrophys. J.* **702**, 1058–1067 (2009).

26. Beers, T. C. et al. Population studies. XIII. A new analysis of the Bidelman–MacConnell ‘weak-metal’ stars confirmation of metal-poor stars in the thick disk of the galaxy. *Astrophys. J.* **794**, 58 (2014).

27. Brook, C. B., Kawata, D., Gibson, B. K. & Freeman, K. C. The emergence of the thick disk in a cold dark matter universe. *Astrophys. J.* **612**, 894–899 (2004).

28. Fuhrmann, K. Nearby stars of the Galactic disc and halo—V. *Mon. Not. R. Astron. Soc.* **414**, 2893–2922 (2011).

29. Recio-Blanco, A. et al. The Gaia-ESO survey: the Galactic thick to thin disc transition. *Astrophys. Rev. Astron.* **567**, A5 (2014).

30. Haywood, M. et al. When the Milky Way turned off the lights: APOGEE provides evidence of star formation quenching in our Galaxy. *Astron. Astrophys.* **589**, A66 (2016).

31. Schönrich, R., McMillan, P. & Eyer, L. Distances and parallax bias in Gaia DR2. Preprint at https://arxiv.org/abs/1902.02355 (2019).

32. Lindegren, L. et al. Gaia data release 2. The astrometric solution. *Astron. Astrophys.* **616**, A2 (2018).

33. Graczyk, D. et al. Testing systematics of Gaia DR2 parallaxes with empirical surface brightness: color relations applied to eclipsing binaries. *Astrophys. J.* **872**, 85 (2019).

34. Hall, O. J. et al. Testing asteroseismology with Gaia DR2: hierarchical models of the red clump. *Mon. Not. R. Astron. Soc.* **486**, 3569–3585 (2019).

35. Khan, S. et al. New light on the Gaia DR2 parallax zero-point: influence of the asteroseismic approach, in and beyond the Kepler field. Preprint at https://arxiv.org/abs/1904.05676 (2019).

36. Leung, H. W. & Bovy, J. Simultaneous calibration of spectro-photometric distances and the Gaia DR2 parallax zero-point offset with deep learning. Preprint at https://arxiv.org/abs/1902.08634 (2019).

37. Muraveva, T., Delgado, H. E., Clementini, G., Sarro, L. M. & Garofalo, A. RR Lyrae stars as standard candles in the Gaia data release 2 era. *Mon. Not. R. Astron. Soc.* **481**, 1195–1211 (2018).

38. Zinn, J. C., Pinsonneault, M. H., Huber, D. & Stello, D. Confirmation of the Gaia DR2 parallax zero-point offset using asteroseismology and spectroscopy in the Kepler field. Preprint at https://arxiv.org/abs/1805.02650 (2018).

39. Luri, X. et al. Gaia data release 2. Using Gaia parallaxes. *Astron. Astrophys.* **616**, A9 (2018).

40. Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G. & Andrae, R. Estimating distance from parallaxes. IV. Distances to 1.33 billion stars in Gaia data release 2. *Astron. J.* **156**, 58 (2018).

41. Lallement, R. et al. Three-dimensional maps of interstellar dust in the local arm: using Gaia, 2MASS, and APOGEE-DR14. *Astron. Astrophys.* **616**, A132 (2018).

42. Casagrande, L. & VandenBerg, D. A. On the use of Gaia magnitudes and new tables of bolometric corrections. *Mon. Not. R. Astron. Soc.* **479**, L102–L107 (2018).

43. Gallart, C., Freedman, W. L., Aparicio, A., Bertelli, G. & Chiosi, C. The star formation history of the Local Group dwarf galaxy Leo I. *Astron. J.* **118**, 2245–2261 (1999).

44. Hernandez, X., Valls-Gabaud, D. & Gilmore, G. Deriving star formation histories: inverting Hertzsprung–Russell diagrams through a variational calculus maximum likelihood method. *Mon. Not. R. Astron. Soc.* **304**, 705–719 (1999).

45. Holtzman, J. A. et al. Observations and implications of the star formation history of the Large Magellanic Cloud. *Astron. J.* **118**, 2262–2279 (1999).

46. Dolphin, A. E. Numerical methods of star formation history measurement and applications to seven dwarf spheroidals. *Mon. Not. R. Astron. Soc.* **332**, 91–108 (2002).

47. Aparicio, A. & Gallart, C. IAC-STAR: a code for synthetic color-magnitude diagram computation. *Astron. J.* **128**, 1465–1477 (2004).

48. Aparicio, A. & Hidalgo, S. L. IAC-pop: finding the star formation history of resolved galaxies. *Astron. J.* **138**, 558–567 (2009).

49. Monelli, M. et al. The ACS LCID project. III. The star formation history of the Cetus dSph galaxy: a post-reionization fossil. *Astrophys. J.* **720**, 1225–1245 (2010).

50. Cignoni, M. & Tosi, M. Star formation histories of dwarf galaxies from the colour-magnitude diagrams of their resolved stellar populations. *Adv. Astron.* **2010**, 158568 (2010).

51. Ruiz-Lara, T. et al. Integrated-light analyses vs. colour-magnitude diagrams. II. Leo A: an extremely young dwarf in the Local Group. *Astron. Astrophys.* **617**, A18 (2018).

52. Tolstoy, E., Hill, V. & Tosi, M. Star-formation histories, abundances, and kinematics of dwarf galaxies in the Local Group. *Annu. Rev. Astron. Astrophys.* **47**, 371–425 (2009).

53. Benedict, G. F. et al. Hubble Space Telescope fine guidance sensor parallaxes of galactic Cepheid variable stars: period–luminosity relations. *Astron. J.* **133**, 1810–1827 (2007).

54. Beaton, R. L. et al. The Carnegie-Chicago Hubble Program. I. An independent approach to the extragalactic distance scale using only population II distance indicators. *Astrophys. J.* **832**, 210 (2016).

55. Hernandez, X., Valls-Gabaud, D. & Gilmore, G. The recent star formation history of the Hipparcos solar neighbourhood. *Mon. Not. R. Astron. Soc.* **316**, 605–612 (2000).

56. Bertelli, G. & Nasi, E. Star formation history in the solar vicinity. *Astron. J.* **121**, 1013–1023 (2001).

57. Cignoni, M., Degl'Innocenti, S., Prada Moroni, P. G. & Shore, S. N. Recovering the star formation rate in the solar neighborhood. *Astron. Astrophys.* **459**, 783–796 (2006).

58. Bernard, E. J. Reconstructing the star formation history of the solar neighbourhood with Gaia. *Proc. Int. Astron. Union* **334**, 158–161 (2018).

59. Bernard, E. J. et al. The spatially-resolved star formation history of the M31 outer disc. *Mon. Not. R. Astron. Soc.* **453**, L113–L117 (2015).

60. Bernard, E. J. et al. Star formation history of the Galactic bulge from deep HST imaging of low reddening windows. *Mon. Not. R. Astron. Soc.* **477**, 3507–3519 (2018).

61. Pietrinferni, A., Cassisi, S., Salaris, M. & Castelli, F. A large stellar evolution database for population synthesis studies. I. Scaled solar models and isochrones. *Astrophys. J.* **612**, 168–190 (2004).

62. Kroupa, P., Tout, C. A. & Gilmore, G. The distribution of low-mass stars in the Galactic disc. *Mon. Not. R. Astron. Soc.* **262**, 545–587 (1993).

63. Evans, D. W. et al. Gaia data release 2. Photometric content and validation. *Astron. Astrophys.* **616**, A4 (2018).

64. Hayden, M. R. et al. Chemical cartography with APOGEE: metallicity distribution functions and the chemical structure of the Milky Way disk. *Astrophys. J.* **808**, 132 (2015).

65. Hidalgo, S. L. et al. The ACS LCID project. V. The star formation history of the dwarf galaxy LGS-3: clues to cosmic reionization and feedback. *Astrophys. J.* **730**, 14 (2011).

66. Sanders, J. L. & Das, P. Isochrone ages for 3 million stars with second Gaia data release. *Mon. Not. R. Astron. Soc.* **481**, 4093–4110 (2018).

67. Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P. & Deng, L.-C. LAMOST spectral survey—an overview. *Res. Astron. Astrophys.* **12**, 723–734 (2012).

68. Buder, S. et al. The GALAH survey: second data release. *Mon. Not. R. Astron. Soc.* **478**, 4513–4552 (2018).

69. Brook, C. B., Kawata, D., Gibson, B. K. & Flynn, C. Stellar halo constraints on simulated late-type galaxies. *Mon. Not. R. Astron. Soc.* **349**, 52–56 (2004).

70. Okamoto, T., Eke, V. R., Frenk, C. S. & Jenkins, A. Effects of feedback on the morphology of galaxy discs. *Mon. Not. R. Astron. Soc.* **363**, 1299–1314 (2005).

71. Font, A. S. et al. Cosmological simulations of the formation of the stellar haloes around disc galaxies. *Mon. Not. R. Astron. Soc.* **416**, 2802–2820 (2011).

72. Pillepich, A. et al. Halo mass and assembly history exposed in the faint outskirts: the stellar and dark matter haloes of Illustris galaxies. *Mon. Not. R. Astron. Soc.* **444**, 237–249 (2014).

73. Rodriguez-Gomez, V. et al. The stellar mass assembly of galaxies in the Illustris simulation: growth by mergers and the spatial distribution of accreted stars. *Mon. Not. R. Astron. Soc.* **458**, 2371–2390 (2016).

74. Tissera, P. B. et al. The central spheroids of Milky Way mass-sized galaxies. *Mon. Not. R. Astron. Soc.* **473**, 1656–1666 (2018).

75. Obreja, A. et al. Introducing galactic structure finder: the multiple stellar kinematic structures of a simulated Milky Way mass galaxy. *Mon. Not. R. Astron. Soc.* **477**, 4915–4930 (2018).

76. Monachesi, A. et al. The Auriga stellar haloes: connecting stellar population properties with accretion and merging history. *Mon. Not. R. Astron. Soc.* **485**, 2589–2616 (2019).

77. Fattahi, A. et al. The origin of galactic metal-rich stellar halo components with highly eccentric orbits. *Mon. Not. R. Astron. Soc.* **484**, 4471–4483 (2019).

78. Wang, L. et al. NIHAO project—I. Reproducing the inefficiency of galaxy formation across cosmic time with a large sample of cosmological hydrodynamical simulations. *Mon. Not. R. Astron. Soc.* **454**, 83–94 (2015).
79. Stinson, G. S. et al. MAGICC haloes: confronting simulations with observations of the circumgalactic medium at $z = 0$. *Mon. Not. R. Astron. Soc.* **425**, 1270–1277 (2012).

80. Miranda, M. S. et al. Origin of the metallicity distribution in the thick disc. *Astron. Astrophys.* **587**, A10 (2016).

81. Stinson, G. et al. Star formation and feedback in smoothed particle hydrodynamic simulations—I. Isolated galaxies. *Mon. Not. R. Astron. Soc.* **373**, 1074–1090 (2006).

82. Stinson, G. S. et al. Making galaxies in a cosmological context: the need for early stellar feedback. *Mon. Not. R. Astron. Soc.* **428**, 129–140 (2013).

83. Haardt, F. & Madau, P. Radiative transfer in a clumpy universe. II. The ultraviolet extragalactic background. *Astrophys. J.* **461**, 20 (1996).

84. Shen, S., Wadsley, J. & Stinson, G. The enrichment of the intergalactic medium with adiabatic feedback—I. Metal cooling and metal diffusion. *Mon. Not. R. Astron. Soc.* **407**, 1581–1596 (2010).

85. Ferland, G. J. et al. CLOUDY 90: numerical simulation of plasmas and their spectra. *Publ. Astron. Soc. Pac.* **110**, 761–778 (1998).

## Acknowledgements

C.G., C.B.B., M.M. and T.R.-L. acknowledge financial support through the grants (AEI/FEDER, UE) AYA2017-89076-P, AYA2016-77237-C3-1-P and AYA2015-63810-P, as well as by the Ministerio de Ciencia, Innovación y Universidades (MCIU), through the State Budget and by the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community, through the Regional Budget. T.R.-L. is supported by a MCIU Juan de la Cierva - Formación grant (FJCI-2016-30342). C.B.B. is supported by a MCIU Ramón y Cajal Fellowship (RYC 2013-12784). S.C. acknowledges support from Premiale INAF 'MITIC' and has been supported by INFN (Iniziativa specifica TAsP). V.H. was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. We used data from the European Space Agency mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; see http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We also used data from the LAMOST and GALAH surveys. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project was provided by the National Development and Reform Commission.

LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. The GALAH survey is based on observations made at the Australian Astronomical Observatory, under programmes A/2013B/13, A/2014A/25, A/2015A/19 and A/2017A/18. We acknowledge the traditional owners of the land on which the Australian Astronomical Observatory stands, the Gamilaraay people, and pay our respects to elders past and present.

## Author contributions

All authors critically contributed to different aspects of the data analysis and model calculation, and to the interpretation of the results. The writing of the manuscript, to which all authors contributed, was led by C.G. and C.B.B. Specifically, C.G. selected the Gaia samples and performed the CMD fitting, E.J.B. helped with the data selection, wrote the CMD fitting software (TheStorm) based on earlier work with M.M. and C.G., and performed the calculation of the three-dimensional interstellar reddening. C.B.B. contributed the galaxy formation models, which were key to the interpretation of the results. T.R.-L. participated in the CMD fitting, contributed key software for various steps such as the error simulation in the synthetic CMDs, and created the figures. S.C. contributed the software to calculate the synthetic CMD, including all the necessary libraries of stellar models and bolometric correction tables for the Gaia photometric passbands. V.H. selected the spectroscopic samples, which were analysed together with T.R.-L.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41550-019-0829-5.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Correspondence and requests for materials** should be addressed to C.G.

**Peer review information:** *Nature Astronomy* thanks Carine Babusiaux and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019