We illustrate the formation and evolution of the Milky Way over cosmic time utilizing a sample of 10 million red giant stars with full chemodynamical information including metallicities and abundances from lowresolution Gaia XP spectra. The evolution of angular momentum as a function of metallicitya rough proxy for stellar age particularly for highFe starsdisplays three distinct phases the disordered and chaotic protogalaxy the kinematically hot old disk and the kinematically cold young disk. The old high disk starts at FeH  1.0 spinning up from the nascent protogalaxy and then exhibiting a smooth cooldown toward more ordered and circular orbits at higher metallicities. The young low disk is kinematically cold throughout its metallicity range with its observed properties modulated by a strong radial gradient. We interpret these trends using Milky Way analogs from the TNG50 cosmological simulation identifying one that closely matches the kinematic evolution of our galaxy. This halos protogalaxy spins up into a relatively thin and misaligned high disk at early times which is subsequently heated and torqued by a major gasrich merger. The merger contributes a large amount of lowmetallicity gas and angular momentum from which the kinematically cold low stellar disk is subsequently born. This simulated history parallels several observed features of the Milky Way particularly the decisive GaiaSausageEnceladus merger that likely occurred at z  2. Our results provide an allsky perspective on the emerging picture of our galaxys threephase formation impelled by the three physical mechanisms of spinup merger and cooldown.
['Milky Way Galaxy', 'Milky Way disk', 'Milky Way dynamics', 'Milky Way formation', 'Milky Way evolution', '1054', '1050', '1051', '1053', '1052', 'Astrophysics - Astrophysics of Galaxies']
The Threephase Evolution of the Milky Way
{'The Three-Phase Evolution of the Milky Way': "Vedant Chandra, 1, 2 Vadim A. Semenov, 1 Hans-Walter Rix, 2 Charlie Conroy, 1 Ana Bonaca, 3 Rohan P. Naidu, 4, ∗ Ren'e Andrae, 2 Jiadong Li ( 李 佳 东 ) , 2 and Lars Hernquist 1  \n1 Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA 2 Max-Planck-Institut fur Astronomie, Konigstuhl 17, D-69117 Heidelberg, Germany 3 The Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA 4 MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Ave., Cambridge, MA 02139, USA", 'ABSTRACT': "We illustrate the formation and evolution of the Milky Way over cosmic time, utilizing a sample of 10 million red giant stars with full chemodynamical information, including metallicities and α -abundances from low-resolution Gaia XP spectra. The evolution of angular momentum as a function of metallicity - a rough proxy for stellar age, particularly for high-[ α /Fe] stars - displays three distinct phases: the disordered and chaotic protogalaxy, the kinematically-hot old disk, and the kinematicallycold young disk. The old highα disk starts at [Fe/H] ≈ -1 . 0, 'spinning up' from the nascent protogalaxy, and then exhibits a smooth 'cooldown' toward more ordered and circular orbits at higher metallicities. The young lowα disk is kinematically cold throughout its metallicity range, with its observed properties modulated by a strong radial gradient. We interpret these trends using Milky Way analogs from the TNG50 cosmological simulation, identifying one that closely matches the kinematic evolution of our Galaxy. This halo's protogalaxy spins up into a relatively thin and misaligned highα disk at early times, which is subsequently heated and torqued by a major gas-rich merger. The merger contributes a large amount of low-metallicity gas and angular momentum, from which the kinematically cold lowα stellar disk is subsequently born. This simulated history parallels several observed features of the Milky Way, particularly the decisive 'GSE' merger that likely occurred at z ≈ 2. Our results provide an all-sky perspective on the emerging picture of our Galaxy's three-phase formation, impelled by the three physical mechanisms of spinup, merger, and cooldown.  \nKeywords: Milky Way Galaxy (1054), Milky Way disk (1050), Milky Way dynamics (1051), Milky Way formation (1053), Milky Way evolution (1052)", '1. INTRODUCTION': "It is a central goal of modern astronomy to piece together a coherent formation history of the Milky Way (MW). More broadly, we wish to understand how much formation memory a galaxy retains over cosmic timescales, and consequently investigate the physical mechanisms that shape disk galaxies. For both aspects, the Milky Way is the ideal 'near-field' laboratory with its detailed star-by-star information. Finally, the Milky Way is set apart from all other galaxies in the Universe  \nCorresponding author: Vedant Chandra  \nvedant.chandra@cfa.harvard.edu  \n∗ NASA Hubble Fellow  \nbecause it hosts us, and it is inherently desirable to understand how our own cosmic home came to be.  \nFrom a theoretical perspective, disk galaxies have long been understood to emerge as a consequence of dissipative gas collapse in dark matter halos, and the retention of angular momentum acquired during large-scale structure formation (Eggen et al. 1962; Mestel 1963; Peebles 1969; White & Rees 1978; Fall & Efstathiou 1980; Mo et al. 1998). These ideas formed the basis of subsequent analytic and semi-analytic models which succeeded in producing reasonably realistic disk galaxies (e.g., Silk 1997; Dalcanton et al. 1997; Dekel et al. 2009b, 2013; Somerville & Dav'e 2015; Stevens et al. 2016). Similar progress has been made in numerical simulations, with a variety of works exploring disk formation in cosmological simulations, from idealized and zoom-in scales to large  \nvolume scales (e.g. Bournaud et al. 2009; Hopkins et al. 2009, 2023; Minchev et al. 2012; Bird et al. 2013; Agertz & Kravtsov 2015; Genel et al. 2015; Schaye et al. 2015; Agertz & Kravtsov 2016; Grand et al. 2017; Naab & Ostriker 2017; Garrison-Kimmel et al. 2018; Jiang et al. 2019; Tacchella et al. 2019; Pillepich et al. 2019; Stern et al. 2021, 2023; Yu et al. 2021, 2023; Agertz et al. 2021; Dillamore et al. 2022; Hafen et al. 2022; Khoperskov et al. 2022a,b; Gurvich et al. 2023; Wang & Lilly 2023; Pillepich et al. 2023; Semenov et al. 2023a,b; Santucci et al. 2023). These simulations showed the crucial importance of feedback-driven galactic outflows-which retain the angular momentum that galaxies inherit from their dark matter halos-as well as the importance of stochastic and occasionally violent processes like mergers and variable gas accretion flows.  \nThe physical mechanism of disk formation is still a subject of active debate. For example, Stern et al. (2021, 2023, see also Hafen et al. 2022) put forward a model where cold disk formation is triggered by the assembly of a hot inner circumgalactic halo, and an associated reduction in the stochasticity of star formation. On the other hand, Hopkins et al. (2023) use controlled numerical experiments based on FIRE-2 cosmological simulations to argue that the fundamental causative variable governing disk formation is a centrally concentrated density profile, which suppresses destructive radial gas motions and permits the coherent growth of angular momentum.  \nRecently, Semenov et al. (2023a, hereafter S23) exploited the diversity of MW-like galaxies in the Illustris TNG50 cosmological simulations (Pillepich et al. 2018; Nelson et al. 2019) to interpret observational evidence that the MW disk formed atypically early (Belokurov & Kravtsov 2022, 2023; Conroy et al. 2022). They corroborate the idea that disk formation is most closely linked with mass assembly, suggesting that the MW may have experienced unusually rapid mass growth at early times. At the same time, the underlying physical driver of disk formation is hard to pinpoint as both hot halo formation and gravitational potential steepening tend to occur concurrently with the disk spinup (Semenov et al. 2023b), and the potential steepening continues throughout the subsequent evolution of the disk. This suggests that a centrally concentrated potential might be a consequence of disk formation rather than its cause (Dillamore et al. 2023).  \nIn the past decade, our perspective of the Milky Way has expanded greatly through the proliferation of photometric, chemical, and kinematic data for millions of stars (see the reviews of Ivezi'c et al. 2012; Rix & Bovy 2013; Bland-Hawthorn & Gerhard 2016). It is reasonably certain that we reside in a barred spiral galaxy, with  \ntwo key characteristics of the MW being a total stellar mass ≈ 5 × 10 10 M ⊙ and a relatively quiescent presentday star formation rate of ≈ 1 . 6 M ⊙ yr -1 (Licquia & Newman 2015; Bland-Hawthorn & Gerhard 2016). The spatial, kinematic, chemical, and temporal structure of the disk has been exquisitely mapped since the advent of all-sky spectroscopic surveys, augmented in the past half-decade with astrometry from the Gaia space observatory (e.g., Bovy et al. 2012b,a; Haywood et al. 2013; Bensby et al. 2014; Hayden et al. 2015; Bovy et al. 2016; Weinberg et al. 2019, 2022; Xiang & Rix 2022; Gaia Collaboration et al. 2023a,b).  \nApart from a potentially early-forming disk, there are several other features that set the MW apart from typical disk galaxies. From an integrated perspective, the MWlies in the intermediate 'green valley' between most star forming and quiescent galaxies (Mutch et al. 2011; Licquia & Newman 2015). Furthermore, it is exceptionally disk-dominated compared to other known galaxies with comparable stellar masses (Kent et al. 1991; Kauffmann et al. 2003). Finally, the MW is not isolated, and belongs to a 'Local Group' of galaxies with a nearby pair of star-forming satellites, the Magellanic Clouds.  \nThe MW disk is known to have a striking bimodality in the abundance of α elements at constant metallicity (e.g., Fuhrmann 1998; Bensby et al. 2014; Hayden et al. 2015). One way to explain this is through the infall of fresh gas from a major merger, which would dilute the overall metallicity of the star-forming gas without significantly changing the abundance ratios (Chiappini et al. 1997; Grisoni et al. 2017; Buck et al. 2020; Buck 2020; Renaud et al. 2021b,a). However, it is also possible to produce this chemical bimodality solely through secular processes of chemical enrichment and radial migration (Schonrich & Binney 2009a,b; Minchev et al. 2013; Loebman et al. 2016; Sharma et al. 2021; Johnson et al. 2021). There is tentative evidence from JWST that such a chemical bimodality is absent from our neighbor, Andromeda (Nidever et al. 2023).  \nRecent work has advanced our understanding of various components of the Galaxy, including the ancient protogalaxy (e.g., Rix et al. 2022; Horta et al. 2023), the birth of the disks (e.g., Xiang & Rix 2022), the in-situ halo (e.g., Bonaca et al. 2020; Belokurov et al. 2020), and the accreted halo (e.g., Ibata et al. 1994; Helmi & White 1999; Majewski et al. 2003; Newberg et al. 2009; Nissen & Schuster 2010; Belokurov et al. 2018; Helmi et al. 2018; Naidu et al. 2020; Bonaca et al. 2021; Malhan et al. 2022). Most recently, Belokurov & Kravtsov (2022) isolate in-situ stars in the APOGEE survey and find a dramatic increase in the rotational support of stars at [Fe/H] ≈ -1 . 0, marking the transition between  \nthe protogalaxy ('Aurora' in their terminology) and the ancient disk (see also Rix et al. 2022). Conroy et al. (2022) study this transition with abundances and ages from the H3 Survey, finding that it coincides with a nonmonotonic rise in [ α /Fe] abundances, suggesting a nearinstantaneous transition in the star formation efficiency or gas inflow (see Chen et al. 2023). Here, we seek to provide a global perspective on these results using the latest data from the Gaia space observatory.  \nThe paper is structured as follows. We summarize the construction of our data set - 10 million stars with full chemodynamical information from Gaia DR3 - in § 2. In § 3, we present a unified perspective on the formation of the Milky Way, focusing on the observed three-phase evolution of angular momentum over cosmic time. Our goal is to present a compact and cohesive set of figures that qualitatively summarize the evolution of the Galaxy, and provide clues about its history. Towards this aim, in § 4 we interpret our data using simulations of MW analogs from the IllustrisTNG50 cosmological simulations, focusing on one particularly similar analog as an illustrative example. We discuss and contextualize the three-phase birth of the MW in § 5, and summarize our conclusions in § 6.", '2. GAIA DR3 DATA': "The third data release of the Gaia space observatory provided low-resolution 'XP' prism spectroscopy for ≈ 200 million stars brighter than G = 17 . 65, and radial velocities for a subset of ≈ 30 million stars (Gaia Collaboration et al. 2022; De Angeli et al. 2022; Montegriffo et al. 2022; Katz et al. 2022). Subsequent work has augmented this dataset by producing precise estimates of the overall metallicity [Fe/H] and alpha-abundance [ α /Fe] using XP spectra (Andrae et al. 2023; Zhang et al. 2023; Leung & Bovy 2023; Li et al. 2023). This sample is powerful not only because of its unprecedented size and homogeneity, but also because of its all-sky coverage and simple selection function, unlike previous pencilbeam spectroscopic surveys. This is therefore the ideal dataset to qualitatively explore a unified picture of our Galaxy's formation.  \nOur final observational sample contains 9,907,115 RGB stars with metallicities, [ α /Fe] abundances, and full 6D phase space coordinates. Here, we outline the construction of this dataset.", '2.1. Metallicities and [ α /Fe] with Gaia XP': "Andrae et al. (2023) construct a data-driven 'XGBoost' model to estimate stellar metallicity [M/H] from low-resolution XP spectra, utilizing a training set from the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017), augmented  \n<!-- image -->  \n<!-- image -->  \nFigure 1. Top and middle: Logarithmic Galactic map of our RGB sample from Gaia DR3, for two Gaia XP metallicity selections. Several selection effect artifacts are visible, primarily due to Gaia 's scanning law; however, these spatial variations are below the ≈ 10% level. Bottom: Galactocentric distribution of our sample, overlaid on an artist's rendering of the Milky Way disk (NASA/JPL-Caltech/R. Hurt).  \n<!-- image -->  \nwith metal-poor stars from Li et al. (2022). These metallicities have been validated to be precise and accurate compared to external spectroscopic surveys. Although the model formally predicts bulk metallicity [M/H], the APOGEE values in the training set are calibrated to match the iron abundance [Fe/H]. We therefore utilize the notation [Fe/H] for metallicity throughout this work.  \nFor more details on feature selection and catalog validation, we refer to Andrae et al. (2023).  \nFrom the metallicity sample of Andrae et al. (2023), we construct a clean sample of RGB stars with Gaia DR3 radial velocities. Specifically, we adopt the following selection criteria:  \n- · phot bp mean mag < 16 . 5\n- · ϖ/σ ϖ > 5\n- · teff xgboost < 5200\n- · logg xgboost < 3 . 5\n- · M W 1 > -0 . 3 -0 . 006 · (5500 -teff xgboost )\n- · M W 1 > -0 . 01 · (5300 -teff xgboost )\n- · ( G -W 2 ) < 0 . 2 + 0 . 77( G BP -W 1 )  \nwhere teff xgboost and logg xgboost are from Andrae et al. (2023), ϖ is the Gaia parallax, and M W 1 = W 1 +5 · log 10 ( ϖ/ 100) is the absolute magnitude in the WISE W1 passband. The first two cuts select bright stars with reliable parallaxes and enough flux in the bluer BP spectrum to derive robust metallicities. The temperature selection removes hot stars that have spuriously low measured metallicities, and the final four cuts isolate red giants on the Hertzsprung-Russel diagram. We excise 6841 stars that have a > 50% probability of residing in a globular cluster according to Vasiliev & Baumgardt (2021). This leaves us with ≈ 11 million RGB stars with Gaia metallicities and radial velocities.  \nLi et al. (2023) present a catalog of [ α /Fe] measurements inferred from Gaia XP spectra through a neural network model. This catalog has been cross-validated to match existing spectroscopic surveys, with median RMS errors of ≈ 0 . 05 dex in [ α /Fe] for stars in our sample. We cross-match our clean RGB sample to this catalog, obtaining our final sample of ≈ 10 million stars with [ α /Fe] information. The spatial distribution of these stars is illustrated in Figure 1. Simply selecting metal-poor stars on this all-sky map reveals the protogalactic population of the Milky Way (Rix et al. 2022), whose evolution is further investigated in our work. Our sample uniformly covers a significant fraction of the Milky Way's stellar disk (bottom panel of Figure 1).  \nFigure 2 summarizes our sample in Galactocentric radius R and midplane height z coordinates. In the top panel (log density of stars), two key selection effects are visible. First, most of our sample is restricted to R ≲ 15 and z ≲ 8 kpc, a consequence of the relative brightness of our sample with Gaia RVs ( G < 16 . 5). This is also the volume within which parallaxes can provide reliable geometric distances, satisfying our selection of ϖ/σ ϖ > 5. Future work could extend this sample to  \nFigure 2. Summary of our Gaia sample in cylindrical Galactocentric coordinates. Top: Logarithmic density of RGB stars with XP metallicities and Gaia RVs in the Galactocentric R -z plane. The survey selection is strongly mediated by dust extinction near the midplane. Bottom: The R -z plane colored by median XP [Fe/H]. There is a strong negative metallicity gradient with increasing R and | z | .  \n<!-- image -->  \nlarger distances by utilizing the full spectrophotometric information content of these stars (e.g., Hogg et al. 2019).  \nFigure 3 shows the XP-inferred stellar abundance plane, in which the two canonical disk populations highα and lowα - are visible. These definitions are often loosely mapped onto the spatial populations of 'thick' and 'thin' disks respectively. As will be demonstrated in this work (and has already been argued, e.g., Rix & Bovy 2013; Bensby et al. 2014; Hayden et al. 2015), the chemical bimodality maps into spatial differences, but there is no clear bimodality in the vertical structure: the high- and lowα disks differ more in their radial structure than their thickness. Therefore, we use this abundance space to differentiate the two populations, and separately study their spatial and kinematic configuration.  \nThe red lines in Figure 3 show the selection cuts applied to select stars from both disk populations. Stars in the ambiguous region between the two selections are omitted, to ensure that we have pure samples of both highα and lowα disk stars. This is particularly important when studying the highα disk, since even modest  \nFigure 3. Logarithmic density of the all-sky Gaia sample in abundance space, using XP-inferred [ α /Fe] abundances from Li et al. (2023). The red band demarcates our highα and lowα samples, with stars lying in the intermediate region being excluded to keep the samples pure.  \n<!-- image -->  \n[ α /Fe] errors can lead to significant contamination from the much more densely populated lowα disk.", '2.2. 6D Phase Space Coordinates': "For each of RGB star in our sample, the Bayesian prior-informed geometric distance from Bailer-Jones et al. (2021) is adopted, which combines the Gaia parallax with a Galactic space density prior. Our results are not significantly affected if we use zero-point-corrected inverse-parallax distances instead, but the Bailer-Jones et al. (2021) distances provide a more principled framework to infer distances across all SNR regimes.  \nArmed with 6D phase-space information for these RGB stars, we compute their Galactocentric coordinates (Figure 1), azimuthal angular momenta L Z , and total energies E . The latter requires an estimate of the Galactic potential, for which we adopt a four-component Milky Way model containing spherical Hernquist (1990) nucleus and bulge components, an (approximate) exponential disk component (Smith et al. 2015), and a spherical NFW halo component. The potential is fit to the same measurements of the enclosed mass as the MilkyWayPotential in gala (Price-Whelan 2017; PriceWhelan et al. 2020), with the added constraint of having a circular velocity at the solar position = 229 km s -1 (Eilers et al. 2019).  \nUsing the adopted potential, we compute the azimuthal angular momentum L c and total energy E of a perfectly circular orbit. We interpolate the resulting  \nFigure 5 illustrates the formation and evolution of the Milky Way, showing the distribution of stars in angular momentum and metallicity, with metallicity serving as an imperfect proxy for cosmic time. At each [Fe/H], Figure 5 shows the column-normalized distribution of orbital circularity η . The top panel shows all stars in our RGB sample, whereas the bottom panels show the stellar populations differentiated by [ α /Fe] following Figure 3. Although it is tempting to read Figure 5 as a temporal evolutionary sequence from left to right, the age-metallicity relation (AMR) of the Milky  \n<!-- image -->  \n≡  \nFigure 4. Distribution of orbital circularity η for Gaia DR3 giants in three different metallicity selections. Stars at solar metallicity [Fe/H] ≈ 0 are on overwhelmingly prograde and circular orbits with η ∼ 1. At [Fe/H] ≈ -0 . 8 there is predominantly a kinematically hot disk component, along with an isotropic 'in-situ' halo. At [Fe/H] ≈ -1 . 5 the isotropic component dominates, albeit with a prominent tail towards prograde orbits.  \nL c ( E ) curve for each observed star's total energy and calculate the 'orbital circularity' L z /L c . The circularity compactly represents the orbital properties of a star, ranging from -1 (perfectly retrograde, in-plane orbit) to 1 (perfectly prograde, in-plane orbit), with intermediate values corresponding to radial and/or polar orbits. The past literature has used various notations ϵ (e.g., Abadi et al. 2003; Yu et al. 2021, 2023) or η (e.g., Naidu et al. 2020) to describe this parameter, and η ≡ L z /L c is adopted in this work. Figure 4 illustrates the overall circularity distribution of our sample.", '3. THE BIRTH OF THE MILKY WAY IN GAIA DR3': 'Figure 1 shows that these metal-poor stars cluster strongly in the Galactic center, as expected for the old-  \n<!-- image -->  \nFigure 5. The three-phase evolution of the Milky Way revealed by our all-sky Gaia giants. The column-normalized distribution of orbital circularity is shown as a function of Gaia XP metallicity - a proxy for cosmic time. The top panel shows all stars in our sample, whereas the bottom panels are split into highα and lowα subsamples following Figure 3. White lines show the metallicity distribution for each sample. [Fe/H] behaves as a particularly good cosmic clock for the highα sequence in the bottom-left panel, since the age-metallicity relation is reasonably monotonic in this regime.  \n<!-- image -->  \nWay is neither linear, nor monotonic (e.g., Strobel 1991; Bensby et al. 2014; Feuillet et al. 2019; Sahlholdt et al. 2022; Xiang & Rix 2022). However, Xiang & Rix (2022) have demonstrated that stars in the highα disk obey a relatively tight and monotonic AMR. Therefore, the highα panel of Figure 5 can be interpreted as a temporal sequence, whereas the lowα panel is complicated by other effects such as the strong metallicity gradient. In the top panel of Figure 5, three distinct regimes are apparent. In metallicity, the first transition is at -1 . 2 ≲ [Fe/H] ≲ -1 . 0, while the second transition has substantive overlap around -0 . 7 ≲ [Fe/H] ≲ -0 . 5. Starting at low metallicity, we proceed with a discussion of these phases and transitions.', '3.1. The Protogalaxy': "Below [Fe/H] ≈ -1 . 2, the Galaxy is dominated by stars with minimal net rotation. Interestingly, these stars are not entirely isotropic, and have a significantly non-zero median circularity of η ≈ 0 . 1 in both the highα and lowα samples. This population is likely dominated by protogalactic Milky Way stars, in addition to a significant fraction of stars accreted from past mergers like the Gaia-Sausage-Enceladus (GSE; Helmi et al. 2018; Belokurov et al. 2018). This interpretation is supported by the marginal distribution of [ α /Fe] abundances at low metallicity, which for example at [Fe/H] ≈ -1 . 6 peaks sharply at [ α /Fe] ≈ 0 . 3 (protogalactic), with a prominent tail towards [ α /Fe] ≈ 0 . 0 (accreted).  \nest stars formed in the MW (Tumlinson 2010; El-Badry & Rix 2018; Rix et al. 2022; Horta et al. 2023; Semenov et al. 2023a). From Figure 6, the circularity distribution of stars with [Fe/H] < -1 . 0 is asymmetric, with a tail towards prograde orbits (i.e., η > 0). This feature was pointed out by Conroy et al. (2022), who noted a net rotation in the H3 Survey's protogalactic population (also seen in the protogalactic population named 'Aurora' by Belokurov & Kravtsov 2022). Whether or not the protogalaxy had significant net rotation when its components first assembled remains uncertain. McCluskey et al. (2023) show that protogalactic populations with zero net rotational velocity v ϕ at formation could acquire a net positive v ϕ by the present day due to dynamical heating. We explore this possibility in § 4.3 using the TNG50 simulations, affirming this mechanism in a MW analog.", '3.2. Spinup and Cooldown: the Birth of the Disk(s)': "At -1 . 2 ≲ [Fe/H] ≲ -1 . 0, there is a dramatic phase transition, with the emergence of a rotation-dominated disk component with η ≈ 0 . 75, along with significant scatter. This is coincident with the 'spinup' transition identified in v ϕ by Belokurov & Kravtsov (2022); Conroy et al. (2022). It is apparent from the lower panels of Figure 5 that the spinup transition is traced almost exclusively by highα stars, whose AMR appears monotonic (Xiang & Rix 2022): more metal rich means younger. The bottom-left panel of Figure 5 therefore illustrates the birth of the Milky Way's disk over cosmic time, from the turbulent protogalaxy to an increasingly ordered disk. As metallicity (and time) increases, the range of circularities spanned by the highα disk narrows, as stars move on increasingly circular orbits. This 'cooling down' of the orbital kinematics is notably smooth and continuous all the way up to solar metallicity, and reminiscent of the disk 'settling' that is observed in external galaxies and simulations (e.g., Kassin et al. 2012; Simons et al. 2017; Ma et al. 2017; Tiley et al. 2021; Gurvich et al. 2023). Note that for the highα disk, solar metallicity was reached about 8 Gyr ago (Xiang & Rix 2022), and this panel therefore represents the first ≈ 5 Gyr of the MW's evolution.  \nThe top panel of Figure 5 shows another distinct phase transition to yet more circular orbits, at [Fe/H] ≈ -0 . 6. This transition was noted by Belokurov & Kravtsov (2022) as a narrowing of the v ϕ velocity distribution as a function of metallicity. Most stars more metal-rich than this 'cooldown transition' are on very cold and nearcoplanar orbits, with η ≳ 0 . 95 and minimal scatter in η . The behavior of this transition becomes clearer when the sample is split by [ α /Fe]: it is a transition from a  \nhighα dominated disk to a lowα dominated disk. Indeed, the highα disk continues to 'cool down' smoothly from -1 . 0 ≲ [Fe/H] ≲ 0 . 0. In contrast, the lowα disk stars have nearly circular orbits across the full metallicity range [Fe/H] ≳ -0 . 6. Here, it is important to emphasize that among lowα disk stars [Fe/H] is not a proxy for temporal evolution. It mostly reflects the radial metallicity gradient: increasingly more metal-poor lowα disk stars reside at larger radii and angular momenta (Figure 2). In this context [Fe/H] may reflect mostly the stars' birth radius rather than their age (e.g. Frankel et al. 2018, 2020).  \nTo zoom into these later phases in the Galactic disk(s) evolution, it is useful to apply a non-linear transformation to the circularity η , using log (1 -η ) instead, as shown in Figure 6. This 'logarithmic circularity' maps monotonically to η but visually expands the dynamic range among stars on nearly circular and co-planar orbits, making it a more suitable metric to chart out the η ≳ 0 . 9 disk regime. Note that Figure 6 is quite sensitive to the dust-modulated selection effect that limits our sample's depth at low Galactic latitudes (see Figures 1,2). That is, stars with η ≳ 0 . 99 likely exist throughout the cold disk of the MW at very low | b | , but are hidden behind dust-especially at low Galactocentric radii-and are consequently missing from Figure 6.  \nFigure 6 illustrates the markedly different distributions of the highα and lowα disks between -0 . 7 ≲ [Fe/H] ≲ 0 . 0, with the lowα disk residing on kinematically colder orbits throughout this range. As the trends in Figure 6 illustrate, the lowα stars show a sharp break in circularity at [Fe/H] ∼ -0 . 8. This seems to imply that there is no lowα disk below [Fe/H] ∼ -0 . 8; the isotropic ( η ∼ 0) stars in this panel reflect the accreted halo population, which is metal-poor and has low α . As noted above, the highα disk has a steady and log-linear circularity evolution all the way from -1 . 0 ≲ [Fe/H] ≲ 0 . 0.  \nThe kinematic distinctions between the highα and lowα disks are emphasized in the top panel of Figure 6, which displays the L Z angular momentum (left) and logarithmic circularity (right) in the [ α /Fe]-[Fe/H] space for our entire sample. The top-right panel of Figure 6 illustrates how the concepts of 'hot disk' and 'cold disk' map onto the 'highα ' and 'lowα ' disks, respectively. As seen in Figure 6, the highα disk exhibits a gradual cooling down of orbits as metallicity increases, whereas the apparent trend in the lowα disk is primarily set by the radial gradient and selection function. The most circular stars in the outer disk are obscured by dust in the mid-plane, decreasing the median circularity of the outer disk. It will be a worthwhile endeavour for fu-  \nFigure 6. A summary of our kinematic results in the [Fe/H]-[ α /Fe] plane. The top panels show the median logarithmic circularity (left) and L Z angular momentum (right) respectively, and the lower panels show column-normalized profiles of these quantities versus metallicity, split by [ α /Fe] following Figure 3. To guide the eye, a dashed white line shows the characteristic circularity and angular momentum for the lowα disk at solar metallicity.  \n<!-- image -->  \nrk to disentangle these effects by modelling the selection function directly.", '3.3. The Heated Old Disk and In-Situ Halo': "The protogalactic and hot disk components overlap between -1 . 2 ≲ Fe / H ≲ -0 . 9 in Figure 5. In this region, there is also a swath of stars at intermediate orbits between the two populations, that visibly extends up to [Fe/H] ≈ -0 . 5. This entire sequence coincides with the known 'in-situ halo' or 'splash' population (Bonaca et al. 2017; Belokurov et al. 2020). To visualize this portion of the distribution more clearly, we display the row-normalized metallicity-circularity distribution in Figure 7. This plot shows the metallicity dis-  \nion function at each value of orbital circularity. The most circular high-alpha stars ( η > 0 . 75 and [Fe/H] > -0 . 5) exhibit a continuous evolution towards lower metallicities, with a break at [Fe/H] ≈ -0 . 5. Below this metallicity, high-alpha stars occupy a broad swath of orbits. These stars have metallicities identical to the thick (or kinematically hot) disk, but extend all the way to isotropic and even retrograde orbits.  \nThe most plausible origin for these stars is that they were born in the early MW disk and were subsequently kicked up onto isotropic and eccentric orbits by a major merger, presumably GSE itself (Bonaca et al. 2017; Belokurov et al. 2020). Indeed remnant stars from the GSE  \nFigure 7. Row-normalized version of the high-alpha panel from Figure 5-the conditional metallicity distribution for high-alpha stars at each orbital circularity. Three key structures are highlighted: the disk itself, the in-situ halo, and accreted stars from the GSE merger that persist even in the highα sample.  \n<!-- image -->  \n|  \n≡  \nmerger are clearly visible in Figure 7 as an independent swath of isotropic-to-retrograde stars at [Fe/H] ≈ -1 . 2.  \nBelokurov et al. (2020) find that the age distribution of isotropic 'splashed' stars truncates sharply ≈ 9 . 5 Gyr ago, potentially timing the merger's impulse on the old disk. This is consistent with the transition metallicity of [Fe/H] ≈ -0 . 5 seen in Figure 7: stars formed before this metallicity are heated into the in-situ halo, whereas stars formed after this metallicity continue the settling of the old disk towards more circular orbits. In § 4.4 we use TNG50 simulations to lend credence to this hypothesis, finding that the kinematic 'hotness' of both the old disk and the in-situ halo population can be created by a single major merger.", '4. BIRTH OF A MW ANALOG IN TNG50': "In § 3, we used Gaia DR3 to unveil a comprehensive perspective on the structure and kinematics of the Milky Way. We are now in a better position to hunt in large-scale cosmological simulations for evolutionary analogs of the MW, rather than simply structural ones that match its present-day properties. As stated earlier, the eventual goal is to produce a plausible biographical story of our own Galaxy's birth and evolution.  \nIn this section, we search the TNG50 cosmological simulation for analogs of the Milky Way that match the Gaia data presented here. We build on the selections of  \nFigure 8. Mass accretion history for 61 MW-like TNG50 halos (gray lines) selected by S23, showing the virial mass (top), gas mass (middle), and stellar mass (bottom). The 5 early-spinup halos analyzed in this work are highlighted with colored lines.  \n<!-- image -->  \nSemenov et al. (2023a, hereafter S23), who isolated MWlike galaxies in TNG50 that exhibit the 'early spin up' in v ϕ as a function of metallicity (Belokurov & Kravtsov 2022). TNG50 is the highest-resolution box from the IllustrisTNG suite of hydrodynamical cosmological simulations (Pillepich et al. 2018; Nelson et al. 2019; Weinberger et al. 2020). The simulation spans 50 3 co-moving Mpc 3 with a dark matter mass resolution 4 . 5 × 10 5 M ⊙ , and a baryon mass resolution of ∼ 8 . 5 × 10 4 M ⊙ . For more details on the initial conditions, sub-grid physics, and implementation, we defer to the above IllustrisTNG references.  \nOur full sample of simulated MW analogs (S23) consists of 61 star-forming galaxies with halo mass (within the spherical overdensity of 200 × critical density of the universe) of M 200c = 0 . 8-1 . 4 × 10 12 M ⊙ and ˙ M ⋆ > 0 . 2 M ⊙ yr -1 that exhibit rotationally supported stellar disks with v rot /σ v > 6 for stars younger than 100 Myr. By examining the evolution of the azimuthal Galactocentric velocity v ϕ as a function of metallicity, S23 identified 11 TNG50 galaxies that matched the 'early spinup' reported in MW data by Belokurov & Kravtsov (2022). These galaxies all have well-settled, extended disks at z = 0, in contrast with the 'late spin-up' galaxies which are less disk dominated at the present day,  \nFigure 9. The evolution of orbital circularity in five 'early spinup' MW analogs from the TNG50 simulations. Only in-situ stars born within 30 kpc (comoving) of the primary halo are shown. Top: Column-normalized distribution of orbital circularity as a function of metallicity. Note the similarity between TNG50 galaxy #519311 and the MW data shown in Figure 5. Bottom: Column-normalized distribution of orbital circularity as a function of stellar age. Vertical lines denote snapshots containing relatively major mergers, where the accreted stellar mass increases by > 5%.  \n<!-- image -->  \nwith less ordered orbits. S23 report that the key evolutionary feature distinguishing these populations is that the early spin-up galaxies undergo more rapid mass assembly at early times. We note that although S23 used a CDF-matching approach to rescale their simulation metallicities to the Belokurov & Kravtsov (2022) data, our Gaia data are subject to a very different selection function. We therefore directly use the metallicities reported by the simulation, which are sufficient for the qualitative analysis presented here. Figure 8 shows the dark matter, gas, and stellar mass accretion history of five of the galaxies that spun up earliest, which are the ones we analyze further in this work.", '4.1. The Evolution of Spin in TNG50 Disks': "The kinematic evolution of our five early-spinup TNG50 galaxies is shown in Figure 9. Only 'in-situ' stars, those born within 30 kpc (comoving) of the galaxy, are shown. The top row corresponds to our Figure 5, showing the column-normalized distribution of circularity η versus metallicity. Remarkably, the galaxy in halo #519311 exhibits the same distinctive three-phase behaviour seen in our MW data (Figure 5). We have verified that #519311 is the most similar to the Milky Way's η -[Fe/H] data among all 61 MW-mass halos identified by S23.  \nThe bottom row in Figure 9 provides the much-needed dimension that our data in Figure 5 are lacking: precise stellar ages. These circularity diagrams, η (Age) can  \ngenuinely be interpreted as evolutionary sequences that show what orbits stars born at different epochs end up on. We emphasize that η here denotes the orbital properties at z = 0, not at birth. As a consequence of the selection in S23, all these galaxies 'spin up' more than 10 Gyr ago, matching MW data. Several sharp spikes towards isotropic orbits are visible in most halos, which coincide with merger events. Dashed lines indicate 'major merger' events, which are selected to increase the fraction of accreted stellar mass in the galaxy by > 5%. Although this metric is imperfect-as evidenced by smaller mergers or gas accretion events without preceding dashed lines-it is sufficient to highlight the most prominent mergers in a galaxy's history.  \nInterestingly, in the ageη space, the three-stage picture of disk evolution (protogalaxy, then old hot disk, then young cold disk) is evident for all halos except #527309, while in the [Fe/H]η space, the distinction between old hot and young cold disk is less apparent. This is a consequence of a significant overlap in the chemical evolution of the two disks (see Section 3.2), implying that, while the three-phase formation picture can be common, its imprint on the chemo-kinematic distribution of stars is less common. The heated old disk population is typically superseded by the cold, young disk across all metallicities. However, in halo #519311 there is enough of a metallicity difference between the disks for the old disk to remain visible in the [Fe/H]-  \nFigure 10. The evolution of TNG50 halo #519311 over cosmic time. 5 key snapshots are temporally ordered from left to right, with the lookback time and redshift indicated at the top. The top panel shows the orbital circularity-metallicity space for all in-situ stars born up to that time, whereas the second panel shows the [Fe/H]-[ α /Fe] plane, with stars younger than 100 Myr in each snapshot highlighted. The orbital circularities are computed in a coordinate frame aligned with the net angular momentum of each snapshot. The third panel shows the logarithmic density of star particles in fixed coordinates oriented edge-on to the present-day disk, and the fourth panel shows the same for gas particles. This halo undergoes a major merger at z ≈ 0 . 7, visible in the central panels. An animated version of this figure is available here.  \n<!-- image -->  \nη space at the present day, matching the three-phase picture observed in the MW.", '4.2. Halo #519311: An Illustrative Milky Way Analog': "As mentioned above, TNG50 halo #519311 exhibits a 'three-phase' structure in the circularity-metallicity space that is remarkably similar to the MW data (compare Figure 9 to Figure 5). This halo is one of the earliest-spinup galaxies in the S23 sample, and also experiences some of the earliest mass growth in the entire TNG50 MW analog sample (Figure 8). This particular halo also experiences a uniquely gas-rich major merger ≈ 8 Gyr ago (bottom panel of Figure 8). Although this merger is slightly later than the hypothesized timing of the Milky Way's own GSE merger (Bonaca et al. 2020; Belokurov et al. 2020; Naidu et al. 2021), the fact that #519311 is the only MW-like TNG50 galaxy to  \nshow the distinctive three-phase kinematic behavior in the circularity-metallicity space warrants further investigation. In this section, we take a closer look at halo #519311 as a Milky Way analog.  \nFigure 10 illustrates the evolutionary history of TNG50 halo #519311. Although the η (Age) shown in the bottom panels of Figure 9 helps to sort the presentday orbits of stars, it lacks information about the birth orbits. In the top row, each panel shows what the η ( z | [Fe / H]) distribution would have looked like at different epochs. The orbital circularities η ( z ) are computed at that epoch for all in-situ stars born to that point. The second row illustrates the [Fe/H]-[ α /Fe] abundance plane for all stars born up to that point (grey) and for young stars born within 100 Myr of the chosen snapshot. The lower two rows show the distribution of stars and gas in each snapshot, in a fixed coordinate frame cho-  \nto be the 'edge-on' orientation of the disk at z = 0. Note that the circularity η ( z ) is not computed in this fixed frame, but in a frame aligned with the net angular momentum vector of the galaxy in each snapshot. Between the second and third panels from the left, this halo experiences a major merger at z ≈ 0 . 7 ( τ ≈ 7 Gyr). The accreted dwarf galaxy is visible in the stellar and gas images in the central panel.  \nAt z = 4 . 2, a billion years after the Big Bang, a nascent rotating population of relatively high [Fe/H] has started to succeed the isotropic protogalaxy. Yet, as the images in Figure 10 show, the galaxy is still spheroidal and clumpy, fed by numerous chaotic streams of gas. By z = 2, the spinup transition is complete and a remarkably thin, kinematically coherent disk has already formed. At the time of its formation, this old disk was kinematically much colder than it is now, with η ≳ 0 . 9 across a wide range of metallicities. This suggests a significant amount of kinematic heating, the source of which is immediately apparent in the next snapshot at z ≈ 0 . 7: the arrival of a major gas-rich merger.  \nThe major merger is visible in the central panels of Figure 10, both in the distribution of accreted stars and in the images of the stellar and gas density. Several interesting features can be gleaned from these panels alone. We emphasize that this is an exceptionally major merger among the early-spinup TNG50 galaxies we analyze. As Figure 8 illustrates, among all our MW analogs, halo #519311 experiences one of the most massive and gas-rich mergers that does not destroy the preexisting stellar disk. At the time of the merger, the distribution of circularities in the old disk thickens considerably. This is not only due to orbital heating by the merger but also due to a rotation of the disk itself: as can be seen by comparing the stellar images, the galaxy has tipped by ≈ 30 · immediately after the merger. Indeed, we find that the merger and subsequent gas accretion contribute significant angular momentum to the galaxy, setting the final orientation of the galactic disk.  \nA comparison between the stellar and gas distributions in the central panel of Figure 10 also reveals that the accreted gas lags behind the accreted stars by several kpc. Although stars from the accreted dwarf have already reached the disk of the primary galaxy, the bulk of the accreted gas is still ≈ 10 kpc above the disk. Although a precise investigation of this offset is warranted, a likely explanation is tidal stripping, along with ram pressure stripping of the accreted dwarf by the circumgalactic medium of the primary galaxy (e.g., Simpson et al. 2018; Rohr et al. 2023). This additional drag force on the gas would naturally cause it to lag behind the collisionless stellar component.  \nMoving forward to one Gyr after the merger, the galaxy remains in a state of significant disequilibrium. Although the stellar component of the accreted dwarf galaxy has started isotropizing into a spheroidal 'in-situ halo,' there are still turbulent flows of gas accreting onto the disk. The rotation and torquing of the galaxy continues, now over ≈ 60 · from the initial orientation of the old disk. This inflow of gas also increases the gas density in the disk, triggering an intense period of star formation. This creates stars on exceptionally circular orbits starting at [Fe/H] ≈ -0 . 5. The z = 1 . 0 snapshot demonstrates that the old disk rapidly enriches to solar metallicities. Subsequently (see the z = 0 . 5 snapshot), it becomes superseded by the young disk at [Fe/H] ≳ -0 . 5, creating the characteristic 'three-phase' appearance of the metallicity-circularity plot.  \nThe last 5 Gyr of this galaxy's evolution are relatively quiescent, with the rest of the kinematically cold young disk assembling and enriching up to super-solar metallicities. The final orientation of the young disk is misaligned from the original old disk by almost 90 · . The final metallicity-circularity distribution of this galaxy is remarkably similar to the MW data shown in Figure 5, a unique feature among all the early-spinup TNG50 galaxies (Figure 9). This feature is rather unusual in the [Fe/H]η plane likely because of the relatively high metallicity of the young disk ([Fe/H] ≳ -0 . 5), such that the remnants of the old disk with [Fe/H] ≲ -0 . 5 remain exposed in this plane. In the other examples in Figure 9, the old disk is buried (in number counts) under the low-[Fe/H] end of the young disk.", '4.3. The Protogalaxy: Past and Present': "Before the galaxy in TNG50 halo #519311 begins the spinup transition to form the old disk, it assembles a sizable spheroidal 'protogalaxy' that already contains 3 × 10 8 M ⊙ of stellar mass by z ≈ 6. Rix et al. (2022) utilized all-sky Gaia data to unveil a substantial population of metal-poor stars residing in the heart of the MW, and used rough scaling arguments to estimate that the MW assembled ≳ 5 × 10 7 M ⊙ in stars by z ≳ 5. Belokurov & Kravtsov (2023) also estimate a stellar mass for the protogalactic 'Aurora' population, finding 5 × 10 8 M ⊙ . This population was also studied by Belokurov & Kravtsov (2022) and Conroy et al. (2022), who found that although this ancient ' in-situ ' population is predominantly isotropic, it exhibits a small amount of net prograde rotation. These results were corroborated in the high-resolution FIRE simulations by Horta et al. (2023), who found that protogalactic populations in all their MW analogs show weak net rotation aligned with the present-day disk.  \nFigure 11. Top: The circularity distribution of all protogalactic stars in TNG50 halo #519311 born up to z = 6 (red), and the circularity distribution of those same stars at the present day (blue). Bottom: Edge-on spatial distribution of protogalactic stars at the present day.  \n<!-- image -->  \nFigure 11 shows the protogalactic stellar population assembled by z = 6 in TNG50 halo #519311. The top panel shows the circularity distribution at z = 6 (red) and - for the same stars - at z = 0 (blue). Again, these circularities are computed relative to the net angular momentum of each snapshot, not relative to the orientation of the present-day disk. Although the kinematics of these protogalactic stars were completely disordered and isotropic in the pre-disk era, they have acquired a prominent tail of prograde orbits by the present. Therefore, the protogalaxy itself experiences 'spin-up' in the true sense of the word, with protogalactic stars being perturbed by the growing disk, bar, or merger around it to preferentially have mildly prograde orbits (e.g., Dillamore et al. 2022). The bottom panel of Figure 11 shows the present-day edge-on distribution of protogalactic ( z > 6) stars. The distribution is visibly oblate, with an ellipticity ϵ ≈ 0 . 7 at z = 0, as compared to ϵ ≈ 0 . 3 at z = 6, matching results from FIRE (Horta et al. 2023).", '4.4. The Highα Disk and In-Situ Halo: Remnants of the Old Disk': "It is apparent from Figures 5 and 6 that the Milky Way's highα disk is kinematically hotter than the younger lowα disk. Over time, secular orbital interactions and mergers are expected to vertically 'heat' stars to kinematically hotter orbits (e.g., Quinn et al. 1993;  \nFigure 12. Top: Circularity distribution for old disk stars selected to have η > 0 . 9 at z = 1 (red), and the distribution for those same stars at the present day (blue). Note that the distribution is logarithmic, so the present day distribution is still disky, albeit with a long tail towards isotropic and retrograde orbits. Bottom: The spatial distribution of the old disk stars at z = 1 (red borders) and at z = 0 (blue borders). The first image is in an edge-on frame relative to the z = 1 galaxy, and the other two are edge-on relative to the z = 0 galaxy.  \n<!-- image -->  \nFont et al. 2001; Kazantzidis et al. 2008; Villalobos & Helmi 2008). It therefore remains an open question to what extent stars in the highα disk were heated after formation, as opposed to, say, being born in more turbulent gas at early times. This question was extensively explored in the high-resolution FIRE simulations by Yu et al. (2023), who found that to first order, stars remain on orbits similar to their birth orbits. Mergers and heating processes only play a secondary role in the  \n12 MW analogs analyzed by Yu et al. (2023), although those particular galaxies have relatively quiescent merging histories over the past 10 Gyr.  \nFigure 12 explores the evolution of the old disk in the TNG50 halo #519311. We identify the 'old disk' as all stars on nearly circular orbits ( η > 0 . 9) that have formed by z = 1, just prior to the last major merger (red dashed histogram in the top panel). Around half of the z = 1 stellar mass of the galaxy is in this highly circular disk component. As the images in the lower panels illustrate, these stars form a coherent and relatively thin disk in the z = 1 'edge-on' frame. However, in the z = 0 frame-in which the present-day disk would be edge-on-this old disk is nearly face-on (see also Figure 10).  \nMoving forward to the present day, these old disk stars have spread out over a large range of orbital circularities, extending all the way down to the retrograde regime (blue histogram in the top panel of Figure 12). The circularity distribution is reminiscent of the observed present-day 'old disk' in the Milky Way (Figures 4 and 5). The spatial distribution resembles a thick disk, ensconced in an oblate halo of scattered stars. The old disk constitutes ≈ 10% of the present-day stellar mass of this galaxy, comparable to the mass fraction of the MW's 'thick' disk (Juri'c et al. 2008; Bland-Hawthorn & Gerhard 2016).  \nThis illustrative example demonstrates how a presentday thick disk and in-situ halo population can be formed by a coherent and thin old disk that is initially misaligned, and significantly torqued by the arrival of a major merger (see also Renaud et al. 2021a; Meng & Gnedin 2021; Dillamore et al. 2022; Dodge et al. 2023). Given observational evidence that the Milky Way spun up its old disk remarkably early, and that it underwent a major merger soon afterwards, this kind of disk tilting is a plausible formation channel for the thick disk in our own Galaxy (Bonaca et al. 2017; Bonaca et al. 2020; Belokurov et al. 2020; Belokurov & Kravtsov 2022; Semenov et al. 2023a).", '5. DISCUSSION': "In the beginning, the nascent Milky Way existed as a chaotic and disordered 'protogalaxy', constantly bombarded by mergers that built up its most ancient stars (e.g., Zolotov et al. 2009; Tumlinson 2010; Pillepich et al. 2015; El-Badry et al. 2018; Belokurov & Kravtsov 2022; Semenov et al. 2023a; Horta et al. 2023). Rix et al. (2022) used an all-sky sample of giants with Gaia lowresolution spectra to reveal this metal-poor population residing in the heart of the MW (building on a rich legacy of observational efforts, e.g., Garc'ıa P'erez et al. 2013; Ness et al. 2013; Schlaufman & Casey 2014; Casey  \n& Schlaufman 2015; Ness et al. 2015; Koch et al. 2016; Garc'ıa P'erez et al. 2018; Lucey et al. 2019; Arentsen et al. 2020a,b; Reggiani et al. 2020; Arentsen et al. 2021; Lucey et al. 2021, 2022). This protogalactic population is prominent in Figure 1, and Figure 4 illustrates how protogalactic stars are net isotropic, albeit with a tail towards prograde orbits. § 4.3 and Figure 11 demonstrate that a protogalactic population in the TNG50 simulations can acquire such a prograde tail after birth, even if it was initially perfectly isotropic (see McCluskey et al. 2023; Horta et al. 2023).  \nWhat followed was the most drastic phase transition in our Galaxy's history, the birth of the Galactic disk. The ancient disk coincides with the present day 'highα ' (also, hot or thick) disk, although there is an imperfect and imprecise correspondence between these various classifications. The highα disk could be 'thick' due to a combination of thicker birth orbits (e.g., due to a more turbulent ISM at high redshift), and subsequent heating by gravitational interactions, with the relative importance of these mechanisms being an active area of research (Forster Schreiber et al. 2009; Dekel et al. 2009a; Bird et al. 2013; Stott et al. 2016; Bird et al. 2021; van Donkelaar et al. 2022; Yu et al. 2023; McCluskey et al. 2023; Hamilton-Campos et al. 2023). Here differences between the adopted physics across simulations also become relevant. For example, galaxies in the FIRE simulations experience stronger feedback and more 'bursty' star formation than TNG (e.g., Gurvich et al. 2023). The close MW analog that we find in TNG50 (#519311) forms a kinematically-ordered albeit misaligned disk prior to its last major merger, which is subsequently torqued and 'thickened' by the merger (Figures 10 and 12). Early results from the first year of JWST observations have revealed a large population of disks at z ≳ 3 -6, supporting the idea that well-ordered disks can assemble in the early Universe (Ferreira et al. 2022; Robertson et al. 2023; de Graaff et al. 2023; Jacobs et al. 2023).  \nAround z ≈ 1, star formation from α -enhanced birth material came to an end, yielding to the lowα (or 'thin') disk. The precise timing of these transitions is uncertain and likely overlapping, with the characteristic observed transition time being around ≈ 8 Gyr ago (Bonaca et al. 2020; Xiang & Rix 2022). This transition timescale coincides with a 'reset' in the metallicity of the disk from [Fe/H] ≈ 0 to [Fe/H] ≈ -0 . 8 at roughly constant [ α /Fe] ≈ 0 . 1, creating the 'bimodal' abundance distribution of the disks (see Figure 10 in Haywood et al. 2013, and also Bensby et al. 2014; Hayden et al. 2015; Martig et al. 2016; Nissen et al. 2020; Conroy et al. 2022). Although it remains to be seen how common such a chem-  \nty is in observed Milky-Way-like galaxies (e.g., Nidever et al. 2023; van de Sande et al. 2023), spatially distinct disks (thick and thin) are known to be ubiquitous, with the prevalence of old thick disks suggesting a merger-related origin (e.g., Dalcanton & Bernstein 2002).  \nThe Milky Way's age pattern in the [Fe/H] - [ α /Fe] abundance space is suggestive of a 'zig-zag' evolution, with the metallicity 'resetting' from the solar-metallicity highα disk to the metal-poor lowα disk, at nearconstant [ α /Fe] (Haywood et al. 2013; Bensby et al. 2014; Conroy et al. 2022). Indeed, recent high-resolution simulations of MW analogs have demonstrated how a major gas-rich merger can produce the chemical bimodality and age trend seen in the Milky Way (Buck et al. 2020; Buck 2020; Agertz et al. 2021; Renaud et al. 2021b,a). However, it has been challenging to associate the GSE merger definitively with the bimodality of the disks, primarily to the lack of precise ages. Different studies have arrived at disparate conclusions about the temporal ordering of the truncation of the highα disk, arrival of the GSE merger, and birth of the lowα disk (Gallart et al. 2019; Bonaca et al. 2020; Xiang & Rix 2022). Improvements in the precision and accuracy of stellar ages will be vital to produce a more fine-grained picture of the early Milky Way.  \nMyeong et al. (2022) find evidence for a new eccentric component of the MW halo, which they name Eos. Eos stars lie broadly in the lowα sequence around [Fe/H] ≈ -1 . 0. At low metallicities, Eos chemically resembles the metal-rich accreted stars from the GSE merger, whereas at higher metallicity Eos connects with the lowα disk. This population therefore exhibits features that we would expect of the earliest lowα disk stars born soon after the GSE merger, immediately following the metallicity 'reset' of the disk. This possibility warrants further investigation, chiefly by agedating Eos stars. Furthermore, it remains to be seen whether detailed hydrodynamic simulations of the MWGSE merger can reproduce such a population in the outer disk, and whether the kinematic properties of Eos match these predictions (see Orkney et al. 2022; Buck et al. 2023; Rey et al. 2023).  \nRegardless of whether a major merger is required to produce the chemical bimodality of the MW disk, it is now quite certain that the MW did undergo such a merger at z ≈ 2. One of the most significant results from the Gaia era was the confirmation that the majority of the MW's stellar halo was likely deposited by a single major merger, the Gaia -Sausage Enceladus (GSE, e.g., Belokurov et al. 2018; Helmi et al. 2018; Haywood et al. 2018; Gallart et al. 2019; Mackereth et al. 2019;  \nNaidu et al. 2020, 2021; Myeong et al. 2022; Chandra et al. 2023). Bonaca et al. (2020) and Giribaldi & Smiljanic (2023) find that star formation in GSE truncated 9 . 5 -10 Gyr ago, preceding the formation of the lowα disk by ≈ 2 Gyr. This merger could also have dynamically perturbed the ancient disk, creating the observed population of stars on halo-like orbits with chemistry that is indistinguishable from the highα disk (the insitu halo or 'splash', Bonaca et al. 2020; Belokurov et al. 2020, see Figure 7). It is tempting to imagine a unified picture in which the GSE merger mediates the transition between the highα and lowα disks, providing the pristine gas required to create a chemical bimodality, and dynamically creating the in-situ halo in the process. Such a scenario is thus far supported-but not demanded-by current MW data.", '6. CONCLUSIONS': "We have assembled a dataset of 10 million red giant stars with 6D kinematics, metallicities, and α -abundances using data from Gaia DR3 (Figures 1-4), and have used it to draw a comprehensive empirical picture of our Milky Way's formation history. This picture builds on a rich legacy of prior work in Galactic archaeology, now synthesized with all-sky Gaia data:  \n- 1. The angular momentum evolution of the Milky Way exhibits three distinct phases as a function of metallicity: the disordered protogalaxy ([Fe/H] ≲ -1 . 2), the now-hot old disk ( -1 . 0 ≲ [Fe/H] ≲ -0 . 7), and the cold young disk ([Fe/H] ≳ -0 . 5). The transitions between these phases are sharp ( ∼ 0.2 dex in [Fe/H]), and can be described as 'spinup' and 'cooldown' respectively. The hot old disk and cold young disk are distinctly separable in [ α /Fe]-[Fe/H] abundance space, mapping onto the highα and lowα disks respectively (Figure 5).\n- 2. The observed protogalaxy is nearly isotropic and mostly confined to the inner heart of the Galaxy, but has mild net rotation even at rather low metallicities (Figures 4 and 5). More detailed individual chemical abundances are required to perform a precise excision of the accreted halo population from the bona-fide protogalaxy in the inner Galaxy.\n- 3. Since the old, highα disk is known to have a relatively monotonic age-metallicity relation, we construct a time-ordered portrait of the birth of the Milky Way's ancient disk (Figure 5, bottomleft). After spinning up from the protogalaxy, the highα disk exhibits a smooth evolution toward  \n- increasingly circular orbits throughout its metallicity range -1 . 0 ≲ [Fe/H] ≲ 0 . 0 (Figure 6). However, at [Fe/H] ≲ -0 . 5, highα stars also occupy a vast swath of orbits extending all the way to the retrograde regime, forming an in-situ halo population (Figure 7).\n- 4. The young, lowα disk occupies highly circular orbits across all metallicities and exhibits a wellestablished radial metallicity gradient. The effects of this gradient are best visualized in the bottomleft panels of Figure 6, which illustrate the strikingly different angular momentum profiles of the high- and lowα disks as a function of metallicity.  \nIn § 4, we interpret these observed trends using MW analogs selected from the IllustrisTNG50 cosmological simulations. We emphasize the importance of finding evolutionary analogs to the Milky Way, instead of only structural analogs at the present day. Our key findings are as follows:  \n- 5. S23 find that out of 61 MW-like galaxies in TNG50, only 10 exhibit an 'early spinup' in rotational velocities that match MW data. We find that although most of these early spinup analogs exhibit a three-stage evolution in the circularityage space, only one (halo #519311) matches our observations of a three-phase behavior in the circularitymetallicity plane (Figure 9). This galaxy is notable for a particularly gas-rich major merger ≈ 7 Gyr ago-reminiscent of (but probably later than) the MW's own merger with GSE-and a relatively high metallicity of the young cold disk (Section 4.2, Figure 8).\n- 6. The time-resolved evolution of this simulated MW analog illustrates that it underwent an early spinup ≈ 11 Gyr ago, and had a kinematically ordered disk well in place by ≈ 8 Gyr ago. The subsequent gas-rich merger tilted and heated the ancient disk, and a new kinematically cold disk subsequently formed, leading to the 'three-phase' kinematic signature (Figure 10).\n- 7. The simulated MW analog exemplifies that protogalactic populations can flatten and acquire net rotation due to external forces from the surrounding disk, even if the protogalaxy was originally completely isotropic (Figure 11). It also demonstrates how an initially misaligned, kinematically cold old disk can be tilted and heated by a major merger, creating both a spatially thick disk and in-situ halo component (Figure 12).\n- 8. Overall, the MW analog from TNG50 demonstrates how a single decisive major merger striking an early-formed galactic disk can reproduce various kinematic properties of the Milky Way. Such a merger-driven formation scenario for the disks is therefore strongly supported-but not yet demanded-by the observed data.  \nOur work demonstrates that the three dominant evolutionary stages of the Milky Way-the protogalaxy, old disk, and young disk-are inextricably linked by the three physical processes of spinup, merger, and cooldown. Whereas there is now significant evidence from the halo that our Galaxy underwent a major merger ≈ 10 Gyr ago, here we show that this merger could have played a decisive role in transforming the Galactic disks as well. Although several pieces of the puzzle are still missing-chiefly precise stellar ages that could support the ordering and causal links between phases-this simple picture does a remarkable job reproducing the overall kinematic structure of the Milky Way. A concert of future surveys and simulations will hopefully further illuminate a plausible biography of the Galaxy we call home.  \nVC gratefully acknowledges a Peirce Fellowship from Harvard University. It's a pleasure to thank Angus Beane, Jonathan Bird, James Bullock, Kareem ElBadry, Michael Hayden, David W. Hogg, Phil Hopkins, Adrian Price-Whelan, Sandro Tacchella, and Turner Woody for insightful conversations and feedback. This work benefited from discussions at the 2022 Disk Formation Workshop at UC Irvine, the 2023 Gaia XPloration workshop at the University of Cambridge, and the 2023 Galaxy Formation × Theory Workshop at Biosphere 2. Support for V.S. was provided by Harvard University through the Institute for Theory and Computation Fellowship. RPN acknowledges support for this work provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51515.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. LH acknowledges support by the Simons Collaboration on 'Learning the Universe'.  \nThis work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www. cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University. This research has made extensive use of NASA's Astrophysics Data System Bibliographic Services.  \nSoftware: numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), gala (PriceWhelan 2017; Price-Whelan et al. 2020)  \nFacilities: Gaia, WISE", 'REFERENCES': 'Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 597, 21, doi: 10.1086/378316 Agertz, O., & Kravtsov, A. V. 2015, ApJ, 804, 18, doi: 10.1088/0004-637X/804/1/18 -. 2016, ApJ, 824, 79, doi: 10.3847/0004-637X/824/2/79 Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826, doi: 10.1093/mnras/stab322 Andrae, R., Rix, H.-W., & Chandra, V. 2023, arXiv, arXiv:2302.02611, doi: 10.48550/arXiv.2302.02611  \nArentsen, A., Starkenburg, E., Martin, N. F., et al. 2020a,  \nMNRAS, 491, L11, doi: 10.1093/mnrasl/slz156  \n-. 2020b, MNRAS, 496, 4964,  \ndoi: 10.1093/mnras/staa1661  \nArentsen, A., Starkenburg, E., Aguado, D. S., et al. 2021, MNRAS, 505, 1239, doi: 10.1093/mnras/stab1343  \nBailer-Jones, C. A. L., Rybizki, J., Fouesneau, M.,  \nDemleitner, M., & Andrae, R. 2021, The Astronomical  \nJournal, 161, 147, doi: 10.3847/1538-3881/abd806'}