Exploring the fate of primordial discs in Milky Way-sized galaxies with the GigaEris simulation

ABSTRACT

Recent observations with JWST and Atacama Large Millimeter/submillimeter Array (ALMA) have unveiled galaxies with regular discs at significantly higher redshifts than previously expected. This appears to be in contrast to our understanding of galaxy evolution at high redshift. Additionally, the stellar populations of the Milky Way suggest that the bulk of the Galactic thin disc formed after |$z=1$|⁠, raising questions about the history, evolution, and survivability of primordial discs in Milky Way analogues. Here, we use GigaEris, a state-of-the-art N-body, hydrodynamical, cosmological ‘zoom-in’ simulation with a billion particles within the virial radius, to delve into the formation of the early kinematically cold discs (KCDs), defined by their ratio between the mean rotational velocity and the radial velocity dispersion, of a Milky Way-sized galaxy at redshifts |$z\gtrsim 4$|⁠. Our analysis reveals a primarily inward migration pattern for disc stars formed at |$z \gtrsim 6$|⁠, turning into a mix of inward and outward migration at later times. Stars migrating outwards undergo minimal kinematic heating, and might be identified as part of the thin disc formed at much later epochs. We find that approximately 76 per cent of all stars formed in the KCD at |$z \sim 7$| become part of a pseudo-bulge by |$z = 4.4$|⁠. This proportion decreases to below 10 per cent for KCD stars formed at |$z \lesssim 5$|⁠. The inward migration of stars born in our KCDs at |$z \gtrsim 4$| deviates from the expected inside-out formation scenario of thin discs at lower redshifts. Our results suggest a novel ‘two-phase’ disc formation process, whereby the early discs transform primarily into the pseudo-bulge within less than a billion years, and the present-day thin disc forms subsequently from higher angular momentum material accreted at later times.

1 INTRODUCTION

Since the initial proposal by Gilmore & Reid (1983) of a dual-component structure within the disc of the Milky Way (MW), numerous observational studies have supported the idea of distinct thin- and thick-disc components. Thin-disc stars, located closer to the plane and with low-eccentricity orbits, exhibit low [|$\alpha$|/Fe] enhancement and a metallicity distribution centred around solar values, indicating a similar abundance of elements heavier than helium to that of our Sun, and ranging from supersolar values down to [Fe/H] |$\sim -0.7$| (e.g. Robin et al. 1996; Norris 1999; Ojha 2001; Larsen & Humphreys 2003; Jurić et al. 2008; Bensby, Feltzing & Oey 2014; Recio-Blanco et al. 2014; Hayden et al. 2015). In contrast, the thick disc extends to larger distances from the plane and its stars exhibit higher [|$\alpha$|/Fe] enhancements and a lower metallicity content than those of the thin disc, with a broader metallicity distribution centred around [Fe/H] |$\sim -0.55$| (e.g. Gilmore, Wyse & Kuijken 1989; Katz et al. 2011). Kinematically, thick-disc stars have larger velocity dispersions than the stars in the thin disc (e.g. Gilmore, Wyse & Norris 2002; Soubiran, Bienaymé & Siebert 2003; Parker, Humphreys & Beers 2004; Wyse et al. 2006; Mackereth et al. 2019). Despite these differences, there is considerable overlap between the two components, making their distinctions less clear-cut (see, e.g. Minchev, Chiappini & Martig 2014; Hayden et al. 2015; Queiroz et al. 2020). Therefore, the origin of the thick and thin discs remains the subject of ongoing debate and uncertainty (see, e.g. for cosmological simulations Grand et al. 2018; Mackereth et al. 2018; Debattista et al. 2019; Buck et al. 2020; Agertz et al. 2021; Yu et al. 2023; Iza et al. 2024; and for isolated simulations Quinn, Hernquist & Fullagar 1993; Abadi et al. 2003; Brook et al. 2005; Villalobos & Helmi 2008; Sales et al. 2009; Schönrich & Binney 2009a, b; Aumer, Binney & Schönrich 2016; Clarke et al. 2019).

The prevailing theory amid the ongoing debate is the formation sequence proposed by Bird et al. (2013), wherein the thick disc forms first, followed by the thin disc, in an ‘upside-down’ and ‘inside-out’ manner. In this sequence, disc stars initially form within a thicker and radially compact structure before being born in progressively thinner and larger structures (see, e.g. Bird et al. 2013, 2021; Wang et al. 2019; Buck et al. 2020; Agertz et al. 2021).

New observations by the JWST (Gardner et al. 2006) have provided an unprecedented view into the early Universe, unveiling well-structured disc galaxies at significantly higher redshifts, |$3 < z < 8$|⁠, than previously thought (e.g. Ferreira et al. 2022, 2023; Kartaltepe et al. 2023; Robertson et al. 2023; Tohill et al. 2024). Studies with the Atacama Large Millimeter/submillimeter Array (ALMA; ALMA Partnership et al. 2015) have also found a large fraction of star-forming dynamically cold discs at as early as |$z \sim 6.8$| (e.g. Smit et al. 2018; Rizzo et al. 2020, 2021; Roman-Oliveira, Fraternali & Rizzo 2023). Moreover, in recent years, scenarios have been proposed where the thin disc begins to form concurrently with the thick disc, suggesting the presence of a substantial population of old stars in the former (e.g. Beraldo e Silva et al. 2021; Zhang et al. 2021; van Donkelaar, Agertz & Renaud 2022, and Kretschmer, Dekel & Teyssier 2022 for the formation of old cold gaseous discs). Tamfal et al. (2022) took it a step further and, using the GigaEris cosmological ‘zoom-in’ simulation, showed that thin discs,¹ namely kinematically cold discs (KCDs), are already present at redshift |$z\sim 7$|⁠. This was further supported by the work of Kohandel et al. (2024), who found dynamically cold gaseous discs using the SERRA (Pallottini et al. 2022) suite of zoom-in simulations.

Despite the new results, it is evident that the present-day MW thin disc comprises predominantly young and kinematically cold stars (e.g. Parker et al. 2004; Wyse et al. 2006; Brook et al. 2012). This prompts questions regarding the evolution and fate of the primordial discs detected by ALMA and JWST, as well as those observed in simulations. The results presented by van Donkelaar et al. (2024) showcase the presence of early KCD stars in the central region of high-redshift disc galaxies, but also suggest that the majority of stars within these primordial discs are unlikely to evolve into the present-day thin disc of spiral galaxies such as our own MW, as both angular momentum transport by a stellar bar and the presence of several massive star clusters that can sink due to dynamical friction underscore the tendency of stars to move inwards. This offers compelling indications that stars within early KCDs experience gradual migration away from their birth locations, a phenomenon known as ‘radial migration’.

Interestingly, recent investigations are uncovering the existence of old metal-poor stars, some even reaching ultra metal-poor levels ([Fe/H] |$< -4$|⁠), with large angular momenta, primarily located close to the Galactic plane and following dynamically cold disc orbits (e.g. Sestito et al. 2019, 2020; Cordoni et al. 2021; Sestito et al. 2021; Dovgal et al. 2024). Notably, observations with GAIA (by, e.g. Nepal et al. 2024) have even come so far to conclude that the MW, similar to the high-redshift galaxies observed by JWST, has an old KCD. The stars in this old disc range from metal-poor to super metal-rich metallicities, suggesting a history of intense star formation (SF) and the early establishment of cold discs (see also Fernández-Alvar et al. 2024). Nevertheless, the existence of an ancient, metal-poor component within the MW’s thin disc remains a highly debated topic. For example, Zhang, Ardern-Arentsen & Belokurov (2024) used a Gaia Radial Velocity Spectrometer sample with a broad metallicity range and found no significant population of very metal-poor stars on disc-like orbits within 2.5 kpc the Galactic plane, suggesting that previous observations might actually reflect the prograde halo component instead. Conversely, Bellazzini et al. (2024) showed a bimodal distribution in the vertical angular momentum distribution of prograde stars, indicating the presence of two disc-like components.

Recent studies utilizing JWST have also demonstrated the presence of bulge-dominated galaxies at high redshift (e.g. Tsukui & Iguchi 2021; Huertas-Company et al. 2024; Tsukui et al. 2024), but the main drivers of bulge growth remain largely unknown. Guedes et al. (2013) employed Eris, a ‘zoom-in’ cosmological simulation of a late-type galaxy, to explore the formation of its pseudo-bulge. They noted that the majority of the pseudo-bulge’s mass originated from a bar configuration at early times, which subsequently transformed into a compact flattened structure and an inner bar. Similarly, Okamoto (2013) examined the formation pathways of two pseudo-bulges within hydrodynamical simulations and showed that both pseudo-bulges formed at high redshift through the accretion of misaligned gas, with secular evolution contributing less than 30 per cent to their final mass. Gargiulo et al. (2019) studied the galactic bulges in the Auriga (Grand et al. 2017) simulations and found that pseudo-bulges are formed predominantly in situ. JWST observations and simulations both seem to suggest that the formation of the primordial discs and bulge occurred within a similar time frame.

In this paper, we examine the fate of primordial KCDs in a MW-sized galaxy. Specifically, we aim to demonstrate that a significant portion of stars born within early primordial thin discs undergo inward migration, resulting in structures with smaller radii, and hinting at the fact that the ‘inside-out’ formation scenario of the present-day MW thin disc does not occur at high redshift. Eventually, we argue that these migrated stars integrate into other substructures within the galaxy, particularly the pseudo-bulge. Additionally, we argue that a subset of stars born at the outer edges of the primordial KCDs may persist within the present-day thin disc. Here, we employ an exceptionally high-resolution, hydrodynamical, cosmological ‘zoom-in’ simulation (GigaEris; Tamfal et al. 2022). The plan is as follows: Section 2 provides a concise overview of the simulation setup. In Section 3, we present the results of the simulation, focusing on the trajectory and final location of stars that were originally ‘born within the KCDs’. Finally, we discuss our findings in Section 4.

2 METHODS

2.1 The GigaEris simulation

The N-body/hydrodynamic cosmological ‘zoom-in’ simulation GigaEris was presented in Tamfal et al. (2022). The simulation was run using ChaNGa (Jetley et al. 2008, 2010; Menon et al. 2015), an N-body smoothed-particle hydrodynamics code. The simulation relies on a low-resolution dark matter (DM)-only run, which was used to identify a Galactic-scale halo at |$z = 0$|⁠. The halo was chosen to closely match the mass of the MW’s halo and exhibit a rather quiet late merging history, following a methodology similar to that of the original Eris suite (Guedes et al. 2011). The halo was then re-simulated at higher resolution. The initial conditions were generated with the music code (Hahn & Abel 2011), with 14 levels of refinement and the cosmological parameters |$\Omega _{\rm m}$| = 0.3089, |$\Omega _{\rm b}$| = 0.0486, |$\Omega _{\Lambda }$| = 0.6911, |$\sigma _8$| = 0.8159, |$n_{\rm s}$| = 0.9667, and |$H_0$| = 67.74 km s|$^{-1}$| Mpc|$^{-1}$| (see Planck Collaboration XIII 2016).

We set the gravitational softening for all particles as a constant in physical coordinates (⁠|$\epsilon _{\rm c} = 0.043$| kpc) for redshifts lower than |$z = 10$|⁠, and scaling as |$11\epsilon _{\rm c}/(1+z)$| at earlier times. The final snapshot, corresponding to |$z=4.4$|⁠, contains the following particle counts: |$n_{\rm DM} = 5.7 \times 10^8$|⁠, |$n_{\rm gas} = 5.2 \times 10^8$|⁠, and |$n_{\star } = 4.4 \times 10^7$|⁠. Star particles, each with an initial mass of |$m_{\star } = 798$| M|$_{\odot }$|⁠, were generated based on gas density and temperature criteria (Stinson et al. 2006). The SF rate depends on the gas density, gas dynamical time, and a fixed SF efficiency. The code calculates non-equilibrium abundances and cooling for hydrogen and helium, including gas self-shielding and a redshift-dependent radiation background (Pontzen et al. 2008; Haardt & Madau 2012). Additionally, cooling from metallic elements’ fine structure lines is computed using rates from Cloudy (Ferland, Fabian & Johnstone 1998; Ferland et al. 2013) and the methodology of Shen, Wadsley & Stinson (2010) and Shen et al. (2013), assuming no self-shielding (for a discussion, see Capelo et al. 2018).

Feedback from Type Ia and Type II supernovae (SNae) injects energy, mass, and metals into the surroundings (Thielemann, Nomoto & Yokoi 1986; Stinson et al. 2006). The SN Type Ia feedback is independent of progenitor mass, whereas SN Type II feedback follows a delayed-cooling recipe (Stinson et al. 2006). For each SN Type II event, gas, oxygen, and iron mass, dependent on star mass, are injected into the surrounding gas (Woosley & Weaver 1995; Raiteri, Villata & Navarro 1996). Stars with masses 8–40 M|$_{\odot }$| undergo SN Type II, whereas those with masses 1–8 M|$_{\odot }$| release mass as stellar winds, with the metallicity of the wind gas being the same as that of the star particle.

2.2 Identifying substructures

2.2.1 Primordial kinematically cold discs

We examine the evolution of the KCDs by initially decomposing the galaxy into its components at different redshifts. In the study by Tamfal et al. (2022), multiple stellar discs were identified within the GigaEris simulation at |$z>4.4$|⁠, emerging at different times and along different, misaligned accretion planes. The discs were identified using the DBSCAN clustering algorithm (Sander et al. 1998) for stars born within the last 17 Myr within a 4 kpc box surrounding the galaxy’s centre.² In this work, we simply define KCDs as those cold stellar components having |$< v_{\phi }>$||$/$||$\sigma _{\mathrm{ R}}$| greater than unity and persisting for at least 150 Myr from the moment their first stars are born. We identify 11 distinct stellar KCDs, each with a unique stellar population, that formed at various times and along misaligned accretion planes (see Fig. 1), the oldest of such discs originating at |$z \sim 7.2$|⁠.

It is worth highlighting the compact radial extent of such discs, spanning no more than 1 kpc (see Section 3.1) at these early epochs. The compact size reflects the small virial radius of the main host halo (⁠|$\sim$|37 kpc at |$z = 4.4$|⁠) which, in turn, imposes a limit on the orbital angular momentum of baryonic matter that can gravitationally bind to the halo and subsequently cool to form a disc (White & Rees 1978; Efstathiou & Jones 1980). Consequently, the overall aspect ratio (scale height to disc scale length) of the identified stellar KCDs is larger compared to their present-day counterparts, exceeding for example the aspect ratio of the present-day MW thin disc, which is estimated to be between 0.1 and 0.25 (e.g. Bland-Hawthorn & Gerhard 2016). We also observe a decreasing trend in the aspect ratio of these structures. For instance, the primordial KCD that formed at |$z \sim 7.2$| exhibits an initial aspect ratio of 0.49, whereas the disc formed at |$z \sim 4.9$| has an initial aspect ratio of 0.32. This trend is accompanied by a shift in kinematics, as depicted in Fig. 1, which shows that the later discs have a higher |$< v_{\phi }>$||$/$||$\sigma _{\mathrm{ R}}$| ratio when born.

2.2.2 Pseudo-bulge

We use a combination of spatial and kinematical criteria, similar to those employed by Sokolowska et al. (2017) and Gargiulo et al. (2019), to define the pseudo-bulge at the final snapshot of the simulation, |$z=4.4$|⁠.

First, we consider particles located inside a sphere of radius |$r_{\rm bulge} = 2 R_{\rm eff}$|⁠, where |$R_{\rm eff}$| is the effective radius of the bulge (enclosing half of the total projected light in a given band), determined here by fitting a Sérsic (1963) profile to the face-on surface brightness in the V band. We employ a non-linear least-squares approach to fit the sum of an exponential and a Sérsic function, expressed in terms of projected radiation intensity profile, as

where |$I_{\rm e}$| is the intensity of the bulge at |$R_{\rm eff}$|⁠, n is the Sérsic index, and |$b_{\rm n}$| satisfies the equation |$\Gamma (2 \rm {\mathit{ n}}) = 2 \gamma (2 \rm {\mathit{ n}}, \mathit{ b}_{\rm \mathit{ n}})$|⁠, where |$\Gamma$| and |$\gamma$| are the Gamma function and lower incomplete Gamma function, respectively (to ensure that |$R_{\rm eff}$| contains half the projected light). Additionally, |$I_{\rm 0}$| stands for the central intensity of the disc component, and h denotes the disc scale length with value |$h=1.12$| kpc. We find a Sérsic index of 0.7 (and |$b_{\rm \mathit{ n}} = 1.08$|⁠). This classifies our bulge as a pseudo-bulge, since pseudo-bulges typically have indices below 2, whereas classical bulges have indices above 2, with minimal overlap between the two categories (see, e.g. Fisher & Drory 2008). We find a value of 0.14 kpc for |$R_{\rm eff}$|⁠.

Secondly, we use a circularity parameter, |$\epsilon$|⁠, derived from the z-component of the angular momentum vector (⁠|$J_z$|⁠) of stellar particles within the galactic disc in the x–y plane, compared to the angular momentum (⁠|$J_{\rm c}$|⁠) of a hypothetical particle on a circular orbit: |$\epsilon = J_z/J_{\rm c}$|⁠. In a typical spiral galaxy, the distribution of circularities – |$f(\epsilon) = \Delta N/\Delta \epsilon$|⁠, where |$\Delta N$| is the number of particles in the circularity bin |$\Delta \epsilon$| around |$\epsilon$| – is expected to exhibit two prominent peaks, one around |$\epsilon \approx 1$|⁠, corresponding to the disc, and the other at |$\epsilon \approx 0$|⁠, for the bulge. To define the pseudo-bulge, we exclude stellar particles with pure disc kinematics as well as those in the outer halo. For that purpose, we consider as pseudo-bulge stars only those particles with circularities |$|\epsilon | \le 0.7$|⁠, consistent with Gargiulo et al. (2019), and within a distance |$R_{\rm eff}$| from the halo centre. The circularity of all the stars within a central cylinder of radius 2 kpc and height 0.5 kpc in the main galaxy halo at the final snapshot is plotted in Fig. 2.

3 RESULTS

To determine the fate of the primordial KCDs, we start by looking into the distances from the galactic centre³ of stars both at birth and in the final snapshot. Specifically, we focus on the discs formed at redshifts |$z \sim 7.2$|⁠, 5.8, 5.3, and 4.9 and plot in Fig. 3 the two distances for all the stars born in such structures. The red dashed lines mark the locus where the birth and final radii are equal, separating outward migration for stars above the line and inward migration for those below it. We utilize 150 bins (evenly sized in logarithmic space) for both radii, spanning from 0.001 to 5 kpc, and implement a strict threshold of at least 50 star particles per bin to make sure we recover general trends and exclude outliers. The feature near |$R_{\rm {birth}} \sim 1.8$| kpc in the top panel of Fig. 3 arises from the binning, due to the oldest KCD forming more stars in the outer regions of the disc compared to younger discs. As this small structure does not impact the understanding of the KCD’s fate, it was not analysed in detail.

Intriguingly, KCD stars born at |$z \sim 7.2$| exhibit predominantly inward migration. This pattern starts shifting at |$z \lesssim 5.8$|⁠, when KCD stars begin displaying a more mixed radial migration activity, involving both inward and outward displacements (see Fig. 4). The peak in concentration of star particles moves closer to the red dashed line for disc stars born at |$z \sim 5.3$| and 4.9. Stellar radial migrations are likely driven by interactions with non-axisymmetric features of the galaxy, such as the bar that is present in GigaEris (as discussed in Tamfal et al. 2022; van Donkelaar et al. 2024). These features may induce ‘churning’ (Sellwood & Binney 2002; see also Willett et al. 2023 for its connection to the thin disc), wherein an exchange of angular momentum redirects stars into new orbits. Additional research is however necessary to ascertain the primary driver behind the radial migration of stars in early KCDs.

3.1 Link to the formation of the pseudo-bulge

We can delve deeper into the information presented in Fig. 3 by examining the distribution of stars within the x–z plane⁴ of primordial KCDs at birth and at the final snapshot of the simulation at |$z = 4.4$|⁠, as depicted in Fig. 5 for the discs born at |$z \sim 7.2$| (top panels) and |$z \sim 4.9$| (bottom panels). Combining the motion of stars shown in Fig. 3 with Fig. 5, we provide the first compelling hints that primordial KCDs at |$z \gtrsim 6$| become part of the central region of the galaxy, whereas later discs contribute more to the thick disc. The upper-right panel of Fig. 5 displays the position of stars originating from a KCD born at |$z \sim 7.2$| in the final snapshot of the simulation (⁠|$z=4.4$|⁠). This panel clearly reveals a clustering of stars at the galaxy’s centre, with the majority congregating within a radius of 250 pc. In contrast, stars stemming from a KCD born around |$z \sim 4.9$| maintain a more disc-like arrangement, though with an increased thickness and shorter scale radius, as shown in the bottom-right panel. This suggests that the primordial KCDs formed at |$z \gtrsim 6$| may play a notable role in shaping central structures like the pseudo-bulge.

However, the primordial stellar KCD at |$z \sim 7.2$| has the lowest mass of all the identified KCDs (⁠|${\approx} 10^{8.79}$| M|$_{\odot }$|⁠).⁵ Therefore, even though 76 per cent of the stellar mass of this structure (⁠|${\approx} 10^{8.67}$| M|$_{\odot }$|⁠) ends up in the pseudo-bulge, only about 7 per cent of the total mass of the pseudo-bulge identified at |$z=4.4$| (⁠|${\approx} 10^{9.85}$| M|$_{\odot }$|⁠) originates from this disc. For all identified primordial KCDs, this mass fraction ranges between 2 and 8 per cent, as shown by the blue crosses in Fig. 6. Consequently, approximately 40 per cent of the stellar mass of the identified pseudo-bulge at |$z=4.4$| comes from stars originating from all primordial KCDs that meet our criteria (right-most light-green square in Fig. 6), with all discs contributing more or less the same amount of stellar mass to the central structure. Still, the early primordial KCDs at |$z \gtrsim 6$| have a significant impact on the early stages of pseudo-bulge formation. Specifically, 38 per cent of the stars that are already formed at |$z \sim 7.2$| and will contribute to the pseudo-bulge originated from the first primordial KCD.

To further investigate this evolutionary scenario, we compute the proportion of the initial KCDs’ stellar mass that will ultimately reside in the pseudo-bulge at different formation redshifts using the identified star particles associated with the pseudo-bulge, as depicted in Fig. 6 (dark-green, dashed line). This figure shows us indeed that most of the stars born in the primordial KCDs at |$z \gtrsim 6$| become part of the pseudo-bulge. It presents a particularly intriguing result, indicating that the vast majority (approximately 76 per cent) of the stellar mass (⁠|${\approx} 10^{8.67}$| M|$_{\odot }$|⁠) in the stellar KCD formed at |$z \sim 7.2$| will transition to the pseudo-bulge by |$z=4.4$|⁠, whereas this proportion decreases to approximately 8 per cent for the stellar mass (⁠|${\approx} 10^{7.68}$| M|$_{\odot }$|⁠) in the KCD formed at |$z \sim 4.8$|⁠. However, as indicated by the light-green dot-dashed line, KCDs account for just under 50 per cent of the total mass in the pseudo-bulge. This suggests that the remaining stellar content in the pseudo-bulge is likely due to other processes, such as mergers or central SF.

It is evident from Fig. 6 that the mass fraction is not a monotonic function of time, as indicated by the peaks around |$z \sim 5.7$| and |$\sim 5.2$|⁠. At these redshifts, a significant portion of the stellar mass formed in the primordial KCDs transitions to the pseudo-bulge by |$z = 4.4$|⁠. Notably, these peaks closely align with two peaks in SF within the main galaxy halo between |$5 \lesssim z \lesssim 6$| (see the SF rate in fig. 3 of Tamfal et al. 2022). One plausible explanation is that stars born in primordial KCDs formed near peaks of SF experience greater radial migration as a result of the formation of the galactic bar. The correlation between the bar and central SF has been illustrated, for instance, by Fanali et al. (2015; see also, e.g. Martinet & Friedli 1997; Sakamoto et al. 1999; Sheth et al. 2005; Ellison et al. 2011; Athanassoula, Machado & Rodionov 2013; Zhou, Cao & Wu 2015; Marioni et al. 2022; Yu et al. 2022; Zee et al. 2023), who demonstrated that as the bar evolves, gas is funneled towards the central regions of the galaxy, thereby enhancing the central SF rate. A higher SF rate could therefore suggest a more pronounced bar and as a non-axisymmetric feature of the galaxy, the bar likely plays a role in radial migration through its interactions (e.g. Sellwood & Binney 2002). However, further investigation is needed to gain a better understanding of this phenomenon and the bar itself, especially as the SF peak around |$z \sim 4.8$| appears to lack a similar impact as the previous SF peaks.

As discussed, approximately 40 per cent of the stellar mass of the identified pseudo-bulge at |$z = 4.4$| comes from stars originating from all primordial KCDs that meet our criteria. Amongst the stars not originating from these discs, the vast majority, around 96 per cent, formed within the pseudo-bulge region, defined as the region within twice the effective radius (⁠|$2 R_{\rm eff}$|⁠), from gas brought in by the bar within the past 500 Myr. Therefore, one could conclude that the pseudo-bulge in the main galaxy halo of GigaEris is mainly formed in-situ, indicating that nearly all stellar particles were born within the virial radius of the main galaxy halo, whether in a KCD or in the pseudo-bulge region itself, aligning with findings from (pseudo-)bulges within the Auriga simulation (see Gargiulo et al. 2019).

3.2 Link to the present day’s discs

It is important to note that Figs 3 and 5 only illustrate the overall trend, as we have not depicted any outliers. The outward migration of stars born within a primordial KCD is common across all redshifts and birth radii, which contrasts with what is shown, for example, in the top panel of Fig. 3. Given that these stars likely experience less kinematical heating at the outer edges of the galaxy, they may still influence the formation of the older component of the thin disc observed today, as e.g. the stars discussed by Nepal et al. (2024).

Fig. 7 illustrates the birth radius plotted against the final radius at |$z=4.4$| for all stars born in primordial KCDs, with the colour bar indicating the properties of these stars. The red dashed lines signify where the birth and final radii are equal, indicating outward migration for stars above the line and inward migration for those below it. The first panel is similar to Fig. 3, yet in this instance, it shows all stars born in all primordial KCDs. This panel indicates that outward migration happens across all birth radii. However, it is notable that the peak of this plot lies below the red dashed line, suggesting that most stars born within the primordial KCDs tend to migrate inwards.

The second panel shows the [O/Fe] ratio⁶ of the stars born in the primordial KCDs. In the MW at |$z = 0$|⁠, thick-disc stars have a larger oxygen abundance than thin-disc stars with the same [Fe/H], as shown by, e.g. Franchini et al. (2021). In this panel, we observe that the [O/Fe] ratio decreases with a higher R|$_{\rm birth}$|⁠. Stars with the lowest [O/Fe] ratios were born outside 1 kpc and then migrated inwards or ended up outside 1 kpc as a result of outwards migration by |$z=4.4$|⁠. In general, a low [O/Fe] ratio indicates that the galaxy has undergone an extended or delayed SF history (e.g. Gallazzi et al. 2005). Consequently, it makes sense that the outer regions of the galaxy exhibit a low [O/Fe] ratio, as these areas typically experience more prolonged and less intense SF compared to the central regions. These low-[O/Fe] stars are in line with the stars observed by Nepal et al. (2024), indicating that these potentially older thin-disc stars exhibit a slightly lower [|$\alpha$|/Fe], akin to [O/Fe], enhancement compared to the thick-disc stars. Previous observations in small, high-resolution spectroscopic samples have also identified old stars with low [|$\alpha$|/Fe] ratios (see, e.g. Anders et al. 2018; Miglio et al. 2021; Gent et al. 2024).

The last two panels show the change in vertical and radial velocity dispersions (⁠|$\Delta \sigma = \sigma _{\rm end}-\sigma _{\rm start}$|⁠). As previously discussed, it is expected that stars in the thin disc have smaller velocity dispersions compared to those in the thick disc (e.g. Gilmore et al. 2002; Soubiran et al. 2003; Parker et al. 2004; Wyse et al. 2006; Mackereth et al. 2019). Therefore, for stars born in the primordial KCDs to remain as thin-disc stars, they should experience minimal kinematical heating. In both panels, we can cleary see that stars reaching the outer radius of 1 kpc, particularly those that migrated there, experience less heating. However, we observe that stars migrating outward to just below 1 kpc experience significant heating in only their vertical velocity dispersion. We do not yet have a clear explanation for this phenomenon, and further investigations in future studies are warranted. The minimal kinematical heating observed at the outer edges of our galaxy could be explained by the reduced turbulence experienced by stars in these regions, such as the declining effect of spiral arm and bar heating with radius (see, e.g. Grand et al. 2016; Mackereth et al. 2019).

This suggests that stars with a low [O/Fe] ratio and minimal kinematical heating within the simulation could be associated with the formation of the older component of the present-day thin disc, possibly representing the recently identified old stars with low [|$\alpha$|/Fe] ratios in the MW. The presence of a kinematically cold component at the outer edge of the disc suggests that mergers with companions may not have been as effective at kinematically heating the stars (see, e.g. Galán-de Anta et al. 2023a, b), given that several mergers occur – at |$z=6.52$|⁠, 6.01, and 4.53, with corresponding mass ratios 0.16, 0.14, and 0.12, respectively (Tamfal et al. 2022) – namely after the formation of the first primordial KCD at |$z \sim 7.2$|⁠. However, we cannot definitively rule out heating via mergers beyond |$z< 4.4$|⁠. None the less, the data from Nepal et al. (2024), who highlight the presence of old thin-disc stars, suggests that heating through mergers may indeed be minimal. In a galaxy characterized by a more turbulent history of mergers, these old thin-disc stars might not always persist till late times. However, our MW is believed to have experienced a relatively quiescent merger history (e.g. Fragkoudi et al. 2020; Kruijssen et al. 2020), making our analysis relevant.

Furthermore, as mentioned earlier, it is important to highlight that the primordial KCDs in our simulations have a characteristic radius of about 1 kpc. Considering that radial migration is expected to still occur beyond |$z \lesssim 4.4$|⁠, it adds uncertainty about where these surviving KCD stars will ultimately settle. It is reasonable to assume that, to preserve their cold kinematics, they would tend to end up near the present-day outer edge of the thin disc, where they are less affected by secular heating due to the galactic bar and spiral structure. Alternatively, as the current GAIA observations of Nepal et al. (2024) suggest, they might settle around 1 kpc from the galactic centre, within the bulge region, where they would experience little secular heating as this region is dominated by the spheroidal bulge potential, which suppresses the propagation of density waves and weakens the effect of the bar, whose major axis is much larger in scale (5–6 kpc).

4 DISCUSSION AND CONCLUSIONS

In this study, we have demonstrated that most stars within primordial KCDs, characterized as stellar discs with |$< v_{\phi }>/\sigma _{\mathrm{ R}}$| larger than unity for at least 150 Myr after birth, undergo migration and ultimately integrate into other substructures within the galaxy, with particular emphasis on the pseudo-bulge. Stars born in KCDs formed at |$z \gtrsim 6$| predominantly migrate to the pseudo-bulge. In contrast, stars born in KCDs formed after this redshift show signs that they might contribute more significantly to the thick disc. However, due to the limitations of the simulations, the fate of these stars at |$z \le 4.4$| remains uncertain. Nevertheless, stars originating from the outer edges of the discs, particularly those that have migrated even farther, exhibit minimal kinematical heating, suggesting a possibility of their survival and integration into the present-day thin disc.

We also remark that the KCDs observed in this work assemble already at |$z \sim 7$| and that the formation of a cold, stable disc at high redshift is possible without relying on the ‘inside-out’ and ‘upside-down’ formation processes described by Bird et al. (2013; see also Agertz et al. 2021; Bird et al. 2021; Belokurov & Kravtsov 2022). Consequently, it will not be surprising to observe stars with cold kinematics in galaxies at high redshift, thereby providing context to high-redshift disc observations, even though they are currently only photometric (e.g. Ferreira et al. 2023). However, as Fig. 7 illustrates, most of the stars in these early discs are expected to heat up over time, causing them to no longer be part of the KCDs. Furthermore, the discs themselves will likely grow as progressively higher angular momentum material is accreted. Alongside the bottom panels in Fig. 5, which demonstrates that later KCDs appear to form a more traditional thick disc, one can expect that both ‘inside-out’ and ‘upside-down’ formation processes play an increasingly significant role in shaping the present-day thin disc. We propose that the prevalence of these mechanisms increases significantly at |$z \lesssim 4$|⁠, extending beyond the final snapshot of our simulation. This results in a ‘two-phase’ scenario: stars formed in early, primordial KCDs primarily migrate inwards, creating new structures with smaller radii than their original formation sites. Nevertheless, at high redshift, these stars remain detectable as disc-like structures, and eventually pseudo-bulges. The present-day KCD, in contrast, forms most likely through an ‘inside-out’ and ‘upside-down’ mechanism following |$z \approx 4$|⁠.

This ‘two-phase’ thin-disc formation process introduces an interesting contrast between the dynamics characterizing the early and low-redshift stages of thin-disc evolution. As a consequence, it becomes important to exercise caution when drawing parallels between high-redshift discs observed with instruments like JWST (see, e.g. Ferreira et al. 2022, 2023; Kartaltepe et al. 2023; Robertson et al. 2023) and their counterparts at lower redshifts. Likewise, since this first phase of disc formation is fast, resulting in a transition to a pseudo-bulge in less than a billion years, the qualitative implication for high-redshift observations is that one should observe a rapid rise of the bulge-to-disc ratio of discy galaxies already at |$z> 4$|⁠. The timing of the transition would depend on the assembly history of the galaxy, in particular how rapidly the stellar mass can grow enough to become susceptible to internal instabilities such as bar formation, which promote radial migration and pseudo-bulge formation. Our single ‘zoom-in’ simulation should thus be seen as a numerical experiment exposing this new evolutionary scenario, while future work will have to assess its quantitative implications and observational diagnostics in a galaxy population at high redshift, by combining information from additional high-resolution ‘zoom-in’ simulations and from cosmological volumes.

However, it is essential to acknowledge again that Figs 3 and 5 present the general trend, without depicting any outliers. The outward migration of stars originating within a primordial KCD is observed across all redshifts, although it may be less frequent as shown by the first panel of Fig. 7. As seen in the third and final panel of Fig. 7, it is plausible to conclude that the stars that migrated to the outer edge of the galaxy experienced less kinematical heating. This suggests that these stars might have had an influence on the formation of the older component of the thin disc observed today and could potentially represent the previously detected old stars with low [|$\alpha$|/Fe] ratios (see, e.g. Anders et al. 2018; Miglio et al. 2021; Gent et al. 2024; Nepal et al. 2024). These stars migrated outwards, aligning with the conventional ‘inside-out’ disc formation theory. This simultaneous outward migration, alongside inward migration, provides additional support for the proposed ‘two-stage’ disc formation process proposed in this paper.

A key concern might be whether the observed changes in radial migration patterns reflect genuine physical properties of galaxy evolution or if they are simply the result of poorer numerical resolution at high redshift. However, despite the challenges of resolving high-redshift discs, the GigaEris simulation has exceptionally high baryonic resolution, i.e. an initial stellar mass of |$m_{\star } = 798$| M|$_{\odot }$|⁠. As mentioned, the lowest mass disc in our sample contains more than |$10^5$| stellar particles, which significantly reduces the risk of resolution-driven biases. Therefore, we interpret the redshift-dependent migration changes as intrinsic to the galaxy’s evolution.

Finally, further investigation is crucial to understand the shift occurring around |$z \lesssim 6$|⁠. Before this time, most of the stellar mass within the primordial KCDs migrates to the pseudo-bulge region, whereas afterwards less than half of the total stellar mass from the formed primordial KCDs goes towards the pseudo-bulge. This trend intensifies around |$z \lesssim 5$|⁠, with only about 15 per cent of the stellar mass of the KCD formed at |$z \sim 5.0$| ending up in the pseudo-bulge. One potential explanation for this trend is that at redshift |$z \lesssim 6$| the atomic neutral hydrogen gas fraction (⁠|$F_{\rm gas, H\,{\small{I}}}$|⁠) in the galaxy reaches approximately 25 per cent. As the redshift decreases further to |$z \lesssim 5$|⁠, |$F_{\rm gas, H\,{\small{I}}}$| in the galaxy halo drops below 20 per cent. Interestingly, this specific gas fraction has been associated with the conditions necessary for a disc to form that exhibits kinematics similar to those observed in the present-day thin disc, as shown using isolated disc galaxy simulations in van Donkelaar et al. (2022). In particular, when the gas fraction is high, newly formed stars in disc galaxies exhibit high vertical velocity dispersions. For instance, they showed that galaxies with a gas fraction |$f_{\rm gas} \gtrsim 30$| per cent have a vertical velocity dispersions of about 30 km s|$^{-1}$|⁠, which aligns more with the kinematics of the MW’s thick disc rather than its thin disc (Mackereth et al. 2019). In contrast, lower gas fractions (⁠|$f_{\rm gas, H\,{\small{I}}}\lesssim 20$| per cent) result in much smaller dispersions, under 10 km s|$^{-1}$|⁠, which are consistent with the kinematics of young stars in the MW’s thin disc. This topic is currently being investigated for an upcoming follow-up paper.

ACKNOWLEDGEMENTS

PRC acknowledges support from the Swiss National Science Foundation under the Sinergia Grant CRSII5_213497 (GW-Learn). LM and FvD acknowledge support from the Swiss National Science Foundation under the grant 200020_207406. Support for this work was provided also by NASA through grant 80NSSC21K027 (PM). The simulations were performed on the Piz Daint supercomputer of the Swiss National Supercomputing Centre (CSCS) under the project id s1014.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES