# The recurrent impact of the Sagittarius dwarf on the star formation history of the Milky Way

Tomás Ruiz-Lara, Carme Gallart, Edouard J. Bernard, Santi Cassisi

► To cite this version:

Tomás Ruiz-Lara, Carme Gallart, Edouard J. Bernard, Santi Cassisi. The recurrent impact of the Sagittarius dwarf on the star formation history of the Milky Way. Nature Astronomy, 2020, 4, pp.965-973. 10.1038/s41550-020-1097-0. insu-03689262

HAL Id: insu-03689262

https://insu.hal.science/insu-03689262v1

Submitted on 2 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
<figure><img src="image_2.png" /></figure>

# The recurrent impact of the Sagittarius dwarf on the star formation history of the Milky Way

Tomás Ruiz-Lara<fnref n="1" />,²✉, Carme Gallart<fnref n="1" />,², Edouard J. Bernard<fnref n="1" />,³ and Santi Cassisi⁴,⁵

Satellites orbiting disk galaxies can induce phase space features such as spirality, vertical heating and phase-mixing in their disks. Such features have also been observed in our own Galaxy, but the complexity of the Milky Way disk has only recently been fully mapped by Gaia Data Release 2 (DR2) data. This complex behaviour is mainly ascribed to repeated perturbations induced by the Sagittarius dwarf galaxy (Sgr) along its orbit, pointing to this satellite as the main dynamical architect of the Milky Way disk. Here, we model Gaia DR2-observed colour-magnitude diagrams to obtain a detailed star formation history of the ~2 kpc bubble around the Sun. It reveals three conspicuous and narrow episodes of enhanced star formation that we can precisely date as having occurred 5.7, 1.9 and 1.0 Gyr ago. The timing of these episodes coincides with proposed Sgr pericentre passages according to (1) orbit simulations, (2) phase space features in the Galactic disk and (3) Sgr stellar content. These findings most probably suggest that Sgr has also been an important actor in the build-up of the stellar mass of the Milky Way disk, with the perturbations from Sgr repeatedly triggering major episodes of star formation.

The accurate photometry and parallaxes delivered by Gaia DR2¹, when combined with well-established techniques of modelling observed colour-magnitude diagrams² (CMD, see Fig. 1), enable the computation³ of detailed star formation histories (SFH) for large volumes within the Milky Way. In this work we use such techniques, applied to Gaia DR2 data, to characterize the stellar content within a bubble around the Sun with a radius of ~2 kpc. Figure 2 shows the derived SFH. We find that the Milky Way has been forming stars continuously throughout its entire history, with a decreasing intensity compatible with other massive spirals⁴. On top of this overall behaviour, it displays three strikingly well-defined star formation enhancements that occurred ~5.7, 1.9 and 1.0 Gyr ago, with decreasing strength. There is also a hint of a fourth possible star formation burst spanning the last 70 Myr. Extensive testing with mock stellar populations on the ability of the method to recover star formation bursts, both in isolation and in combination with an exponential decline, confirms their detection with a significance above 3–4σ, and a maximum duration of 2.0, 0.6 and 0.4 Gyr, respectively (see Fig. 3 and Methods for a description of the tests). The episodic nature of the recovered SFH for the Milky Way disk is remarkable and although hints of star forming bursts in the Milky Way disk at specific times have been suggested before, they have been detected and characterized with precision and significance in the present work.

Previous determinations of the SFH of the Milky Way disk, using Hipparcos⁵ data⁶⁻⁸ or different methodologies³,⁹, exist. However, they either concentrate on the solar neighbourhood within 250 pc and therefore provide a local view of the SFH of our Galaxy disk, or lack the age resolution needed to isolate individual star formation events. The combination of the large volume probed by the Gaia DR2 data and our state-of-the-art methodology allowed us to address both limitations. Such previous works show a relatively constant SFR with a clear enhancement in approximately the past 4 Gyr (refs. ³,⁷,⁸,¹⁰). The differences between our derived SFH and that of previous works are considerably minimized when our analysis is

restricted to a similar region in the solar neighbourhood or when our age resolution is reduced. However, while previous works are unable to detect an enhanced star formation around 5–6 Gyr ago, this is a robust feature of our SFH. We conclude that it is the combination of improved data (Gaia DR2) and an improved methodology that allows us to clearly detect this elusive (but important) star forming burst, previously barely hinted at⁹. The age-metallicity relation derived in this work is also compatible with that found from the study of individual stars within the Milky Way¹¹,¹².

The analysed volume (a bubble of ~2kpc in radius) necessarily comprises stars of very different chemical and dynamical properties¹³,¹⁴. To understand how the SFH varies for different kinematic components of the Milky Way disk, we divide our original sample into thin- and thick-disk stars on the basis of their kinematics and following the same criterion used in ref. ¹⁵, based on the tangential velocities of individual stars. Figure 4 shows the recovered SFHs for the thin (light blue) and thick (dark blue) disks, respectively. As expected, the thick disk presents a larger contribution of old stars than the thin disk and, although the youngest bursts are also visible in this component, they are clearly less prominent than in the thin disk. However, the burst at ~5.7 Gyr is prominent in both kinematic components (16% of the thick stars, 24% of the thin disk stars), suggesting that the event giving rise to such burst should have been more powerful and have a more global effect than those producing the more recent bursts (mainly observed in the thin disk). Pinpointing the causes of these star formation episodes will boost our knowledge of our own Galaxy history, providing crucial information on the processes driving star formation on galactic scales.

Galaxy major mergers are thought to be one of the main factors triggering star formation in galaxies¹⁶⁻¹⁸, and are predicted to have clear effects on their chemical evolution¹⁹. In particular, the current cosmological paradigm predicts that massive mergers play a crucial role in the formation of galaxies²⁰, including the Milky Way²¹,²². However, we do not have any indications of such massive mergers happening in the late evolution of our Galaxy²⁰,²³.

<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, OCA, CNRS, Lagrange, France. ⁴INAF-Astronomical Observatory of Abruzzo, Teramo, Italy. ⁵INFN - Physics Department, University of Pisa, Pisa, Italy.
✉e-mail: tomasruizlara@gmail.com</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | CMD of the - 2-kpc bubble around the Sun. **a**, Observed $M_G$ versus $G_{BP} - G_{RP}$ Gaia DR2 CMD using dereddened absolute magnitudes. Colour is expressed on the x axis as the difference of the Gaia magnitudes in the BP ($G_{BP}$) and RP ($G_{RP}$) passbands. Absolute magnitudes in the Gaia G passband are indicated on the y axis ($M_G$). Lines are representative of the average run in the colour-magnitude plane of theoretical stars belonging to each star formation enhancement. **b**, The best-fitting model to the observed CMD. In both panels we use a Hess diagram representation of the CMDs with the number of stars colour coded on a logarithmic scale. **c**, Position in the CMD of synthetic stars linked to each of the four (including also the youngest one) star formation enhancements using different hues of purple. The bundle strategy chosen for the CMD fitting is overlaid in black. The sizes of the boxes (ΔColour - ΔMagnitude) in each bundle are given in the inset table. Note that the highest weighting applied in the fitting corresponds to the more densely sampled bundle 2. **d**, Residual difference between the observed and the best-fit model CMDs in Poissonian $\sigma$ (colour scale). Note that in the great majority of the CMD, the residuals are lower than $2\sigma$. **e,f**, $M_G$ (e) and $G_{BP} - G_{RP}$ (f) distributions of the observed and best-fit model stars, respectively.</figcaption></figure>

On the other hand, Milky Way halo substructures or stellar streams such as those linked to Sgr$^{24,25}$ do suggest that our Galaxy has experienced minor satellite encounters in the past few billion years. We compared the times at which we find star formation enhancements with those of the proposed first infall and following pericentre passages for Sgr from numerical simulations. The coincidence is striking (see Fig. 2c). Mimicking the observed angular positions, distances and radial velocities of tidal streams from Sgr, simulations of Sgr tidal disruption suggest close pericentric passages (~10kpc) around 6.5, 4.5, 2.75, 1, and 0.1 Gyr ago$^{26}$ (secondary pericentres at 5.5, 3.7, and 1.9 Gyr ago). Other attempts to reconstruct the orbital history of Sgr, some including the possible effect of the Magellanic Clouds, suggest that Sgr could have

been affecting the Milky Way disk since 6–7 Gyr ago, dating the first Sgr infall to around 5 Gyr ago, with subsequent passages also compatible with the observed SFH enhancements$^{27,28}$. Further constraints on the orbit of Sgr come from the observed peculiarities in the phase space and stellar streaming motions in the Milky Way disk$^{29–31}$, which have been explained as the result of perturbations excited by pericentre passages of Sgr$^{32}$. Such observations can be explained if Sgr has a total mass of $M_{\text{total}} \approx 2.5 \times 10^{10} M_{\odot}$ and passed close to the Milky Way disk 2.2 and 1.1 Gyr ago. These two recent pericentric passages also agree closely with two of the bursts of star formation derived with our techniques. New studies based on Gaia DR2 data have unveiled even more detailed kinematic features that suggest that the Galaxy is experiencing an ongoing phase-mixing
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | SFH representative of the ~2-kpc-radius bubble around the Sun. **a**, Normalized star formation rate (SFR) as a function of look-back time for the different tests under analysis. **b**, The youngest star formation enhancement spanning approximately the last 100 Myr. Solid and dashed purple lines display the results from the different tests using two different extinction coefficients. The solid orange line depicts the overall star formation history averaging all the tests (see Methods). Shaded areas represent the errors in the SFH determinations, and orange error bars represent the errors in the peak significance and age resolution, respectively. These error bars do not account for possible degeneracies in the solution or systematic errors induced by the choice of input parameters of our method (see Methods). Green areas highlight the recovered star formation enhancements. **c**, Comparison between our results and knowledge regarding the orbit of Sgr and stellar content from literature. Asterisks and line intervals correspond to proposed pericentric passages (grey asterisks mark secondary pericentres). The arrow indicates the period when interactions have been traced back according to ref. ²⁸ (6–7 Gyr). Lines in greyscale mimic variable star formation from black (star forming) to white (no star formation) for the Sgr main body³⁸ or stream³⁷.</figcaption></figure>

process<sup>14,33,34</sup>. Comparison of the observed distribution of stars in the vertical position and velocity (Z–V<sub>Z</sub>) plane with N-body models of galaxy dynamics suggests that the Milky Way disk was perturbed just 300–900 Myr ago, in reasonable agreement with the most recent Sgr pericentre passage<sup>26–28,35</sup>. Such findings reinforce the idea of Sgr being the main perturbing agent during the last 6 Gyr of the evolution of the Milky Way<sup>34</sup>. In fact, a recent work combining APOGEE data with Gaia DR2 claims that satellite perturbations are an important heating mechanism of the outer disk that could explain the high vertical velocity dispersion of ~6-Gyr-old stars as well as the relationship between [Fe/H] and the radial age–velocity dispersion<sup>36</sup>. According to our derived SFH, Sgr is also an important actor in the build-up of the stellar content of the Milky Way disk.

So far, we have discussed only the effect of Sgr on the Milky Way. However, such interactions should also have had an effect on Sgr, which is considerably less massive than our Galaxy<sup>27</sup>. Analyses of the stellar content of the Sgr stream via CMD fitting infer a typical dwarf-like star formation and chemical enrichment history until 5–7 Gyr ago<sup>37</sup>. At that time, star formation in the Sgr stream halted, giving way to a residual activity. This drop in star formation, together with the presence of the first

enhancement on the SFH reported in this work, strongly suggest that the first infall of Sgr into the Milky Way potential took place around 5.5 Gyr ago. The event forced the stripping of gas in Sgr, quenching star formation, and triggered the formation of new stars in the Milky Way due to induced shocks and perturbations and, possibly, gas accretion. However, analysis of just the stream, associated with the outer parts of the original dwarf galaxy and most probably stripped during this first infall, prevents the reconstruction of more recent events. Fortunately, high-quality Hubble Space Telescope data of the core of the Sgr main body are also available<sup>38</sup>. Its CMD is characterized by the presence of several main-sequence turn-offs (the most robust age diagnostics), indicating distinct star forming epochs. Apart from a metal-poor population older than 8 Gyr (globular cluster-like), careful examination of the stellar content of the Sgr core is indicative of three distinct stellar populations: (1) a prominent and extended intermediate-age population (4–8 Gyr old), with the bulk of stars being around 5–6 Gyr old; (2) a younger, slightly more metal-rich population 2.3 Gyr old; and (3) a metal-rich (above solar), young population ~0.9 Gyr old. The tight correlation between the stellar content of Sgr and the Milky Way disk SFH further supports the
<figure><img src="image_5.png" /></figure>

**Fig. 3** | Testing the robustness of the SFH recovery with mock stellar populations. The recovered SFH for the 2-kpc volume analysed from Gaia, the input SFH for the mock stellar populations and the recovered SFH after mimicking errors in the associated CMD and applying our methodology are shown. The shaded grey envelope indicates errors in the SFH determinations (see Fig. 2). **a,b,** Mock population including only narrow (a) and wide (b) bursts. **c–f,** Mock populations following exponential declines plus bursts of different strengths: 0 (c), $1\sigma$ (d), $2\sigma$ (e) and $3\sigma$ (f). The legend in **b** applies to all panels. See Methods for a detailed explanation of the mock stellar populations test and the contents of this figure.

hypothesis of the link between these newly reported bursts of star formation in the Milky Way and interactions with Sgr.

However, at this point, we should note that Sgr is not the only Milky Way satellite that might be affecting our Galaxy at present. Among all of them, the Magellanic Clouds, owing to their mass and close pericentre, stand out. Several works suggest that they are in first infall towards the Milky Way$^{39}$. As a consequence, we cannot discount a possible effect of the Magellanic system on the Milky Way SFH in approximately the last 1–2 Gyr, and particularly around the time of the last enhancement ~0.1 Gyr ago, when the first pericentric approach of the Large Magellanic Cloud may have occurred$^{40}$. Given the observational nature of this work, we cannot rule out an effect from a combination of pericentre passages of less massive and more distant satellite galaxies on the observed SFH. However, evi-

dence piles up to suggest that Sgr is indeed the main actor behind these conspicuous enhancements in the Milky Way disk SFH.

It is worth noting that there might be a number of factors that could delay star formation linked to predicted pericentre passages. Among the possible reasons we can name the likely time lag between close passages and star formation enhancements or limitations on the ability of dynamical models to accurately predict orbital histories (for example, dynamical friction). Still, the episodic nature of the SFH recovered for the Milky Way disk clearly implies that drastic events should have occurred around 6, 2 and 1 Gyr ago. As detailed above, a wealth of evidence has accumulated to indicate that Sgr (first infall and subsequent pericentric passages) might be behind the reported star formation enhancements. These results perfectly supplement a recent characterization of the stellar content
<figure><img src="image_6.png" /></figure>

**Fig. 4** | SFH in the -2-kpc-radius bubble around the Sun distinguishing between the thin and thick disks. The SFR is shown as a function of look-back time for the thin and the thick disks (selected on the basis of tangential velocity criterion given in ref. ¹⁵). These SFHs have been obtained in a similar fashion to the orange line in Fig. 2. Green-shaded areas highlight the location of the detected star-forming bursts. The shaded envelopes represent errors in the SFH determinations. The inset shows the spatial distribution of the stars belonging to each substructure in the x-z plane. Some regions close to the Galactic plane have not been included in the analysis due to the high extinction in those directions.

of the Anticentre Stream, in which claims are made that this stream is the outcome of the interaction between the Milky Way and Sgr in its first infall¹¹. Of all of the pericentric passages of Sgr, the first infall is the one that seems to affect the SFH of the Milky Way disk more widely. Figure 4 clearly shows that, while the first Sgr passage triggered star formation in both disk kinematic components, subsequent passages were mainly able to produce substantial star formation in the thin disk.

The possibility that a satellite dwarf galaxy considerably less massive than its host is able to induce repeated massive events of star formation involving an important part of its disk is striking. Star formation induced by tidal perturbations has been reported before in galaxy pairs with lower mass ratios. An example in our vicinity is that of the M31–M33 pair (mass ratio ≈ 7:1). N-body simulations indicate that M31 and M33 had their latest pericentric passage ~2.6 Gyr ago⁴², and analysis of deep Hubble Space Telescope CMDs⁴³ reveals a clear episode of star formation peaking ~2 Gyr ago, in the outer regions of both M31 and M33, in close agreement with the timing of the suggested pericentre. In more distant galaxies, only ongoing star formation bursts associated with current interactions can be reasonably identified due to the difficulty of obtaining detailed time-resolved SFHs for these systems¹⁶.

Theoretical works have historically asserted that relatively massive mergers should enhance star formation in the embroiled systems²³,⁴⁴,⁴⁵. In particular, star forming bursts have been found within a cosmological framework in interacting systems with mass ratios of up to ~10:1 (ref. ¹⁷). Some other works claim that the star formation efficiency should decrease at high mass ratios, and indicate that little enhancement of star formation can be expected at the Sgr/Milky Way mass ratio⁴⁶. Other models of the effect of mergers on the SFH of interacting galaxies⁴⁴,⁴⁵ indicate that the induced star formation should be concentrated in the inner parts, as a consequence of gas

infall in response to tidal forces. This is expected even in the case of minor mergers (mass ratios of ~10:1 or lower). Some other simulations of major mergers (up to ratios of 2.5:1) even find that such star formation enhancements in central regions are accompanied by a moderate suppression of star formation at larger galactocentric distances⁴⁷, in clear contradiction to the global star formation episodes that we find.

However, modern attempts, including more complex and detailed modelling of star formation, are starting to cast doubts on these ideas. A proper understanding of the underlying physical processes behind star formation in galaxies and its inclusion in models is one of the major challenges for galaxy formation and evolution theory²³. The increasingly detailed modelling of the interstellar medium, including refinements such as multiphase structure, has provided results that have begun to change our understanding of star formation episodes in merging systems. Rather than inducing large-scale gas inflows towards the central parts, certain mergers can result in gas fragmentation into massive and dense clouds⁴⁸. Other aspects such as modifications of the modelling of star formation itself (based on density or shock-induced) and the effect of a hot gaseous halo (usually not included) also influence the way star formation is enhanced in simulations of interacting galaxies⁴⁹,⁵⁰. In all cases, implementation of more-realistic and higher-resolution modelling affects the starburst efficiency and, more importantly for this discussion, results in star formation enhancements in regions at larger galactocentric radii as opposed to the centrally concentrated star formation suggested by previous studies. Unfortunately, to our knowledge, there are no theoretical works where a single satellite could induce several global star formation enhancements as Sgr seems to be causing in the Milky Way disk. Therefore, although we can easily link the reported enhancements with possible peri-
centric passages of Sgr, we cannot pinpoint what exact physical mechanisms are triggering such events.

In this work we obtained a detailed SFH of the 2-kpc local volume. Strikingly, we find that such volume is characterized by an episodic SFH, with clear enhancements of star formation ~5.7, 1.9 and 1.0 Gyr ago. All evidence seems to suggest that recurrent interactions between the Milky Way and Sgr dwarf galaxy are behind such enhancements. These findings imply that low mass satellites, not only affect the Milky Way disk dynamics, but also are able to trigger notable events of star formation throughout its disk. The precise dating of such star forming episodes provided in this work sets useful boundary conditions to properly model the orbit of Sgr and its interaction with the Milky Way. In addition, this work provides important constraints on the modelling of the interstellar medium and star formation within hydrodynamical simulations, manifesting the need of understanding physical processes at subresolution scales and of further analysis to unveil the physical mechanisms behind global and repeated star formation events induced by satellite interaction.

## Methods

The European Space Agency (ESA) mission Gaia is providing an unprecedented view of our immediate cosmic neighbourhood through mapping of the stellar content of vast volumes around the Milky Way. In particular, the precise parallaxes to individual stars delivered in Gaia DR2¹ are among the most relevant breakthroughs. They allow deep CMDs in the absolute magnitude plane to be obtained for large volumes within the Milky Way for the first time¹⁵. The analysis of such CMDs enables the computation of detailed and quantitative SFHs⁵ extended from the very earliest times in the evolution of the Milky Way to the present. We should emphasize here that the quality and volume of the Gaia data in combination with advanced analysis techniques and a robust stellar evolution theoretical framework are crucial and essential ingredients in unveiling the SFH of the Milky Way with unprecedented detail and resolution.

**Gaia data and sample selection.** We make use of Gaia DR2 data¹ to infer the SFH characteristics of the disk of the Milky Way using the well-tested technique of CMD fitting. The original sample is composed of all Gaia DR2 stars within a sphere of ~2 kpc in radius (distance obtained by simply inverting the parallax, that is $\text{parallax} \ge 0.5$, not considering parallax offset, see first point in the list below) and brighter than $M_G=7$ (absolute magnitude). From the nearly 82 million stars fulfilling these initial criteria, we then restrict our sample to those displaying a relative uncertainty in the parallax below 20% (parallax\_over\_error > 5 (ref.⁵¹), keeping around 90% of the stars). The analysed sample of stars is drawn from this initial sample after some extra technical aspects are taken into account:

* The parallaxes provided by the Gaia Collaboration are known to be underestimated¹². Several works in the literature studying the Gaia DR2 parallax zero-point suggested different values ranging from 0.029 to 0.082 mas (refs. ⁵³⁻⁶⁰). We corrected the parallaxes of each and every star in our sample using a single value of 0.054 mas, not considering possible spatial variations, following ref. ⁵⁶. We have repeated the whole analysis using parallax offset values near the extremes of the proposed range, namely 0.029 and 0.082 mas. This allowed us to verify that the particular parallax offset adopted within reasonable values has little impact on the results presented in this work, although the percentage of old stars tends to be higher when using large values (that is, 0.082 mas).

* After considering this parallax offset, we have computed absolute magnitudes using the inverse of the parallax as a distance estimator to individual stars²². Different works in the literature provide Bayesian geometric distance estimates³⁰,⁶¹,⁶² (among other quantities) refining the Gaia DR2 data products that are pertinent for stellar population analysis. Although highly advisable for studies involving stars at large distances from the Sun, given the large parallaxes of our sources and their small parallax errors, the simple approximation followed in this work to compute distances is valid (see fig. 6 in ref. ⁶¹, fig. 13 in ref. ⁵⁶ or fig. 18 in ref. ⁶², to name but a few).

* We correct the Gaia magnitudes for interstellar extinction interpolating the colour excess from a three-dimensional map valid in the analysed volume⁶³. Different extinction coefficients are used¹⁵,⁶⁴ with little or no effect on the final results (see Fig. 2). SFH represented by solid lines are computed using the Casagrande coefficients⁴¹, whereas those depicted using dashed lines make use of the Babusiaux coefficients¹⁵. In the case of the Casagrande coefficients, temperatures for individual stars were estimated using a correlation between $G_{BP} - G_{RP}$ colour and effective temperature determined from a large sample of Gaia DR2 stars with measured temperatures (dependence with metallicity was not considered). Highly reddened stars (extinction in the Gaia G magnitude above 0.5, $A_G > 0.5$) are not included in the final sample (40% of

the stars are removed by applying this criterion). This cut introduces a spatial bias in the sample, as stars close to the plane and at larger distances from the Sun are affected by a larger amount of reddening. The inset of Fig. 4 shows the spatial location of the stars in the final sample. From the different SFH for the kinematic thin and thick disks shown in Fig. 4, it can be deduced that the spatial bias introduced by the reddening cut probably results in underestimated values of the young SFR, possibly also underestimating the heights of the youngest bursts.

* Given our knowledge of the stellar content of the Milky Way, we expect that its oldest main-sequence turn-off is located at around $M_G=4$. To ensure that we reach at least 1 magnitude deeper than this critical position, after the parallax offset is taken into account and the magnitudes are corrected by interstellar extinction, we keep for further analysis only those stars brighter than $M_G=5$. Almost half of the stars are removed in this step.

* Finally, following refs. ¹⁵,⁵², stars not fulfilling the following criterion were discarded (less than 1% of the stars): astrometric\_chi2\_al/(astrometric\_n\_good\_obs\_al-5) < 1.44×max(1, exp(-0.4×(phot\_g\_mean\_mag-19.5)))

Stars with no $G_{BP} - G_{RP}$ colour determination available are by definition not included in our sample (less than 0.2%).

At this point, we should highlight the fact that our sample is essentially complete, as all the analysed stars are brighter than $G=17.0$ (in particular, brighter than 16.5, apparent magnitude corresponding to a star of $M_G=5$ located at a maximum distance of 2 kpc, taking into account reddening). In the brightest end ($G<12$), Gaia DR2 is estimated to be 97% complete¹, so our sample is also basically complete at the bright end. This is important as any incompleteness or bias in the coverage of the observed CMD could affect the results and thus should be properly replicated in the synthetic CMD, which is the essential ingredient in the CMD fitting technique that we use in this work (see SFH from Gaia CMDs).

**SFH from Gaia CMDs.** The analysis of deep CMDs has long been used and demonstrated as an excellent way to study SFHs in dwarf galaxies within the Local Group²⁵,⁶⁶⁻⁷⁰. In this work, we derive SFH of the Milky Way disk by comparing Gaia DR2-observed CMDs (drawn from the original sample described previously) to synthetic ones following the methodology carefully outlined in ref. ⁶⁹. We used a synthetic CMD containing 200 million stars with ages and metallicityes ranging from 0.03 to 14 Gyr and 0.0001 to 0.05, respectively, calculated by adopting a constant SFR in the age-metallicity plane within the quoted ranges. This synthetic CMD was computed using the BaSTI stellar evolution library⁶⁹ in its solar-scaled version, with a Reimers mass loss parameter ($\eta$) of 0.2, assuming a Kroupa initial mass fraction⁷⁰, a binary fraction of 70% ($\beta$) and a minimum mass ratio for binaries of 0.1 ($q$). Regarding the treatment of binaries, we should clarify here that the stars in a binary system are drawn from the same (single-star) stellar evolutionary tracks, and their luminosities are combined. This should mimic the fact that Gaia DR2 contains some non-resolved binary systems.

The predictions from the stellar evolutionary models were translated to the Gaia photometric system by using bolometric correction tables following ref. ⁷¹ for Gaia DR2 (revised passbands corresponding to the data plotted in their fig. 21 (ascii data downloaded from the official Gaia home page) and the corresponding zero points in their table 1, column 4). Before comparing observed and synthetic CMDs, we need to simulate the effect of observational errors in the synthetic diagrams following the procedure explained in ref. ²². Given the small photometric errors of Gaia⁷¹, the major sources of uncertainty in the magnitude and colour of the stars stem from the error in the distance estimation used to transform apparent to absolute magnitudes using the parallax, and from the errors in the reddening correction (which dominate the error in the $G_{BP} - G_{RP}$ colour). We mimic observational errors in $M_G$ on a star-by-star basis using the equation $\Delta M_G = 2.17/\text{parallax\_over\_error}$ (see ref. ²² for more information). To this end, we recreate the observed parallax\_over\_error distribution in the synthetic CMD (including parallax offset). Then, each theoretical $M_G$ is modified by adding a correction following a Gaussian distribution centred at 0 with a standard deviation of $\Delta M_G$. Small errors in the Gaia photometry itself⁷¹, as well as in the reddening correction, are introduced in a similar fashion by adding a dispersion of 0.01 in the $G_{BP} - G_{RP}$ colour (see ref. ⁷¹ and especially ref. ²² for further details).

The SFH derivation was done using the code THESTORM²⁷,³⁷ (Tracing the evolution of the star formation rate and metallicity). In this step, the combination and relative weights of simple stellar populations (with a small range of age and metallicity) that best fit the observed CMD are determined. We note that both the star formation and the chemical enrichment rate are derived simultaneously and that no assumptions are adopted for any of them. These simple stellar populations are drawn from the synthetic CMD after simulation of observational errors. The comparison is carried out by counting the number of stars in colour-magnitude boxes of different sizes depending on their location within the CMD, where a number of bundles have been defined (see Fig. 1). Bundles are regions in the CMD sampled with an homogeneous box size, chosen according to the reliability of the stellar evolution models in that region, or the density of stars in them. We have shown that the exact bundle and box size configuration has little effect in the final SFH⁷⁴. Uncertainties on the SFH were computed as described in ref. ⁷⁵. They are assumed to be a combination in quadrature of the uncertainties due to the effect of binning in the colour-magnitude and age-metallicity planes (obtained from the
dispersion of different solutions obtained after shifting the bin limits), and those due to the effect of statistical sampling in the observed CMD (calculated from the dispersion of several solutions derived after resampling the observed CMD following Poissonian statistics). The applied methodology in its current version does not take into account the effect on the error computation of secondary factors such as a possible degeneracy in the solution due to secondary solutions (close to the best fit), or the choice of the stellar population parameters adopted to calculate the synthetic CMD (initial mass fraction, binary simulation and so on), which do not exactly match those of the observed population. For a more detailed explanation on the fitting procedure and error calculation, we refer the reader to refs. 68,73,75.

For a proper computation of the SFH coming from a given sample of observed stars, the synthetic CMD to which we compare should greatly outnumber the observed stars. For this reason, and to enable the computation of the SFH of such a vast volume (~24 million stars in the ~2-kpc-radius sphere, after quality cuts have been applied), we randomly select several subsamples from the original one, containing up to ~250,000 stars. Figure 1 shows one of these fits (Fig. 1b) compared with the observed CMD (Fig. 1a). Figure 1d shows the residual difference between the observed and the best-fit model CMD in Poissonian σ. Note that in the great majority of the CMD the residuals are lower than 2σ. The largest difference between the model and the observed CMD are observed in the red clump and in the fainter part of the main sequence. However, these differences are usually within 2–3σ and have understood origins. (1) The presence of higher than average residuals in the region of the CMD corresponding to the core He-burning stage, is mainly due to the fact that the adopted version of the BaSTI library relies on the use of conductive opacity prescriptions that have been improved in the new BaSTI library<sup>76,77</sup> and that have the effect of slightly decreasing the brightness of the horizontal branch stellar models, which would contribute to improving the match between the synthetic and the observed CMD; additionally, the exact morphology—in particular the colour distribution—of the red clump strongly depends on the adopted mass loss efficiency during the previous RGB stage, for which we have adopted a (fixed) average efficiency. (2) In the case of the main sequence below the old main-sequence turn-off, the positive residuals are related to the stars for which the reddening correction was underestimated, and that, therefore, could not be reproduced in the model. These stars are shifted to a redder position, and thus are absent from their true, blue position in the main sequence (thus the negative residual there). Figure 2 shows the recovered SFH for all the trials we have performed. Different lines represent the SFHs from the analysis of different random selections and different extinction coefficients to ensure the robustness of the result and the negligible effect of reducing the sample under analysis. The final solution (orange line in Fig. 2) is computed after averaging all the solutions. Vertical error bars encompassed the 16th and 84th percentiles of all solutions.

The ability of this method to precisely recover the age and duration of well-defined and isolated bursts of star formation has been extensively tested<sup>75,78</sup> using mock stellar populations. These works (and references therein) show that (1) a high precision (ability to pinpoint the location of the peak of the burst) is achieved in the recovery of the age of the star formation events (typically within few hundred million years for older events, and even under a hundred million years for younger events) and that (2) the age resolution (how narrow a burst should be to be detected) gets worse for older stellar populations, and may range from a few hundred million years for young and intermediate ages, to 1–3 Gyr for old stars. However, our experience recommends performing new tests with mock populations for each dataset, carefully mimicking all experimental conditions. Therefore, we have performed similar tests tailored for the particular case of Gaia DR2 data, and focusing on the results presented here. We created a series of mock stellar populations: a set of ~200,000 synthetic stars from the BaSTI models following a particular SFH (exponential decline, burst of star formation and combinations of exponential decline and bursts of different widths and strengths and so on). We then simulated the observational errors to such synthetic mock CMDs as explained above to finally obtain a recovered SFH using THESTORM (same methodology as that described above). This experiment is shown in Fig. 3 and allows us to reach conclusions about the significance of the star formation bursts and their maximum duration. They also provide information about the precision of the age determination and the age resolution at different ages of the derived SFH.

First, Fig. 3a,b,e,f confirms that a very high precision in the age recovery is achieved for the current dataset, both when the simulated star formation bursts are discrete and well separated (Fig. 3a,b), or in presence of a background star formation (Fig. 3e,f). Figure 3a shows the case of very narrow bursts, whereas Fig. 3b shows the case of wide bursts. In the latter, the star formation peak is recovered very approximately at the mean age of the wide burst. Second, Fig. 3a illustrates the age resolution of the SFH presented in this paper, and how it changes with age: it shows that we can recover very narrow bursts with our methodology and present dataset. The mock bursts shown have a width of 0.1 Gyr, 0.2 Gyr and 0.5 Gyr at the ages of the bursts of star formation measured in the Gaia data (from young to old, respectively). As can be seen, we clearly reproduce their presence and input width (a small tail towards older ages can be seen in the burst at 5.5 Gyr ago). Horizontal error bars in Fig. 2 indicate our age resolution in the three main star forming bursts detected. Third, Fig. 3b allows us to estimate the burst duration. The input mock CMD had bursts ranging from 0.7–1.1, 1.6–2.2 and 4.5–6.5 Gyr ago, respectively, and the recovered SFH closely mimics the width of the bursts that characterize the SFH derived from the Gaia data. Therefore, we can conclude that the durations of the different bursts are compatible with being (at most) 0.4, 0.6 and 2.0 Gyr, respectively (from young to old). Finally, Fig. 3c–f (exponential plus bursts of different strengths) allows us to conclude that the bursts that we recover for the ~2-kpc-radius around the solar neighbourhood (see Fig. 2) are significant above 3–4σ with respect to the background formation.

## Data availability

All data analysed in this paper are publicly available from the Gaia DR2 archive (http://gea.esac.esa.int/archive/). The datasets containing information not available in public catalogues that is necessary to reproduce this work and the corresponding figures are available as Supplementary Data 1. An explanatory README file is included. Supplementary Data 1 consists of the following files: (1) For Figs. 1 and 2, two examples of randomly selected samples of ~250,000 stars (with different extinction coefficients), together with the solution from the THESTORM code. This includes two tables directly retrieved from the Gaia archive as described in the Methods, supplemented by extinction information on a star-by-star basis. (2) For Fig. 3, six mock population CMDs together with the outputs from THESTORM. (3) Examples of randomly selected samples of 250,000 stars for thin and thick disk stars, together with the solutions from the THESTORM code (for Fig. 4). Other data not included in the above-mentioned link are available from the corresponding author on reasonable request.

## Code availability

The code used to interpolate the three-dimensional extinction maps<sup>63</sup> can be retrieved from https://github.com/edober/dustMAPS_3d.

Received: 20 September 2019; Accepted: 26 March 2020; Published online: 25 May 2020

## References

1. Gaia Collaboration et al. Gaia Data Release 2: summary of the contents and survey properties. *Astron. Astrophys.* **616**, A1 (2018).

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. Mor, R., Robin, A. C., Figueras, F., Roca-Fàbrega, S. & Luri, X. Gaia DR2 reveals a star formation burst in the disc 2-3 Gyr ago. *Astron. Astrophys.* **624**, L1 (2019).

4. Heavens, A., Panter, B., Jimenez, R. & Dunlop, J. The star-formation history of the Universe from the stellar populations of nearby galaxies. *Nature* **428**, 625–627 (2004).

5. Perryman, M. A. C. et al. The HIPPARCOS Catalogue. *Astron. Astrophys.* **323**, L49–L52 (1997).

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

7. Vergely, J.-L., Köppen, J., Egret, D. & Bienaymé, O. An inverse method to interpret colour-magnitude diagrams. *Astron. Astrophys.* **390**, 917–929 (2002).

8. 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).

9. Haywood, M. et al. Phylogeny of the Milky Way's inner disk and bulge populations: implications for gas accretion, (the lack of) inside-out thick disk formation, and quenching. *Astron. Astrophys.* **618**, A78 (2018).

10. Bernard, E. J. in IAU Symposium 334: *Rediscovering Our Galaxy* (eds Chiappini, C. et al.) 158–161 (IAU, 2018).

11. Feltzing, S., Holmberg, J. & Hurley, J. R. The solar neighbourhood age-metallicity relation—does it exist? *Astron. Astrophys.* **377**, 911–924 (2001).

12. Bergemann, M. et al. The Gaia-ESO Survey: radial metallicity gradients and age-metallicity relation of stars in the Milky Way disk. *Astron. Astrophys.* **565**, A89 (2014).

13. 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).

14. Gaia Collaboration et al. Gaia Data Release 2: mapping the Milky Way disc kinematics. *Astron. Astrophys.* **616**, A11 (2018).

15. Gaia Collaboration et al. Gaia Data Release 2: observational Hertzsprung-Russell diagrams. *Astron. Astrophys.* **616**, A10 (2018).

16. Kennicutt, R. C. Jr., Keel, W. C., van der Hulst, J. M., Hummel, E. & Roettiger, K. A. The effects of interactions on spiral galaxies. II: disk star-formation rates. *Astron. J.* **93**, 1011–1023 (1987).

17. Tissera, P. B., Domíguez-Tenreiro, R., Scannapieco, C. & Sáiz, A. Double starbursts triggered by mergers in hierarchical clustering scenarios. *Mon. Not. R. Astron. Soc.* **333**, 327–338 (2002).

18. Ellison, S. L., Mendel, J. T., Patton, D. R. & Scudder, J. M. Galaxy pairs in the Sloan Digital Sky Survey–VIII. The observational properties of post-merger galaxies. *Mon. Not. R. Astron. Soc.* **435**, 3627–3638 (2013).
19. Tinsley, B. M. & Larson, R. B. Stellar population explosions in proto-elliptical galaxies. *Mon. Not. R. Astron. Soc.* **186**, 503–517 (1979).

20. Freeman, K. & Bland-Hawthorn, J. The new galaxy: signatures of its formation. *Annu. Rev. Astron. Astrophys.* **40**, 487–537 (2002).

21. Meza, A., Navarro, J. F., Abadi, M. G. & Steinmetz, M. Accretion relics in the solar neighbourhood: debris from *ω* Cen's parent galaxy. *Mon. Not. R. Astron. Soc.* **359**, 93–103 (2005).

22. Gallart, C. et al. Uncovering the birth of the Milky Way through accurate stellar ages with Gaia. *Nat. Astr. J.* **3**, 932–939 (2019).

23. Naab, T. & Ostriker, J. P. Theoretical challenges in galaxy formation. *Annu. Rev. Astron. Astrophys.* **55**, 59–109 (2017).

24. Ibata, R. A., Gilmore, G. & Irwin, M. J. A dwarf satellite galaxy in Sagittarius. *Nature* **370**, 194–196 (1994).

25. Belokurov, V. et al. The field of streams: Sagittarius and its siblings. *Astrophys. J.* **642**, L137–L140 (2006).

26. Law, D. R. & Majewski, S. R. The Sagittarius dwarf galaxy: a model for evolution in a triaxial Milky Way halo. *Astrophys. J.* **714**, 229–254 (2010).

27. Purcell, C. W., Bullock, J. S., Tollerud, E. J., Rocha, M. & Chakrabarti, S. The Sagittarius impact as an architect of spirality and outer rings in the Milky Way. *Nature* **477**, 301–303 (2011).

28. Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N. & Besla, G. The influence of Sagittarius and the Large Magellanic Cloud on the stellar disc of the Milky Way Galaxy. *Mon. Not. R. Astron. Soc.* **481**, 286–306 (2018).

29. Eggen, O. J. The motions of the A Stars at the North Galactic Pole. *Publ. Astron. Soc. Pac.* **81**, 741 (1969).

30. Siebert, A. et al. Detection of a radial velocity gradient in the extended local disc with RAVE. *Mon. Not. R. Astron. Soc.* **412**, 2026–2032 (2011).

31. Williams, M. E. K. et al. The wobbly galaxy: kinematics north and south with RAVE red-clump giants. *Mon. Not. R. Astron. Soc.* **436**, 101–121 (2013).

32. Gómez, F. A. et al. Vertical density waves in the Milky Way disc induced by the Sagittarius dwarf galaxy. *Mon. Not. R. Astron. Soc.* **429**, 159–164 (2013).

33. Antoja, T. et al. A dynamically young and perturbed Milky Way disk. *Nature* **561**, 360–362 (2018).

34. Laporte, C. F. P., Minchev, I., Johnston, K. V. & Gómez, F. A. Footprints of the Sagittarius dwarf galaxy in the Gaia data set. *Mon. Not. R. Astron. Soc.* **485**, 3134–3152 (2019).

35. de la Vega, A., Quillen, A. C., Carlin, J. L., Chakrabarti, S. & D'Onghia, E. Phase wrapping of epicyclic perturbations in the Wobbly Galaxy. *Mon. Not. R. Astron. Soc.* **454**, 933–945 (2015).

36. 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).

37. de Boer, T. J. L., Belokurov, V. & Koposov, S. The star formation history of the Sagittarius stream. *Mon. Not. R. Astron. Soc.* **451**, 3489–3503 (2015).

38. Siegel, M. H. et al. The ACS survey of galactic globular clusters: M54 and Young populations in the sagittarius dwarf spheroidal galaxy. *Astrophys. J.* **667**, L57–L60 (2007).

39. Besla, G. et al. Are the Magellanic Clouds on their first passage about the Milky Way? *Astrophys. J.* **668**, 949–967 (2007).

40. Laporte, C. F. P., Gómez, F. A., Besla, G., Johnston, K. V. & Garavito-Camargo, N. Response of the Milky Way's disc to the Large Magellanic Cloud in a first infall scenario. *Mon. Not. R. Astron. Soc.* **473**, 1218–1230 (2018).

41. Laporte, C. F. P., Belokurov, V., Koposov, S. E., Smith, M. C. & Hill, V. Chemo-dynamical properties of the Anticenter Stream: a surviving disc fossil from a past satellite interaction. *Mon. Not. R. Astron. Soc.* **492**, L61–L65 (2020).

42. McConnachie, A. W. et al. The remnants of galaxy formation from a panoramic survey of the region around M31. *Nature* **461**, 66–69 (2009).

43. Bernard, E. J. et al. The star formation history and dust content in the far outer disc of M31. *Mon. Not. R. Astron. Soc.* **420**, 2625–2643 (2012).

44. Mihos, J. C. & Hernquist, L. Triggering of starbursts in galaxies by minor mergers. *Astrophys. J.* **425**, L13–L16 (1994).

45. Hernquist, L. & Mihos, J. C. Excitation of activity in galaxies by minor mergers. *Astrophys. J.* **448**, 41 (1995).

46. Cox, T. J., Jonsson, P., Somerville, R. S., Primack, J. R. & Dekel, A. The effect of galaxy mass ratio on merger-driven starbursts. *Mon. Not. R. Astron. Soc.* **384**, 386–409 (2008).

47. Moreno, J. et al. Mapping galaxy encounters in numerical simulations: the spatial extent of induced star formation. *Mon. Not. R. Astron. Soc.* **448**, 1107–1117 (2015).

48. Teyssier, R., Chapon, D. & Bournaud, F. The driving mechanism of starbursts in galaxy mergers. *Astrophys. J.* **720**, L149–L154 (2010).

49. Chien, L.-H. & Barnes, J. E. Dynamically driven star formation in models of NGC 7252. *Mon. Not. R. Astron. Soc.* **407**, 43–54 (2010).

50. Moster, B. P., Macciò, A. V., Somerville, R. S., Naab, T. & Cox, T. J. The effects of a hot gaseous halo in galaxy major mergers. *Mon. Not. R. Astron. Soc.* **415**, 3750–3770 (2011).

51. Luri, X. et al. Gaia Data Release 2: using Gaia parallaxes. *Astron. Astrophys.* **616**, A9 (2018).

52. Lindegren, L. et al. Gaia Data Release 2: the astrometric solution. *Astron. Astrophys.* **616**, A2 (2018).

53. Stassun, K. G. & Torres, G. Evidence for a systematic offset of -80 µas in the Gaia DR2 Parallaxes. *Astrophys. J.* **862**, 61 (2018).

54. Riess, A. G. et al. Milky Way Cepheid standards for measuring cosmic distances and application to Gaia DR2: implications for the Hubble constant. *Astrophys. J.* **861**, 126 (2018).

55. 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. *Astrophys. J.* **878**, 136 (2019).

56. Schönrich, R., McMillan, P. & Eyer, L. Distances and parallax bias in Gaia DR2. *Mon. Not. R. Astron. Soc.* **487**, 3568–3580 (2019).

57. Khan, S. et al. New light on the Gaia DR2 parallax zero-point: influence of the asteroseismic approach, in and beyond the Kepler field. *Astron. Astrophys.* **628**, A35 (2019).

58. 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).

59. 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).

60. Leung, H. W. & Bovy, J. Simultaneous calibration of spectro-photometric distances and the Gaia DR2 parallax zero-point offset with deep learning. *Mon. Not. R. Astron. Soc.* **489**, 2079–2096 (2019).

61. 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).

62. Anders, F. et al. Photo-astrometric distances, extinctions, and astrophysical parameters for Gaia DR2 stars brighter than $G = 18$. *Astron. Astrophys.* **628**, A94 (2019).

63. 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).

64. 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).

65. 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).

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

67. 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).

68. 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).

69. 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).

70. Kroupa, P. On the variation of the initial mass function. *Mon. Not. R. Astron. Soc.* **322**, 231–246 (2001).

71. Evans, D. W. et al. Gaia Data Release 2: photometric content and validation. *Astron. Astrophys.* **616**, A4 (2018).

72. 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).

73. 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).

74. 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).

75. 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).

76. Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M. & Salaris, M. Updated electron-conduction opacities: the impact on low-mass stellar models. *Astrophys. J.* **661**, 1094–1104 (2007).

77. Hidalgo, S. L. et al. The updated BaSTI stellar evolution models and isochrones. I. Solar-scaled calculations. *Astrophys. J.* **856**, 125 (2018).

78. Rusakov, V. et al. The bursty star formation history of the Fornax dwarf spheroidal galaxy revealed with the HST. Preprint at https://arxiv.org/abs/0202.09714 (2020).

79. Hunter, J. D. Matplotlib: a 2d graphics environment. *Comput. Sci. Eng.* **9**, 90–95 (2007).
80. Astropy Collaboration et al. Astropy: a community Python package for astronomy. Astron. Astrophys. 558, A33 (2013).

81. Astropy Collaboration et al. The Astropy Project: building an open-science project and status of the v2.0 core package. Astron. J. 156, 123 (2018).

## Acknowledgements

We thank K. C. Freeman, C. Bailer-Jones, G. Battaglia, M. Beasley, C. Brook, C. Dalla Vecchia, J. Falcón-Barroso, R. Leaman and I. Pérez for useful discussions. T.R.-L. and C.G. acknowledge financial support through grant numbers (AEI/FEDER, UE) AYA2017-89076-P, AYA2016-77237-C3-1-P (RAVET project) and AYA2015-63810-P, as well as from the Ministerio de Ciencia, Innovación y Universidades (MCIU) through the State Budget and the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community through the Regional Budget (including IAC project, TRACES). T.R.-L. is supported by a MCIU Juan de la Cierva-Formación grant (FJCI-2016-30342). S.C. acknowledges support from Premiale INAF "MITIC" and grant number AYA2013-42781P from the Ministry of Economy and Competitiveness of Spain; he has also been supported by the INFN (Iniziativa specifica TAsP). 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. This research makes use of Python (version 3.6.7, http://www.python.org); Matplotlib (version 3.0.0)⁷⁹, a suite of open-source python modules that provide a framework for creating scientific plots; and Astropy (version 3.0.5), a community-developed core Python package for astronomy⁸⁰,⁸¹.

## Author contributions

The manuscript was written by T.R.-L. and C.G. T.R.-L. and C.G. defined the final samples under analysis and extracted the SFH presented in this work. The software to analyse Gaia DR2 data was written by T.R.-L. and E.J.B. S.C. contributed to the tools used to generate the synthetic CMDs and evolutionary model predictions in the Gaia photometric system. All authors contributed to the interpretation and analysis of the results.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41550-020-1097-0.

**Correspondence and requests for materials** should be addressed to T.R.-L.

**Peer review information** *Nature Astronomy* thanks Michele Cignoni and Roger Mor for their contribution to the peer review of this work.

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

**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 2020