spheroids (Oort 1922 ;Lindblad 1925 ;de Vaucouleurs 1959 ;Burstein 1979 ;Tsikoudi 1979 ;Yoachim & Dalcanton 2006 ;van der Kruit & Freeman 2011 ).Thin discs are systematically younger and more metal-rich than both spheroids and thick discs (Yoachim & Dalcanton 2008 ;Comer ón et al. 2011 ). By spheroid, we refer quasi-isotropic populations, which could include inner stellar haloes (Carollo et al. 2007 ;Bonaca et al. 2017 ) and isotropic bulge components. Note that bulge components in galaxies are usually classified in one of two w ays: flattened, mildly rotating pseudo-bulg es or more spheroidal classical bulges (Kormendy & Kennicutt 2004 ;Gao et al. 2020 ). Most stars in classical bulges appear to be formed at high redshift, with only minor growth at late times (Renzini 1999 ).

Among the most important ideas to emerge in galaxy formation theory is that rotationally supported disc galaxies form naturally as a result of angular momentum conservation, with the angular momentum source driven by gravitational collapse in an expanding universe (Peebles 1969 ;Fall 1979 ). Rotating galaxies may then be seen as a natural starting point for subsequent morphological evolution (Hodge et al. 2019 ;Rizzo et al. 2020 ;Kretschmer, Dekel & Teyssier 2022 ;Tamfal et al. 2022 ). In the simplest version of this scenario, thin discs form early and continuously, with thick discs and bulges arising only from secular evolution and/or mergers. A more nuanced version of this idea, moti v ated by CDM galaxy formation simulations, is that the first discs to emerge in the early universe are thick and turbulent, with thin-disc galaxies developing only late in cosmic assembly (Brook et al. 2012 ;Wuyts et al. 2012 ;Bird et al. 2013Bird et al. , 2021 ; ;Park et al. 2019 ;Yu et al. 2021 ; Se go via Otero, Renaud & Agertz 2022 ).

Observationally, disc galaxies at higher redshift do appear to be more disordered and clumpy (Elmegreen et al. 2007 ;Shapiro et al. 2008 ;Genzel et al. 2008 ;Overzier et al. 2010 ;Elmegreen et al. 2017 ;Osborne et al. 2020 ) and only later does star formation begin to occur primarily in extended thin discs (Kassin et al. 2012 ). This behavior is consistent with a picture where the discs evolve in quasi-stable equilibrium, with higher ISM velocity dispersion at early times owing to an increasing gas fraction (Faucher-Gigu ère, Quataert & Hopkins 2013 ; Wisnioski et al. 2015 ;Ceverino et al. 2017 ;Bird et al. 2021 ). Interestingly, some of the first observations with JWST suggest that discs may be more common in the early universe than previously believed (Ferreira et al. 2022 ;Robertson et al. 2022 ).

The rotational and structural properties of pseudo-bulges are indicative of a formation channel linked closely to disc formation and/or disc/bar evolution (Okamoto 2013 ;Kormendy 2015 ;Spiegel & Polyachenko 2019 ;Devergne et al. 2020 ). Classical bulges, on the other hand, appear to be less connected to disc properties and it remains unclear if their formation is linked to or completely disjoint from disc assembly. Traditional scenarios posit that classical bulges emerge from the mergers of discs (van den Bosch 1998 ; Bois et al. 2010 ;Hopkins et al. 2010 ;Kannan et al. 2015 ), though there is theoretical and observ ational e vidence that bulges can form in multiple ways (e.g. Obreja et al. 2013 ;Seidel et al. 2015 ). Interestingly, direct measurements suggest that bulge mass is not linked closely to galaxy merger history at the MW scale (Bell et al. 2017 ). Secular process, like rapid gas inflow to the galaxy center (Scannapieco et al. 2009 ), or the migration of giant gaseous clumps (Minchev & F amae y 2010 ;Cev erino et al. 2015 ), could naturally produce quasi-isotropic bulge formation.

Similarly, the inner stellar haloes of galaxies may be populated by multiple channels. Today there is general agreement that accretion is responsible for the majority of stars in the outer stellar haloes of galaxies (Bullock, Kravtsov & Weinberg 2001 ;Bullock & Johnston 2005 ;Bell et al. 2008 ;De Lucia & Helmi 2008 ;McConnachie et al. 2009 ), though a fraction could originate in outflows from the main galaxy (Yu et al. 2020 ). Inner stellar halo, on the other hand, is likely populated by both accreted stars (e.g. Helmi et al. 2018 ;Mackereth et al. 2019 ;Simion, Belokurov & Koposov 2019 ) and stars that were born within the inner galaxy (Carollo et al. 2007 ;Cooper et al. 2015 ). This in situ inner stellar halo is often believed to consist of stars that were born on orbits initially confined them to a disc but became heated to eccentric orbits by mergers (Zolotov et al. 2009 ;Purcell, Bullock & Kazantzidis 2010 ) or by potential fluctuations from e xplosiv e feedback ev ents (El-Badry et al. 2018 ).

Alternatively, inner halo or spheroid formation could arise from an earlier, very turbulent phase of galaxy assembly with no need for coherent rotation. Interestingly, this possibility bares considerable resemblance to a conjecture put forth by Larson ( 1976 ), who used numerical experiments to suggest that spheroidal/bulge components arise in an early stage of rapid star formation, while the formation of discs requires a later stage of slower star formation to allow the residual gas to settle into a disc before forming stars (see also Eggen, Lynden-Bell & Sandage 1962 ;Gott & Thuan 1976 ). Larson concluded that 'Spheroidal systems may, thus form from gas which experiences strong turbulence or cloud collisions, and disc systems may form from more quiescent residual gas in which collisions are less important'.

The MW provides a detailed touchstone for testing these ideas. Our Galaxy has both thin and thick disc components (Gilmore & Reid 1983 ;Juri ć et al. 2008 ;Bensby et al. 2011 ;Bovy & Rix 2013 ;Hayden et al. 2015 ). Its inner stellar halo appears to be populated by both in situ and accreted components (Carollo et al. 2007 ;Belokurov et al. 2018 ;Helmi et al. 2018 ). The Galaxy has a small classical bulge (Kunder et al. 2016 ) and a more dominant pseudo-bulge (Gonzalez & Gadotti 2016 ), which appears to resemble the thick disc in chemical properties and formation time (Alves-Brito et al. 2010 ;Haywood et al. 2018 ;Di Matteo et al. 2019 ).

Recently, chemo-dynamical data sets have started to uncover exciting clues to the origin of the rotating components of the Galactic disc. Belokurov & Kravtsov ( 2022 ) used APOGEE and Gaia data to identify a characteristic 'spin-up' metallicity for MW stars. Specifically, in situ stars show a rapid increase in net rotation as a function of metallicity at [Fe/H] -1, from a median tangential velocity of ∼0 km s -1 (typical of a spheroid) to ∼100 km s -1 (typical of a thick disc). This feature may point to transition epoch from disordered kinematics to increasingly coherent rotation. Similarly, Conroy et al. ( 2022 ) have used H3 Survey spectroscopy and Gaia astrometry to identify a transition time where star formation efficiency rapidly increased while simultaneously stellar kinematics become more disclike. Stars that formed before this time retain an isotropic velocity distribution. There are features of this picture, where there is an early 'spin-up' phase that follows a less well-ordered phase, that are quite similar to those seen in cosmological simulations of disc galaxy formation (e.g. Park et al. 2021 ).

In what follows, we use FIRE-2 (Hopkins et al. 2018 ) simulations to study the orbital properties of stars formed in the main progenitors of MW size galaxies, from birth to the present day, in order to gain insight into the origin of thin discs, thick discs, and in situ isotropic spheroids. It extends work from a series of FIRE-2 papers that hav e e xamined the relationship between star formation burstiness, g alaxy kinematics, g alaxy morphology, metallicity gradients, and the development of a hot gaseous haloes around galaxies (Ma et al. 2017 ;Stern et al. 2021 ;Yu et al. 2021 ;Bellardini et al. 2022 ;Gurvich et al. 2023 ;Hafen et al. 2022 ). In particular, these simulations have revealed a correlation between internal galaxy properties and the mode of gas deposition from the circum-galactic medium (CGM) into the interstellar medium (ISM) (Stern et al. 2021 ). As galaxy haloes evolve from low mass to high, the inner CGM virialization, star formation transitions from 'bursty' to 'steady', and stellar-driven galaxy-scale outflows are suppressed. The bulk of thick-disc stars form prior to this transition, and this gives rise to a tight correlation between the ages of thick-disc stars and the end of the bursty phase (Yu et al. 2021 ). CGM virialization also drives an abrupt change in the angular momentum coherence of accreting gas (Hafen et al. 2022 ). Only after this time do stars form along a single long-lived plane in circular orbits, making possible the formation of a thin disc (Yu et al. 2021 ;Hafen et al. 2022 ). Conversely, during the earliest epochs, the ISM has a quasi-spheroidal morphology, and negligible rotation support (Gurvich et al. 2023 ). Stars formed during this phase may naturally produce a population of centrally concentrated old stars that are isotropic and spheroidal in nature, with qualitative resemblance to stars that may contribute to classical bulges and/or inner stellar haloes today.

In Section 2 , we provide an overview of our simulations and the kinematic definitions (thin, thick, and spheroid) we adopt in our analysis. In Section 3 , we present results on the dynamical evolution Table 1. The five simulations we employ for the bulk of this work are summarized in the top section of this table. The second set of seven are used only in Section 3.5 . We list the following: the name of the zoom-in target halo, the stellar mass ( M ) within the central 20 kpc of the halo at z = 0, the radius ( R 90 ) enclosing 90 per cent of M , the halo virial mass using the Bryan & Norman ( 1998 ) definition ( M halo ), the halo virial radius ( R halo ), the resolution of each simulation quantified by the initial baryonic particle mass ( m i ), and the reference that first introduced each halo at the quoted targeted resolution. The remaining columns present derived quantities: the lookback time to the end of the bursty phase/onset of the steady phase ( t B ), the mass-weighted thin-disc fraction ( f thin disc m ), and the luminosity-weighted thin-disc fraction ( f thin disc l ). Hosts with names starting with 'm12' are isolated configurations selected from the Latte suite, whilst the rest are in LG-like pairs from the ELVIS on FIRE suite. The four galaxies marked with an asterisk correspond to minor mergers taking place after the onset of the steady phase. The haloes in each list are ordered by decreasing t B .

weighted) ( L weighted) Romeo 7.4 × 10 10 13 .3 1.0 × 10 12 317 3500 6 .52 0.45 0.70 A m12b * 8.1 × 10 10 9 .8 1.1 × 10 12 335 7070 6 .32 0.37 0.64 A m12i 6.1 × 10 10 12 .8 9.2 × 10 11 318 7070 5 .11 0.33 0.64 C m12f * 8.6 × 10 10 11 .0 1.3 × 10 12 357 7070 4 .36 0.33 0.62 B Juliet 4.2 × 10 10 9 .6 8.5 × 10 11 302 3500 4 .40 0.30 0.62 A Remus 5.1 × 10 10 12 .3 9.7 × 10 11 320 4000 5 .88 0.36 0.62 D Louise 2.9 × 10 10 12 .0 8.5 × 10 11 310 4000 5 .56 0.32 0.65 A Romulus 1.0 × 10 11 14 .2 1.5 × 10 12 375 4000 4 .90 0.37 0.69 D m12m 1.1 × 10 11 11 .3 1.2 × 10 12 342 7070 3 .81 0.34 0.58 E m12c * 6.0 × 10 10 9 .7 1.1 × 10 12 328 7070 3 .70 0.32 0.62 A Thelma * 7.9 × 10 10 12 .4 1.1 × 10 12 332 4000 2 .57 0.27 0.57 A m12w 5.8 × 10 10 8 .7 8.3 × 10 11 301 7070 0 .0 0.24 0.43 F Note. The references are: (A) Garrison-Kimmel et al. ( 2019a ), (B) Garrison-Kimmel et al. ( 2017 ), (C) Wetzel et al. ( 2016 ), (D) Garrison-Kimmel et al. ( 2019b ), (E) Hopkins et al. ( 2018 ), and (F) Samuel et al. ( 2020 ).

of galaxy populations with time. Section 4 is reserved for discussion and conclusions.

Our analysis utilizes cosmological zoom-in simulations performed with the multimethod gravity plus hydrodynamics code GIZMO (Hopkins 2015 ) from the Feedback In Realistic Environments (FIRE) project. 1 We rely on the FIRE-2 feedback implementation (Hopkins et al. 2018 ) and the mesh-free Lagrangian Godunov (MFM) method. The MFM approach pro vides adaptiv e spatial resolution and maintains conservation of mass, energy, and momentum. FIRE-2 includes radiative heating and cooling for gas across a temperature range of 10 -10 10 K. Heating sources include an ionizing background (Faucher-Gigu ère et al. 2009 ), stellar feedback from OB stars, AGB mass-loss, type Ia and type II supernovae, photoelectric heating, and radiation pressure, with inputs taken directly from stellar evolution models. The simulations self-consistently generate and track 11 elemental abundances (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe), and include sub-grid diffusion of these elements in gas via turbulence (Hopkins 2016 ;Su et al. 2017 ;Escala et al. 2018 ). Star formation occurs in gas that is locally self-gravitating, sufficiently dense ( > 1000 cm -3 ), Jeans unstable and molecular (following Krumholz & Gnedin 2011 ). Locally, star formation efficiency is set to 100 per cent per free-fall time, i.e. SFR particle = m particle • f mol / t ff with gas particles stochastically converted to stars at this rate (Katz, Weinberg & Hernquist 1996 ). Note that this does not imply that the global efficiency of star formation is 100 per cent within a giantmolecular cloud (or across larger scales). Self-regulated feedback limits star formation to ∼1 per cent-10 per cent per free-fall time (Faucher-Gigu ère et al. 2013 ;Hopkins 2017 ;Orr et al. 2018 ).

1 ht tps://fire.nort hwestern.edu/

Most of this paper relies on a detailed analysis of five MWmass galaxies, which are summarized in the top section of Table 1 . These zoom simulations are initialized following O ˜ norbe et al. ( 2014 ). Three of these galaxies (with names following the convention m12 •) are isolated and part of the Latte suite (Wetzel et al. 2016 ;Garrison-Kimmel et al. 2017, 2019a ;Hopkins 2017 ). Two, with names associated with the famous duo ( Romeo and Juliet ), are part of the ELVIS on FIRE project (Garrison-Kimmel et al. 2019a , b ) and are set in Local-Group-like configurations, as in the ELVIS suite (Garrison-Kimmel et al. 2014 ). This suite includes three simulations in total, containing two MW/M31-mass galaxies each. The main haloes were selected so that they have similar relative separations and velocities as of the MW-M31 pair in the Local Group (LG). Table 1 lists the initial baryonic particle masses for each simulation. Latte gas and star particles have initial masses of 7070 M , whilst ELVIS on FIRE has ≈2 × better mass resolution ( m i = 3500 M ). Gas softening lengths are fully adaptive down to 0.5 -1 pc. Star particle softening lengths are 4 pc physical and the dark matter force softening is 40 pc physical. The last set of seven galaxies is used only in Section 3.5 to demonstrate sample-wide trends.

This analysis focuses on in situ stars that were born within 10 per cent of the virial radiusfoot_2 of the most massive progenitor of each galaxy o v er time. As in Yu et al. ( 2021 , hereafter Y21 ), we classify star particles using their orbital circularity, = j z / j c ( E ), defined as the ratio of each particle's angular momentum in the ˆ z direction to that of a circular orbit with the same energy (e.g. Abadi et al. 2003 ). The angular momentum direction ˆ z is set by total stellar angular momentum within 10 kpc of each galaxy's center at the lookback The cyan block shows our definition of thick-disc stars, which we set to be those with 0.8 > ≥ 0.2. The yellow block marks spheroid stars, which we define as those with < 0.2. The percentage of stars in each block is shown in the legend. We see that these definitions produce components that qualitatively resemble geometrically defined discs and spheroids. Note that the inner ∼1 kpc of the 'thick-disc' component does contain some stars that would likely be identified as a bar or perhaps pseudo-bulge material in a more detailed study. time of interest.foot_3 A star with = 1 is on a circular orbit in the plane of the disc; < 0 implies counter-rotation. We categorize star particles with = 0.8 -1 as thin-disc stars, those with = 0.2 -0.8 as thick-disc stars, and those with = -1.0 -0.2 as spheroid stars. As we show below, we find that circularity appears to be a useful parameter for keeping track of the dynamical evolution of stars o v er time.

The images shown in Fig. 1 illustrate how these classifications manifest morphologically. While we have chosen to classify these components using familiar names, it is important to emphasize that those names are traditionally assigned to populations using morphological decomposition rather than dynamical assignment. Never the less, as can be seen in the left-most image, our thin-disc component does indeed resemble thin discs as usually conceived (see Yu et al. 2021 , where we discuss scale heights, and so on). Our thick-disc component consists of a clear vertically extended disc-like structure, again in keeping with traditional expectations for that name. Ho we ver, it also contains a central, bright, and mildly rotating distribution of stars that qualitatively resembles the observed characteristics of pseudo-bulges (Kormendy & Kennicutt 2004 ). Our spheroid stars (third from left image) have an isotropic configuration that bare qualitative resemblance to a classical bulge and/or inner stellar halo component (e.g. Gao et al. 2020 , and references therein). While we make no direct comparisons to observations in what follows, it is worth keeping in mind that our spheroid component is most aptly associated with an isotropic stellar component that could contain some stars that would be associated with a classical bulge and/or inner stellar halo.

While the orbital circularity works decently as a parameter for classification, we also examine two additional parameters in order to understand how different aspects of the kinematics change o v er time: the 3D orbital circularity, 3D = j / j c ( E ), and alignment angle, θ = arccos( j z / j ). The 3D circularity 3D is defined as the ratio of each particle's total angular momentum to that of a circular orbit with the same energy. Since a circular orbit has maximal angular momentum for a given energy, 3D circularity 3D ranges between 0 and 1, with 3D = 1 corresponding to perfectly circular orbits, and 3D = 0 to purely radial orbits. Alignment angle θ is defined as the angle between each particle's angular momentum to the rotation axis of the galaxy. It describes how aligned the orbit is with respect to the Galactic disc. θ = 0 • corresponds to orbit that lies perfectly in the disc plane, θ = 90 • means that star has a orbit perpendicular to the disc plane, and θ = 180 • indicates that it is counter rotating.

We measure the approximate 'birth' circularity and 3D circularity of each star particle using a post-process analysis of snapshots saved from the simulation. Specifically, we define birth and 3D birth at the first snapshot available after each star particle is formed. The time spacing between snapshots ranges from 16 to 25 Myr, which is small compared to the time-scales of interest ( > 100 Myr). Note that the angular momentum direction ˆ z of the galaxy (which affects birth but not 3D birth ) is set by the total stellar angular momentum within 0.1 R vir at the time of the snapshot immediately following the star's birth (not z = 0). This method has been verified to produce a steadily evolving reference frame that changes from snapshot-to-snapshot typically on the order of a degree or less in orientation (Gurvich et al. 2023 ).

We also find it useful to define a transition time between an early bursty phase of star formation and a later steady phase of star formation. We define the bursty phase to end at a lookback time t B when the standard deviation in 'instantaneous' star-formation rate (SFR) first falls below B = 0.2 times the time-averaged SFR:

Here, the 'instantaneous' SFR is defined to be the rate measured o v er 10 Myr intervals and the time-average rate is measured o v er a 500 Myr interval, as in Y21 . We use this definition to assign a specific bursty-phase time-scale to each galaxy's star formation history. Gurvich et al. ( 2023 ) uses the running scatter in the SFH in 300 Myr windows, σ 300 Myr (log(SFR)), to quantify the fluctuations and define t B as the time after which this quantity remains below 0.3 dex. The difference between different definitions of burstiness results in variations in t B of a few hundreds of Myr, which is relatively small compared to the cosmological time-scales of the transition (see the t B column in Table 1 ). Finally, we quantify inner CGM virialization using the ratio of the cooling time of shocked gas t (s) cool to the free-fall time t ff at an inner radius r = 0.1 R vir . This parameter was introduced by Stern et al. ( 2021), who used it to show that the bursty to steady transition in galaxy star formation coincides with virialization of the inner CGM.foot_4 When t (s) cool /t ff 1, the inner CGM is smooth and largely supported by thermal pressure. In contrast, when t (s) cool /t ff 1, the inner CGM has large pressure fluctuations and is highly dynamic. The bursty-to-steady transition, as well as a transition from thickdisc to thin-disc formation, coincides with the time when the ratio first crosses t (s) cool /t ff 2 (Stern et al. 2021 ;Yu et al. 2021 ).

As mentioned in the introduction, FIRE-2 galaxies display coevolution in a number of parameters related to star formation activity, stellar kinematics, and inner CGM properties. Fig. 2 provides an illustrativ e e xample of this co-evolution in m12i . Similar figures for the four other primary simulations are provided in Fig. A1 . The top panel in Fig. 2 shows the star formation history of m12i as a function of lookback time. The SFR displayed is averaged over both a short time-scale of 10 Myr (SFR 10 , blue) and a longer time-scale of 500 Myr (SFR 500 , red). The relati ve v ariance is much larger at early times than at late times. This is consistent with previous examinations (e.g. Muratov et al. 2015 ;Sparre et al. 2017 ;Flores Vel ázquez et al. 2021 ;Stern et al. 2021 ;Gurvich et al. 2023 ;Hafen et al. 2022 ) that have shown that star formation in massive FIRE galaxies tends to transition from bursty to steady as we approach the present day. The vertical dashed marks the bursty-phase lookback time as defined in the previous section (equation ( 1)), and as described in Y21 . We refer to times later than this transition as the 'steady phase' and times prior to this transition as the 'bursty phase'

The second panel in Fig. 2 presents the 3D circularity, 3D = j / j c ( E ), of newly formed stars (ages < 100 Myr) as a function of lookback time. The black line marks the median value while the shaded region shows the 16th-84th percentile range. Note that the circularity distribution of young stars has a sharp transition once the steady phase begins, at a lookback time of t B 5 Gyr. During the steady phase, stars are born on quite circular orbits, very close to 3 D = 1. About 2.5 billion years before the steady phase begins, we see a gradual 'spin-up' phase, where the orbits of young stars become more circular. Prior to a lookback time of about 7.5 billion years, young stars are born on fairly elliptical orbits. We see qualitatively similar behaviors for all the galaxies in our sample. Note that at early times, the 3D circularities of young stars have typical values near ∼0.7. This is typical of median j / j c ( E ) values one often finds for isotropic orbits in spherically symmetric (non-rotating) systems (e.g. van den Bosch et al. 1999 ).

The third panel in Fig. 2 presents the evolution of the ratio t (s) cool /t ff , which tracks the propensity of the inner CGM to be virialized. The ratio t (s) cool /t ff = 2 is marked. At early times, t (s) cool /t ff 1, the inner CGM is clumpy and dominated by the supersonic infall of cold gas. At late times, t (s) cool /t ff 2, and the GCM becomes hot, smooth, and largely supported by thermal pressure. This transition also coincides with the time that the star formation transitions from bursty to steady 2021), who used it to show that the bursty to steady transition in galaxy star formation coincides with virialization of the inner CGM. The thin dashed line is a smoothed fit to the thick line measured in the simulations. When t (s) cool /t ff 2 (horizontal dashed line) the inner CGM is smooth and largely supported by thermal pressure. In contrast, when t (s) cool /t ff 2, the inner CGM has large pressure fluctuations and is highly dynamic. Bottom row: age distribution of stars that have orbital circularities classified thin disc (magenta), thick disc (cyan), and spheroid (yellow). The solid lines show the distribution for stars classified by their circularity at birth, birth = j z / j c ( E ). The dashed lines show the distribution of stars classified by circularity measured at z = 0. The vertical dashed red line indicates the start of the steady phase in star formation. After this time, stars form on very circular orbits in a thin-disc configuration, and the inner CGM is virialized. The dashed cyan and magenta lines show that some of the stars that form in the steady phase are heated enough to be classified at thick disc at z = 0. Ho we ver, heating appears to be a secondary effect because the age distributions (dashed) show a similar peak as the formation time distributions (solid). Only the ∼ 10 per cent youngest thick disc stars, were formed 'thin' in the steady phase. and the time when stars begin to form on very circular orbits. The circularity of young stellar orbits is enabled by the ability of accreting gas to become coherently aligned in angular momentum space prior to deposition into the galaxy only after the inner CGM becomes smooth and hot (Hafen et al. 2022 ).

The bottom panel in Fig. 2 presents the age distributions of different components identified using circularity = j z / j c ( E ). The solid line shows the distribution for the stars classified using birth circularity birth . We find that almost all the thick-disc and spheroid stars form in the earliest periods of galaxy assembly, whilst thin-disc stars form later after the star formation settles down. The difference between thick-disc and spheroid stars is subtle; o v erall, our identified bulge stars are a bit older than the thick-disc stars. The distribution of thick-disc stars peaks near the time of transition, while bulge star ages peak ∼2.5 Gyr prior.

The dashed lines show the distribution of different populations classified using the circularity measured at z = 0, similar to the method adopted in Y21 . The two bulge distributions are almost the same. For thick-disc stars, the truncation of the age distribution is more abrupt when using birth , while there is an extended tail towards younger age for population classified by . This is likely due to secular disc heating effects that allow stars born with thin-disc like orbits ( birth > 0.8) into the thick-disc regime (with = 0.2 -0.8). The o v erall time sequence, from bulge formation to thick-disc formation to thin-disc formation, stays the same. While there is some degree of disc heating, as we quantify in Section 3.3 , the effect is relatively small compared to the birth-orbit trend.

The main tak eaw ay from this subsection is that levels of starformation burstiness, inner CGM virialization, and birth circularities of new stars appear to be coupled across time. Fig. A1 shows a similar result for the other four galaxies in our main sample. At early times, stars are formed in irregular/more radial orbits; the SFR is quite bursty, and the inner CGM is clumpy, cool, and not virialized.

The red regions in the top panels of Fig. 3 show the time evolution of the 3D orbital circularity, 3D , of young stars (age < 100 Myr) as a function of lookback time for Romeo (left) and Juliet (right). The solid red lines are median values at fixed formation time and the shaded regions plot the 16th-84th percentile range. The black solid lines and shaded regions show for the same quantities for the same stars at z = 0. The vertical red dotted line marks the bursty phase lookback time t B in each galaxy. At lookback times smaller than this, the star formation is steady. The images in the middle panels show edge-on and face-on images 5 of the young ( < 100 Myr) stellar populations at five specific times in the past in each galaxy: 2.7, 4. 7, 6.3, 8.4, and 10.3 Gyr. The images in the bottom panels show the distributions of those same stars today.

For Romeo , we see that 10.3 Gyr ago, the young stars had an irregular morphology (middle panel, far right) with a range of 3D circularities that, in the median, are fairly radial ( ∼0.6). Those stars today are arranged in an isotropic, bulge-like configuration (bottom panel, far right). Just ∼2 Gyr later, the young stars have begun to show some coherent rotation and take the form of a clumpy/irregular disc. Those stars today are in the form of a smooth, thick disc, with median 3D circularity ∼0.8 and significant scatter. At 6.3, 4.7, and 2.7 Gyr (after the steady phase has commenced), Romeo 's young stars are situated in thin discs (middle, left three panels) and and have 3D circularities tightly peaked at 3D ∼ 1.0. These stars remain in relatively thin configurations at z = 0, though there has been some heating/elongation of their orbits o v er time (the grey-shaded band is thicker than the red). We explore this evolution later in Section 3.3 .

The right-hand panels in Fig. 3 show similar behavior for Juliet , but shifted later in time. This galaxy ended its bursty phase more recently ( t B = 4.40 Gyr) than Romeo ( t B = 6.52 Gyr). This results in the later emergence of thin-disc formation. While Romeo had a pronounced thin-disc component 6.3 Gyr ago, Juliet has no thin disc at that time. Only in the most recent image (2.7 Gyr) does Juliet start to have young stars forming in a thin-disc like configuration.

The two sample galaxies we show are representative of our larger simulated sample. Specifically we find that there are three stages in the evolution of the FIRE-2 MW-mass systems o v er time: (1) a v ery early, chaotic bursty phase; (2) a later, quasi-stable, bursty-disc 'spinup' phase; and (3) a stable, thin-disc, 'cool-down' phase. The three stages are marked in the bottom of Fig. 3 . Stage 1 is associated with very bursty star formation. Young stars resemble irregular, chaotic systems in Stage 1. Those stars evolve into an isotropic, classicalbulge-like configuration at z = 0. Later, the host galaxy enters Stage 2, a 'spin-up' phase, where a clumpy and disordered disc begins to emerge. Stars born at this stage end up in a thick-disc configuration at z = 0. As star formation settles down, we enter Stage 3, when thindisc formation occurs. During this time, young stars are forming on extremely circular orbits ( 3D 1.0). Their morphology today also looks much thinner than stars born in Stage 2. Because of the longer steady star-formation phase in Romeo , it also experiences a much shorter stages 1 and 2, resulting in an older bulge, an older thick disc, and also a higher thin-disc fraction ( Y21 ). Juliet , instead, stays in Stages 1 and 2 for a longer time. This results in a smaller thin-disc fraction and a younger thick disc. 5 2D luminosity-weighted in Sloan r band. 3 , the circularity shown here ( j z / j c ) encodes information on the shape of stellar orbits and their orientation with respect to the evolving galaxy's total angular momentum. The vertical red-dashed line marks the bursty phase lookback time t B . The solid red line and shaded area represent the median and onesigma distributions of at birth, while the the black line and the grey areas are the same quantities for the same stars at z = 0. For reference, the horizontal grey-dashed line marks the threshold abo v e which stars are classified as thindisc stars. At early times, stars tend to be born with spheroidal-type orbits ( < 0.2) and to evolve into orbits that are similarly radial with respect to the z = 0 disc ( ∼ 0 in the median). At late times, after the bursty phase, stars form on highly aligned, thin-disc like orbits (red). Those thin-disc stars do get heated o v er time (gre y), but in the median, stars formed after the bursty phase remain in a thin disc. A transition between these two extremes-a 'spin-up' or thick-disc phase-occurs between ∼5 and 8 Gyr ago. During this time, young stars begin to show an increasing level of aligned/coherent rotation, though without the very tight thin-disc alignment we see after the bursty phase has ended.

could be due to both the chaotic, bursty nature of star formation and the swift change in the galaxy's orientation (e.g. Dekel et al. 2020 ;Santiste v an et al. 2021 ). Their current circularity distribution, after a long time of interaction and evolution, becomes centered around ∼ 0, similar to the distribution of an isotropic, spheroidal component. This corresponds to the early chaotic bursty phase we have discussed in Section 3.2 as 'Stage 1'. Approximately 2 Gyr before star formation settles down, the coherence in spin starts to build up and the median value of rises to ∼0.4, with some of the stars already surpassing our threshold for thin-disc kinematics. During this 'Stage 2' phase, the majority stars are born on orbits that we classify as thick-disc. After the star formation becomes very steady, all the young stars are formed in an extremely coherent manner with ∼ 1. The distribution of their current circularity (grey) is wider, especially for stars that formed at earlier times, and becomes more narrow at late times, indicative of a fairly steady heating rate. This 'thickening' appears to be a result of a combination of both dynamical heating and vertical torquing, which we explore in more detail below.

Fig. 5 illustrates how stellar v alues e volve in m12i in a slightly dif ferent way, no w with emphasis on our classification categories. The left-hand panel is similar to the most right-hand panel of Fig. 1 , but now we show the distribution of birth circularity birth for all stars in the galaxy, instead of the circularity today. Using our standard definition, stars are classified as having thin-disc orbits at birth (magenta), thick-disc orbits at birth (cyan), and spheroid-like orbits at birth (yellow). The mass-weighted fractions of the different components are listed on the plot. Note 43 per cent of stars are born on thin disc orbits. This is higher than the fraction of stars that are on thin-disc orbits today (33 per cent, see Fig. 1 ). This is consistent with the observation that some stars that are born in thin disc orbits become heated o v er time.

Figure 5. How do the orbital classifications of stars evolve from birth to the present day? Left-hand panel: the birth circularity birth distribution for stars that currently reside within R 90 (12.8 kpc) in m12i . The magenta block marks stars that formed with thin-disc kinematics, with birth ≥ 0.8. The cyan block marks thick-disc stars, which we define to be those with 0.8 > birth ≥ 0.2. The yellow block marks spheroid stars with birth < 0.2. The mass-weighted fraction of stars in each block is shown in the legend. Right-hand panel: the present-day circularity distributions for each of the three populations defined by their birth circularities in the left-hand panel. The colours in each of the right-hand panel match the components they track as defined on the left-hand panel: thin-disc, thick-disc, and spheroid, from top to bottom. Note that the vertical (mass-density) axis for each panel is different. Top right: the current circularity distribution for stars born with thin-disc like orbits (magenta). Most of these stars retain thin-disc orbits. The distribution does have a tail to wards lo wer circularity, which is due to disc heating, but 65.8 per cent of the stars still have current circularity greater than 0.8. Middle right: the distribution of the stars identified as thick disc (cyan block) in the left-hand panel evolve to have a distribution that looks like a combination of a thick disc and isotropic bulge population. Bottom right: stars born with spheroid-like orbits (yello w) e volve into a nearly isotropic distribution at the present day, though there is a slight preference for prograde orbits. The red-shaded region shows what happens to an isotropic spheroid population that is forced to be symmetric around birth = 0 (dashed grey, lower left). We see that this choice has a minimal impact on the shape of the distribution today.

The right-hand panels of Fig. 5 display the distribution of current circularities for each of the birth-based categories. Stars that are born with thin-disc circularities ( birth > 0.8, magenta) end up in clearly prograde, disc-like configuration today, with an distribution that peaks at 0.9 and with 66 per cent of the distribution still in a thin-disc (with > 0.8). There is a tail less well-aligned orbits, likely as a result of heating or torquing o v er time (see below). For stars born with thick-disc orbits (cyan), their current circularity distribution resembles a combination of a prograde/thickdisc component (peaking at ∼ 0.6) and a broad isotropic/bulge component centered on ∼ 0. Stars born with spheroidal orbits (yellow) tend to stay in an isotropic distribution today, with a mild prograde asymmetry.

Ideally, an isotropic spheroid population would have symmetric circularity distribution around 0. For the simplicity of this work, we have adopted a sharp cut in circularity to define spheroids ( < 0.2).

We hav e e xplored the impact of imposing symmetric distribution, mirroring the observed < 0 distribution (hashed yellow in the lefthand panel) for > 0 (marked by grey-dashed line in the left-hand panel). In this dashed-line selection, we have simply chosen particles randomly at fixed such that they inhabit the symmetric distribution sho wn. The dif ference between including the stars with = 0.0 -0.2 at birth is minimal since statistically only a small fraction of stars fall in this region (yellow block abo v e the gre y-dashed line). Ov erall, we observed a trend that stars born with bulge-like orbits retain spheroid-like orbits today (dashed grey distribution on the right).

Fig. 6 shows the evolution of circularity for different mono-age stellar populations in Romeo , m12b , m12i , and m12f , from top to bottom, respectively. Here we have binned star particles based on their birth ages in 100 Myr increments and have calculated the median value of circularity in each age bin as a function of lookback time. Lines are coloured by the stellar age of the population (see colour bars) and are dashed when the host galaxy is in its bursty phase. They are solid during the steady phase for each galaxy.

The left-hand column of Fig. 6 tracks the evolution of each monoage population o v er all time. The colour bars are scaled differently for each galaxy in order to emphasize the regularity when normalized with respect to the bursty phase lookback time. Specifically, colour bars are set so that the white colour corresponds to the transition time from bursty to steady star formation. Populations that formed Figure 6. Evolution of the population-median circularity = j z / j c ( E ) of mono-age stellar populations for Romeo , m12b , m12i , and m12f . Each curve is coloured by stellar age of the population. Left-hand panel: these figures span all of cosmic time. The colour bars are slightly different for each galaxy, such that the bursty-phase transition time is al w ays white. Red lines track the evolution of mono-age populations that formed during the bursty phase. Blue lines track the evolution of mono-age populations that formed during the steady phase. The older populations (red lines) tend to evolve to spheroid-like orbits today ( ࣠ 0.2) while younger populations (blue lines) form with much higher circularities ( ∼ 1). Right-hand panel: zoomed-in plots for only the populations that have formed during the steady phase for each host galaxy. Note that the time-axis is slightly different for each galaxy because each galaxy has a different bursty-phase lookback time. The stellar ages are again mapped to the colour bar. The key takeaway is that, in the steady phase, new stars are born with thin-disc orbits and remain abo v e the 0.8 thin-disc threshold in most cases. That is, most stars that form during the steady-star-formation period have thin-disc kinematics today. Note that, for Romeo and m12i , all populations form during the steady phase experience an initial drop in j z / j c but quickly plateau to nearly nearly constant values. Both m12b and m12f undergo ∼1/10 merger events during the steady phase and this appears to drive slightly different behaviour, with less clear plateauing. The time of each merger is indicated by an arrow, and there is an obvious dip at the time of merger. See the text for a discussion of these dips.

during the bursty phase are coloured red and populations formed during the steady phase are coloured blue. The two grey-dashed lines mark our thresholds for the classification of three components (thin disc, thick disc, and spheroid). We observe the same behavior in all four galaxies: stars that form during the early bursty phase (red) mostly have orbits characteristic of spheroid stars ( < 0.2) while populations born during steady phase (blue lines) have median circularities characteristic of thin-disc stars. While the red (earlyforming) populations show large fluctuations, they all approach ∼ 0 (characteristic of an isotropic population) at late times. Interestingly, there is a brief period in m12i at a lookback time of ∼11 Gyr where the young stars have more coherent, thick-disc like orbits ∼ 0.4; though as in all other cases they evolve to have isotropic orbits with respect to the final system at z = 0.

The right-hand panels of Fig. 6 are 'zoomed-in' versions of the left-hand panels and concentrate e xclusiv ely on the steady phase. The lines are again coloured by the stellar age. Stars formed during this period are born along a tight plane with extremely circular orbits ( ∼ 1). Their circularities decrease slowly with time, but in almost all cases retain orbits characteristic of thin discs to the present day, with 0.8 at z = 0. The steady-phase mono-age populations for Romeo , m12i , and m12f all show more rapid evolution in orbital circularity within the first ∼500 -1000 Myr after being born, and begin to flatten somewhat at late times. This behavior is most extreme for Romeo where the evolution plateaus to near constant circularity at recent lookback times.

The tracks in m12b show somewhat different behavior, with all populations evolving slowly but continuously to lower circularity until the present day. This is likely because the disc experiences a merger with an LMC-size satellite at a lookback time of approximately 2.9 Gyr (marked by the arrow). This merger, which comes in on a polar orbit, produces a sharp feature in all of the mono-age tracks (and also drives a small starburst; see Y21 ). This feature is almost certainly an an artefact of our definition of spin axis. When the merger happens, the satellite galaxy gets close to the host galaxy. This affects the ˆ z direction of the galaxy, since we use all the stars within 0.1 R vir in calculating the direction of total angular momentum. While the mass of stars in the accreted object is small, they are all moving coherently at the time of the merger and this changes the bulk angular momentum direction dramatically. After the merger, the accreted stars become phasemixed and no longer contribute substantially to the total angular momentum.

Even with the fairly significant merger in m12b , the median circularity of the oldest post-bursty-phase population remains within our thin-disc classification ( > 0.8). Note that m12f also experiences a similarly sized merger, at a lookback time of ∼300 Myr (see arro w). The ef fect on e volution is less dramatic, but still visible: the 'plateauing' in circularity stops after the merger, and we see a mild drop afterwards. The weaker effect seen with this merger is possibly because the orbit is prograde (though interestingly, this merger drives a larger starburst, see Y21 ).

Figs 7 and 8 complement Fig. 6 by tracking the evolution of 3D circularity 3D and the alignment angle θ of the orbits with respect the existing disc. The 3D circularity 3D quantifies how circular the orbits are regardless of orientation with respect to the disc plane. The alignment angle θ quantifies the level of alignment of angular momentum in the stars with respect to the rotation axis of the galaxy disc. While provides a measure of thickness because it depends on orbital orientation and elongation, these figures provide insight on how the change in circularity is driven separately by orbital elong ation (reg ardless of orientation of the orbit) or vertical torquing of the orbit out of the disc plane.

In the left-hand panels of Fig. 7 , we see that stellar populations born during the bursty phase typically start off with fairly radial orbits ( 3D 0.6). The left-hand panels of Fig. 8 show that these early-forming stars also have significant misalignment, with median θ ∼ 90 • in most cases. Note that this does not mean that most stars are actually orbiting on orbits perpendicular to the disc. Rather, the distribution of orbital planes is quite random, such that in the median the angle is θ ∼ 90 • .

As mentioned abo v e, m12i deviates some from the trends we see in the other cases. In particualr, there is a brief phase at ∼11 Gyr lookback time where populations of young stars have 3D circularities as large as ∼0.8. At the same time, the median orbital alignments are slightly tighter, with values sometimes as small as ∼45 • . However, the angular momentum direction of the system varies quickly enough that we also see median orbits counter-rotating with θ ∼ 135 • o v er short time interv als. This le vel of irregularity in total angular momentum direction at early times helps explain why most of the ∼11 Gyr-old stars in m12i have evolved to inhabit an isotropic 'spheroid' population at late times. This can be seen by comparing the red (at birth) and grey (now) distributions o v er the same time period in Fig. 4 .

The right-hand panels of Figs 7 and 8 focus on the evolution of mono-age populations since the bursty phase has ended in each galaxy. There are a few interesting trends visible here. First, in all galaxies except m12b , we see that stellar populations stars born with gradually increasing degree of circularity and alignment within the plane as we approach the present day. That is, even the youngest stars become increasingly thin-disc like as time progresses. This is not the case in m12b , where the merger (marked by the arrow) resets the trend: post-merger, the orbits steadily become more aligned with the plane and have more circular orbits.

After birth, the stars' orbits become gradually less circular and more misaligned, though the level heating/torquing is mild and in most cases the mono-age lines plateau to near-constant values after a short ∼500 Myr period just after the stars form, where the elongation/torquing behavior is quickest. Interestingly, a late-time prograde merger in m12f (see arrow) causes a fairly significant drop in 3D circularity but not in the alignment angle.

The preceding discussion has shown that while there is some secular or merger-driven heating of orbits after formation, the degree of heating is small compared the broader trends with birth orbits o v er cosmic time. The scale of the v ertical ax es in the right-hand column of Figs 6 -8 is quite narrow compared to those of the left, and in most cases even the first populations to form after the bursty phase ends do not evolve enough to become thick-disc stars by our classification.

Fig. 9 explores the question of heating rate in more detail by plotting the average rate of change in orbital properties of all mono-age stellar populations between the time of their birth and today. The left-hand panels show the average change in circularity per Gyr, calculated by dividing nowbirth by the stellar age. Red dots represent populations that form during bursty phase and blue dots show the ones born in the steady phase.

The grey-shaded region provides a sense of the rate 'thickening' that would be required to change a star born in a perfectly thin disc ( = 1.0) to one that has just enough orbital circularity to inhabit the thick disc at z = 0 by our definition ( = 0.8). The grey band rises towards the current epoch because stars born later have less time to be heated and thus need a higher heating rate to join the thick disc before z = 0. Note that this comparison only makes sense for stars born in steady phase (blue) since stars born during the bursty phase are typically born with thick-disc or spheroid-like orbits already. We see that in almost all cases, the blue points sit below the gre y re gion. This means that the av erage 'thickening' since birth is simply not enough to turn thin-disc stars into thick-disc stars in most cases.

The middle and right-hand columns explore the rate of change in 3D circularity 3D and alignment angle θ . We show these results in order to disentangle the effects of orbital elongation and torquing. As we can see, both mechanisms work to thicken the newly born stars, but the amount of heating and torquing is minimal. The change in 3D is smaller than 0.05 per Gyr and the change in θ is smaller than 5 • per Gyr.

By comparing the magnitude of the left and middle columns, we see the thickening rate (left-hand column) very similar in magnitude to the heating rate (middle column), suggesting that heating (rather than torquing) is dominating the thickening evolution.

Note that the average rate of change typically shows an upward trend in stars with the youngest age ( < 1 Gyr), which is driven by the fact that the denominator in the rate calculation is al w ays the time since birth. We find that the heating rates for stellar populations are usually highest within the first ∼500 Myr of formation. This makes sense in our simulations because star Figure 7. Evolution of the median 3D circularity 3D = j/j c (E) for mono-age stellar populations in Romeo , m12b , m12i , and m12f . As in Fig. 6 , each curve is coloured by stellar age of the population. Left-hand panel: each line tracks the evolution of 3D o v er cosmic time, from birth (right) to the present day (far left). For each galaxy, red lines correspond to populations that formed during the bursty star-formation phase, while blue lines track stars that formed during the steady phase. The older populations (red lines) tend to form with less circular orbits ( 3D 0.6 -0.8) while younger populations (blue lines) form on quite circular orbits ( 3D ∼ 1). Bursty-phase stars tend to remain on fairly radial orbits today. Steady-phase stars evolve to slightly less circular orbits today but remain more circular than an of the bursty-phase populations. Right-hand panel: zoomed-in plots for only the populations that have formed during the steady phase for each host galaxy. These stars are born on very circular orbits, and become slightly less circular o v er time, with the sharpest evolution seen within the first ∼500 Myr after formation. The arrows in m12b and m12f indicate times that these galaxies experience ∼1/10 mergers. Note that the merger in m12b , which is polar, does not cause a sharp feature in the 3D circularity of orbits; this suggests that the distinct feature we see in Fig. 6 is driven by an artefact of how we define the orientation of the galaxy (see te xt). Conv ersely, the merger in m12f , which is prograde, does cause a sharp feature in 3D circularity, suggesting that this merger has heated existing stars to some extent and driven their orbits to become less circular. formation occurs only in o v erdense, self-gravitating structures, where also the strongest dynamical perturbations after formation will occur. This can be seen in the right-hand columns of Fig. 7 , for e xample. We further e xplore these effects in Appendix B (Figs B1 and B2 ).

The evolutionary trends we have presented in this section show that stars born on circular, disc-like orbits remain mostly in thin discs today and that stars born on elliptical orbits remain on similar orbits today. Broadly speaking, stars that currently exist in kinematic/morphological classes today were born that way. Though we do see some dynamical heating and torquing o v er time, the level of these evolutionary effects are secondary compared to the birth properties. Of course, given our small sample size and the fact that we have focused on MW size galaxies, our results may not be universal for all galaxy classes. In Y21 , we explored the importance of mergers in shaping thin/thick-disc formation in 12 MW size simulations. Of these galaxies, four had significant mergers after star formation settles down, including two that are discussed abo v e, m12b , and m12f . These mergers do not destroy or disrupt the thin disc and contribute only in a second-order way to populating the thick disc. The role of mergers in shaping activity during the bursty phase will be the topic of future work.

Figure 8. Evolution of the median alignment angle θ = arccos ( j z /j ) of mono-age stellar populations with respect to the evolving galaxy's total angular momentum direction in Romeo , m12b , m12i , and m12f . As in Fig. 6 , each curve is coloured by stellar age of the population. Left-hand panel: the colour bars are set so that all the populations form during the bursty phase are coloured red while the populations form during the steady phase are coloured blue. Orbits of older populations (red lines) tend to be in random directions with respect to the orientation of the galaxy while younger populations (blue lines) form more aligned with the disc plane. Right-hand panel: zoomed-in plots for only the populations form during the steady phase for each host galaxy. All populations form fairly aligned with the disc ( θ ࣠ 10 • ). Soon after formation, most populations experience some evolution initial change in θ , but later plateau at a nearly constant v alue. The arro ws in m12b and m12f indicate times that these galaxies experience ∼1/10 mergers. Note that the merger in m12b , which is polar, causes a sharp feature in the relative alignment of stellar orbits; this suggests that the distinct feature we see in Fig. 6 is driven by an artefact of how we define the orientation of the galaxy (see te xt). Conv ersely, the merger in m12f , which is prograde, does not cause a sharp feature in orientation of the orbits, suggesting that this merger has driven orbits to be come more elongated (see Fig. 7 ) but has not affected their orientations significantly.

In this subsection, we briefly explore the age distributions for our spheroid, thick-disc, and thin-disc components in all 12 of the simulated galaxies listed in Table 1 . Our aim is to illustrate how the ages of these components scale systematically with the burstyphase lookback times for each galaxy.

Every line colour in the three panels of Fig. 10 corresponds to a different simulated galaxy. The colour code is mapped to the bursty-phase lookback time of each galaxy. These have values that range from t b = 6.5 Gyr ( Romeo , purple) to t b = 0 ( m12w , yellow). Note that we have included m12w for complete-ness, even though it is still (barely) in its bursty phase by our definition.

The top panel shows the stellar age distribution of thin-disc stars at z = 0, classified by , as illustrated in Fig. 1 . We have smoothed all lines with a Gaussian filter. 6 The middle and bottom panels show thick-disc stars and spheroid stars, respectively . Importantly , all ages are plotted with respect to the bursty-phase lookback time. That is, we show stellar age -t b instead of stellar age, in order to highlight . Blue dots show rates for populations born after the star formation settles down. The grey bands in the left-hand panel estimate the amount of 'thickening' needed for stars born thin ( birth = 1.0) to become thick ( now ≤ 0.8). Note that, this region is only rele v ant as a comparison for population born in steady phase since stars born in bursty phase have now < 0.8 at birth. The middle panels and right-hand panels show the kinematic heating and vertical torquing rates, respectively. The changes in both quantities are relatively small compared to the trends we observe in birth orbital properties with cosmic time. Note that the blue points rise sharply at the youngest stellar ages in m12f . This is because this galaxy experiences a late-time merger. The other galaxies show a much milder rise in heating rate for the youngest galaxies. This is an artefact of our measuring the rate of heating since birth. All populations show their most dramatic heating within 1 Gyr of formation. The youngest stars are still in their rapid heating phase when the time-average heating rate is measured. the correspondence between the transition from bursty to steady star formation and the formation of different components in the galaxy.

When offset by t b , all galaxies and galaxy components show a remarkable similarity in their relative age distributions, especially the spheroid populations. Almost all the spheroid stars formed during the bursty phase (age -t b > 0), while thin-disc stars dominate after star formation has settled down (age -t b < 0). As discussed abo v e, one galaxy in our sample, m12w , never really settles down, so we have defined its busty phase to end at the present day ( t b = 0 Gyr; yellow lines in Fig. 10 ). This likely explains why it is an apparent outlier. This galaxy does have some stars with orbits that fall within our thin-disc category, but morphologically, its disc is much thicker than all of the other runs with a relatively small thin-disc fraction ( Y21 ).

The thick-disc stars show an age distribution intermediate to thindisc and spheroid stars. The distributions all peak prior to the end of the bursty phase but slightly closer to the end of the bursty phase than do the spheroid age distributions. Though the distributions are similar, thick-disc stars have a tendency to be slightly younger than spheroid stars in most galaxies. In addition, while the majority of thick-disc stars form during the bursty phase, there is a tail of slightly younger thick-disc stars that formed after the bursty phase ended. This is not seen in the spheroid population. Based on the analysis presented abo v e, the tail of younger, post-bursty-phase, thick-disc stars is likely populated by stars that were born thin and subsequently heated to have tilted and more elongated orbits.

The abo v e discussion is consistent with the broad picture suggested by the analysis in previous sections: the orbital properties of z = 0 stars in MW-size galaxies track an age sequence linked closely to the transition between bursty and steady star formation in galactic evolution. Isotropic spheroid populations in MW-mass galaxies are formed early, from stars that are born on fairly radial orbits when star formation was irregular and bursty. Thin-disc populations form late, after the bursty phase has ended. Thick-disc stars are an intermediate age population, the bulk of which form prior to the end of the b ursty phase, b ut slightly later than the isotropic spheroid. If these simulations are correct, then spheroid age and thick-disc age should be tightly correlated with the bursty-phase lookback time.

We have used FIRE-2, MW-mass galaxy simulations to track the orbital evolution of stars from birth to the present day. As illustrated in Fig. 1 , we use stellar orbital properties at late times to classify stars with thin-disc , thick-disc , and isotropic spheroid . All of our galaxies have an early, bursty star formation period that transitions to a steady, more uniform SFR at late times (Figs 2 and A1 ). During the earliest period of the bursty phase, galaxies hav e clumpy, irre gular morphologies (Fig. 3 ) and stars are born on fairly radial orbits. These stars evolve to inhabit an isotropic spheroid population today (Figs 3 and 4 ).

At late times, after galaxies transition from bursty to steady star formation, new stars are born on extremely circular orbits along a narrow plane and the galaxy begins to build a substantial thin disc (Fig. 4 ). Prior to the transition from bursty to steady star formation, new stars begin to show more substantial rotation than in the bursty phase. During this 'spin-up' phase, young stars have more circular orbits in the median, though there is a fair amount of scatter and orbital misalignment. An example illustrating this is shown Fig. 4 at a lookback time of ∼7 Gyr. Stars formed during this phase contribute substantially to thick-disc components at the present day.

We hav e e xplicitly use the term 'spin-up' to make connection with the observational findings of Belokurov & Kravtsov ( 2022 ), who used it in connection with a transition epoch they identified in MW stars at metallicity [Fe/H] -1. Below this metallicty, they find that Milk Way stars have median tangential velocity of ∼0 km s -1 (typical of a spheroid). At [Fe/H] -1, stars quickly transition to median tangential velocities of ∼100 km s -1 which his typical of a thick disc.

In summary, and as illustrated in Fig. 3 , our galaxies progress through three stages in cosmic evolution: (1) a very early, bursty-irregular, chaotic phase; (2) a later, bursty-disc, 'spin-up' phase; and (3) a late-time, thin-disc, 'cool-down' phase. Stars born during the first two stages evolve into isotropic spheroid and thick-disc populations today. One corollary to this time sequence is that our thick-disc and spheroid populations are mostly formed at early times prior to the formation of the thin disc. As detailed in Figs 6 -8 , stellar populations have median orbital properties at z = 0 that are quite similar to their median orbital properties at birth.

While birth orbit appears to be the most important driver in predicting where stars end up at z = 0, mergers and secular heating do affect the dynamics of stars at late times in our simulations. Figs 9 , B1 , and B2 explore the evolutionary 'thickening' after birth. Most mono-age populations experience the most orbital heating/torquing right after they form. After ∼0.5 Gyr, the rate of heating drops and the median orbits become fairly stable. The majority of-though not all-stars born with thin-disc orbits remain in thin-disc orbits at the present day. There is some degree of heating that occurs after the steady/thin-disc phase begins, and these heated stars populate the youngest-age distribution tail of the thick disc.

Of course, the scenario we discuss here, especially as it relates to thin or thick-disc orbital evolution, would change in the event of major late-time mergers. Ho we ver, for galaxies that are MW size and smaller, we expect that the major merger rate will be quite rare at late times (Stewart et al. 2009 ;Fakhouri, Ma & Boylan-Kolchin 2010 ;Rodriguez-Gomez et al. 2016 ;Hu ško, Lacey & Baugh 2022 ). While our sample of galaxies is small, there was no specific selection for systems without major mergers. Half of the 12 galaxies in our sample experience mergers with stellar mass ratios greater than 0.15 since z = 3. All of these lar ge mer gers occur before z = 1. None of our 12 galaxies have a merger with mass ratio larger than 0.15 since the onset of the thin-disc phase. These statistics are consistent with the expectations from larger samples. For example, Rodriguez-Gomez et al. ( 2016) used used a cosmological sample of haloes to show that that the major merger rate (involving stellar mass ratios > 0.25) for MW size galaxies at late times is only ∼0.02 Gyr -1 . This suggests that, while rare, significant late-time mergers should occasionally occur at this mass scale. In those rarer cases, especially those involving close to one-to-one mass ratio mergers, we would expect orbits to be randomized and destroy discs.

Two of the five primary simulations we have analysed experience substantial, though still minor (LMC-size). late-time mergers during the steady phase (see the arrows in right-hand panels of Fig. 6 ). These mergers, with stellar mass ratios ∼0.1, neither destroy thin discs nor significantly alter the orbital properties of young stars. Mergers do appear to heat some thin-disc stars enough to populate thick-disc components, but only by populating the tail of youngest thick-disc stars. An analysis like ours performed on a much larger sample of galaxies that include significantly larger later-time mergers will be required to understand how and how often late-time major mergers affect the picture presented here.

As discussed in the introduction, there have been a series of FIRE-2 papers that have examined the relationship between star formation burstiness, galaxy kinematics, galaxy morphology, and the development of a hot gaseous haloes around galaxies (Ma et al. 2017 ;Stern et al. 2021 ;Yu et al. 2021 ;Gurvich et al. 2023 ;Hafen et al. 2022 ). In particular, all galaxies of sufficiently high mass in these simulations experience a fairly sharp transition from bursty to steady star formation. The lookback time to this transition, t B , is different for each galaxy, coincides with the virialization of the inner CGM (Stern et al. 2021 ;Yu et al. 2021 ), and also correlates with distinct changes in the angular momentum properties of both the CGM and ISM (Gurvich et al. 2023 ;Hafen et al. 2022 ).

The lookback time to the bursty-to-steady transition, t B , sets a characteristic age for morphological components in our galaxies at z = 0 (Fig. 10 ). Specifically, virtually all in situ isotropic spheroid stars form prior to t B . Moreo v er, when stellar ages are shifted with respect to t B for each galaxy, their age distributions are remarkably similar. This means, for example, that the youngest stars in the in situ spheroid should have ages that map to the bursty-to-steady transition time, which corresponds to the time of inner CGM virialization. As shown in Fig. 10 and Yu et al. ( 2021 ), the median ages of thick-disc stars are also tightly correlated with t B . The youngest tail of heated thick-disc stars are formed after t B . The majority of thin-disc stars form after t B . These trends provide a potential avenue for linking the observed age distributions of stars in the MW and similar local galaxies to the nature of in situ star formation o v er cosmic time.

Future work that compares predicted distributions of galaxy morphologies and kinematic properties to observed characteristics at intermediate and high redshft (Smit et al. 2018 ;Wisnioski et al. 2019 ;Zhang et al. 2019 ;Ferreira et al. 2022 ;Robertson et al. 2022 ;Wu et al. 2022 ) will provide useful tests of these simulations.

Each panel of Fig. A1 is similar to Fig. 2 in the main text, but shows results for the other four galaxies in our primary sample. The top panels track the SFR as a function of lookback time, with a clear transition from bursty to steady star formation marked by the vertical dashed red line (at t B ). The middle panel shows the median (solid) and one-sigma distribution (grey-shaded band) of 3D circularities for all young stars (age <100 Myr) as a function of lookback time. The third panel shows the evolution of the ratio t (s) cool /t ff , which tracks the propensity of the inner CGM to be virialized. The ratio t (s) cool /t ff = 2 is marked by a horizontal dashed line. During the early, bursty phase, we see that t (s) cool /t ff 1. When this condition is met, the inner CGM is dominated by the supersonic infall of cold, often clumpy, gas. At late times, t (s) cool /t ff 2, and the GCM becomes hot, smooth, and largely supported by thermal pressure.

During the steady phase, stars are born on quite circular orbits, very close to 3 D = 1, with very little scatter. During the early bursty phase, the scatter is much higher. A few billion years before the bursty phase ends, we see a gradual 'spin-up' in the orbits of young stars, as the circularities rise towards unity. Note that m12b and m12f (bottom two panels) Both undergo ∼1/10 merger events during the steady phase, at lookback times of ∼3 Gyr and ∼0.3 Gyr, respectively. At these times, the distribution of young stellar orbits in both galaxies becomes slightly more radial, but the period is fairly brief. The feature is more pronounced in m12b , where the merging orbit is polar. The merger in m12f is prograde.

In Section 3.4 , we discussed the evolution of the orbital properties of mono-age stellar populations with cosmic time, and presented the average rate of change in orbital properties from birth to the present day in Fig. 9 . We specifically explored the change in circularity , 3D circularity 3D , and alignment angle θ , per Gyr.

Here, in Figs B1 and B2 , we present the rate of change measured o v er shorter time-scale to understand how the heating rates evolve o v er time. Instead of averaging over the entire age of a stellar population, we select 500-Myr time windows right after stars form and right before z = 0. Since the time period we choose is relatively small compared to the cosmological time-scale of the transition in the galaxy evolution, it is reasonable to think of these two values as the 'instantaneous' thickening/heating/torquing effect.

Red dots represent stellar populations born during the bursty phase and blue dots show for stars that form after star formation settles down. The grey-dashed lines align with zero, to guide the eye. Dots that lie on this line have no change in median , 3D , or θ during the time period explored.

Generally speaking, the instantaneous rates are lowest near z = 0 and higher just after birth. Stars born during the bursty phase (red dots) have significantly higher rates of orbital evolution just after birth.

Together with the results from Section 3.3 , these results suggests that disc heating/torquing is not significant enough to turn most thindisc stars into thick-disc stars, and thus is unlikely to be the primary formation mechanism for the bulk of the thick-disc stars in the MWmass systems. 9 , here we quantify the amount of thickening/heating/torquing happening for different populations at birth and now by calculating the change in circularity , 3D circularity 3D , and alignment angle θ 500 Myr after alignment angle r birth and 500 Myr before z = 0. For all three quantities, the change happens mostly right after birth while for the period of time before z = 0, the change for all populations, especially the ones born in steady phase, is almost zero. This could also be seen in Fig. 9 as the lines are getting flatter around 0 Gyr.