EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 709 (May 2026) A&A, 709 (2026) A131 Full HTML
Open Access
Issue
A&A
Volume 709, May 2026
Article Number A131
Number of page(s) 14
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202557661
Published online 08 May 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

The Milky Way (MW) stellar disc exhibits a well-defined chemical bimodality, with stars occupying two distinct sequences in the [α/Fe]-[Fe/H] plane (Fuhrmann 1998; Bensby et al. 2014; Hayden et al. 2015; Imig et al. 2023; Guiglion et al. 2024). The stars in each of these sequences display a range of kinematic and spatial trends that represent the end product of the Galaxy’s formation and evolution. Understanding the origin of these features can therefore constrain the physical processes that shape the evolution of galactic discs, including gas flows, radial migration, and feedback mechanisms.

Several scenarios have been proposed to explain the observed bimodality. For instance, the results of infall chemical evolution models indicate that distinct episodes of gas accretion and star formation have the potential to establish different abundance sequences (e.g. Chiappini et al. 1997; Spitoni et al. 2022).

However, recent works have found that infall scenarios of pristine gas face strong constraints, as they must reproduce the α-bimodality without any excessive dilution of the metallicity of the interstellar medium, since this type of dilution would prevent these models from matching the observed age-abundance trends (Dubay et al. 2025). Another hypothesis emphasises the role of stellar radial migration, in which stars originating in regions with different chemical enrichment histories can mix throughout the disc due to changes in the stars’ specific angular momentum through non-axisymmetric forces, such as those imparted by the bar or spiral arms (e.g. Sellwood & Binney 2002; Schönrich & Binney 2009; Minchev et al. 2013). A third scenario explains the chemical dichotomy as the outcome of two distinct modes of star formation. On the one hand, in the model of Clarke et al. (2019), the formation of the high and low-α sequences is not sequential. The former forms in high-density star formation rate (SFR) clumps that dominate the early galaxy, while the latter forms in spatially extended regions of low-SFR density. On the other hand, Khoperskov et al. (2021) showed that two star formation regimes arise naturally in isolated MW–like simulated galaxies as a consequence of an initially rapid dissipative gas collapse that drives a high SFR, builds the high-α sequence, and triggers strong stellar feedback that expels gas from the disc. The subsequent re-accretion of this gas on longer time-scales sustains a more quiescent star formation phase, giving rise to an extended low-α sequence. These two regimes are supported by observational reconstructions of the MW’s star formation history (e.g. Snaith et al. 2014).

Furthermore, the prevalence of α-bimodalities in external disc galaxies is also unknown. For instance, the question of whether a chemical bimodality exists in M31 remains a matter of debate, as the currently available abundance measurements are sparse and provide inconclusive or model-dependent evidence for different α sequences (Kobayashi et al. 2023; Nidever et al. 2024). From a theoretical point of view, the identification of bimodalities in simulations depends on how the bimodality is defined in the abundance space and on the adopted analysis methodology, with no currently agreed-upon, consistent criterion applied across different simulation suites. As a result, it remains unclear whether this feature is a universal outcome of MW-mass galaxy formation models (e.g. Grand et al. 2018; Mackereth et al. 2018; Parul et al. 2025). All this raises the question of whether the α-bimodality is a ubiquitous feature of disc galaxies or an outcome of the particular formation history of the MW, highlighting the importance of establishing its presence in external systems (e.g. Pinna et al. 2024).

In this paper, we analyse the role of gas accretion and radial gas flows in creating a chemical dichotomy in a simulated hydrodynamical, zoom-in MW-mass galaxy. We quantify the fluxes of gas across the disc, examine their relation to changes in the star formation and merger history, and their origin and interplay with the galactic corona. We investigate how these properties contribute to the observed chemical and kinematic evolutionary trends of the stellar disc. For this purpose, we adopted a simulated MW-mass galaxy with a clear α-bimodality, which has no bar and is characterised by a quiet merger history. The layout of this paper is as follows.

In Section 2, we describe the general properties of the simulation we use in this paper. Section 3 details the spatial, kinematical, and chemical properties of its stellar disc. Section 4 then discusses the physical mechanisms that drive these properties and, specifically, the role of gas accretion and radial gas flows, along with the chemical evolution of the galactic corona and Romeo’s merger history. Lastly, we discuss our results and present the conclusions in Section 5.

2 Romeo in context

2.1 Simulation details

We analysed the simulated Romeo galaxy (Garrison-Kimmel et al. 2019), part of the Exploring the Local Volume in Simulations (ELVIS) FIRE-2 cosmological zoom-in suite (Hopkins et al. 2018), which models galaxy pairs similar to those residing in the Local Group. The centre of Romeo is defined as in Wetzel et al. (2023). An iterative zoom-in approach is applied to the star particles: starting from their mean centre-of-mass, the enclosing sphere is repeatedly shrunk by 50% until its radius falls below ∼10 pc. The systemic velocity is then taken as the centre-of-mass velocity of star particles within 8 kpc.

Romeo is the more massive of the pair and its disc formed earlier than those of the other galaxies in the suite. In particular, it started forming at 11 Gyr in look-back time, at odds with other FIRE-2 galaxies, where disc settling took place from about 9 Gyr ago onwards (McCluskey et al. 2024), but comparable with the MW disc, which is thought to have settled 10–13 Gyr ago (Bensby et al. 2003; Haywood et al. 2013; Belokurov & Kravtsov 2022; Xiang & Rix 2022; Gallart et al. 2024). As we will see, the stellar disc in Romeo exhibits an α-bimodality, which is not unusual within FIRE-2 galaxies, since 8 (including Romeo) of the 11 simulated galaxies analysed by Parul et al. (2025) were identified as exhibiting bimodality. We note, however, that the number of galaxies identified as exhibiting α-bimodality in Parul et al. (2025) might be higher, since the definition of bimodality adopted is conservative. In addition, Romeo only develops a bar during a short period at the very end of its evolution (Ansar et al. 2025), after the onset of the α−bimodality. This means that most of its evolution is unaffected by bar-driven radial mixing of stars, allowing us to isolate the role of gas flows in shaping the chemical evolution of its stellar disc.

The simulation was run using the GIZMO code in mesh-less finite-mass (MFM) mode (Hopkins 2015) and using the FIRE-2 galaxy formation model (Hopkins et al. 2018). An overview of the details of this model can be found in Hopkins et al. (2018); Wetzel et al. (2023). Romeo was simulated in a flat ΛCDM cosmology with the following parameters: Ωm = 0.31, ΩΛ = 0.69, Ωb = 0.048, h = 0.68, and σ8 = 0.82, ns = 0.97. Dark matter particles have a mass of MDM = 1.9 × 104 M⊙, while the initial mass of star particles and gas cells is Mbar = 3.5 × 103 M⊙. Star particles were set to lose mass over time through stellar winds, transferring mass to the gas. Dark matter and stars have gravitational softening lengths ϵDM = 32 pc and ϵstar = 4.4 pc, respectively. The softening length of the gas is adaptive and reaches a minimum value of ϵgas,min = 0.7 pc. The simulation consists of 600 snapshots in the range z = 99, t = 0 Gyr to z = 0, t ≃ 13.73 Gyr, with a time spacing of ≲25 Myr. In our analysis, we used the 39 snapshots that are publicly available and have an average and maximum time spacing of ∼0.3 Gyr and 1.57 Gyr, respectively.

2.2 General properties

Table 1 displays several global characteristics of the Romeo simulated galaxy and the MW. While the virial DM masses1 of both galaxies are similar within observational uncertainties, Romeo exhibits greater stellar and gas masses compared to the MW. However, it is important to note that the MW DM halo mass remains uncertain, within a factor of at least two (Bland-Hawthorn & Gerhard 2016; Karukes et al. 2020).

Romeo recently experienced a brief bar episode from roughly 900 to 100 Myr ago (Ansar et al. 2025), but since the α-bimodality is already established before the bar’s onset, our analysis focuses on the period preceding its formation, that is in the first 12.39 Gyr of the galaxy up to redshift z ∼ 0.1. Table 1 includes values for both z ∼ 0.1 and z = 0, showing that the global properties of Romeo remain similar between these two epochs. This indicates that extending the time frame to include the full 13.73 Gyr evolution of the galaxy does not alter the overall picture. In addition, Romeo has had a quiet merger history with no major mergers. Specifically, it did not undergo any mergers with galaxies exceeding 1% of its own stellar mass after z = 3, with the ratio calculated as the peak satellite stellar mass over the stellar mass of the host at the time of the merger.

The top panel of Fig. 1 compares the stellar surface density profile of Romeo with those of the Vintergatan simulated galaxy (Agertz et al. 2021) and the observed MW. Romeo exhibits the most extended disc with an exponential scale-length of 4.5 kpc at z = 0 and 4.1 kpc at a lookback time of 1.34 Gyr (right before the bar period). This is larger than the MW’s estimated value of 2.6 kpc (Bland-Hawthorn & Gerhard 2016), or the value of 3.26 ± 0.25 kpc, most recently estimated using Gaia DR3 (Khanna et al. 2025). Romeo is also less concentrated as evident from Romeo’s lower surface density up to ∼15 kpc. Similarly to the MW (see solid orange line in top panel of Fig. 1 and Lian et al. 2024), Romeo has a broken surface density profile at 1.34 Gyr in lookback time, although the break radius occurs at larger galactocentric distances.

The bottom panel of Fig. 1 depicts the global star formation history of Romeo. We determined the SFR by summing the masses of stars born within a cylinder of a radius of R < 40 kpc and height of |z| < 10 kpc at a chosen timestep, ∆t, with the mass of each star measured at the moment of its birth. To evaluate the burstiness of the SFR, we used ∆t = 10 Myr, as this corresponds to the typical snapshot time spacing (Wetzel et al. 2023), providing the best quantification of the instantaneous SFR without being affected by stochasticity. For a direct comparison with the Vintergatan simulated galaxy (Agertz et al. 2021), we adopted ∆t = 100 Myr. Finally, to smooth the SFR and capture the overall trend, we used ∆t = 500 Myr. The star formation history (SFH) of Romeo exhibits an early bursty phase characterised by an overall rise in the SFR, followed by quiescent SF at later times when the SFR declines, with a mild increase closed to the present epoch. These phases in the SFH are directly connected to the different evolutionary stages of the disc and the development of the α-bimodality, as we discuss below, and they broadly resemble the behaviour seen in the Vintergatan simulation.

Table 1

Global properties of the Romeo simulated galaxy at z ∼ 0.1, just before the start of its bar episode, and at z = 0, compared with those of the MW.

3 Disc properties

3.1 Gaussian Mixture modelling of components

We note that the term disc is used at two levels in this work. At a global level, it refers to the entire stellar disc of the simulated galaxy, Romeo. At a more fine-grained level, the global stellar disc is subdivided into distinct components purely based on differences in the age and kinematical properties of its star particles. To identify these disc components, we first performed a Gaussian mixture model (GMM) analysis, which allows for stellar particles to be classified in a probabilistic manner, avoiding hard cuts and being agnostic about the number of distinct components (see also Nikakhtar et al. 2021). Thus, can carry out a direct comparison of the best-fit model with traditional kinematic and stellar-abundance divisions.

We applied a GMM using 3D cylindrical velocities and the ages of star particles within R < 40 kpc and |z| < 10 kpc at z = 0. Specifically, the probability density of these four properties for each individual star particle in the volume of interest, di = (vR, vϕ, vz, A)i ∈ ℝ4, is given by a sum of N Gaussians, p(di|{ewn,μn,Σn}n=1N)=∑n=1NewnN(di|μn,Σn),Mathematical equation: $p({d_i}|\{ {e^w}_n,{\mu _n},{\Sigma _n}\} _{n = 1}^N) = \sum\limits_{n = 1}^N {{e^w}_n{\cal N}({d_i}|{\mu _n},{\Sigma _n}),} $(1)

with a weight of ew n ∈ ℝ4, mean µn ∈ ℝ4, and covariance matrix Σn ∈ ℝ4×4. Including the stellar age alongside cylindrical velocities allows particles with similar present-day orbital properties to be distinguished by their formation time, thereby enabling the separation of populations formed by different channels. We allowed N to vary from 1 to 16 and used the elbow rule and the deceleration in the decline of Bayesian information criterion (BIC) values to identify N = 4 as the optimal number of components (see Appendix A). For each GMM-identified component, we assigned star particles based on their probability of belonging to a specific Gaussian cluster. In particular, star particles in a cluster, k, satisfy ew^kN(di|μ^k,Σ^k)∑n=14ew^nN(di|μ^n,Σ^n)≥0.68,Mathematical equation: ${{{e^{\hat w}}_k{\cal N}({d_i}|{{\hat \mu }_k},{{\hat \Sigma }_k})} \over {\mathop \sum \nolimits_{n = 1}^4 {e^{\hat w}}_n{\cal N}({d_i}|{{\hat \mu }_n},{{\hat \Sigma }_n})}} \ge 0.68,$(2)

where ew^nMathematical equation: ${e^{\hat w}}_n$, μ^nMathematical equation: ${{\hat \mu }_n}$, and Σ^nMathematical equation: ${{\hat \Sigma }_n}$ are the estimated weight, mean and covariance matrix of each of the four Gaussian components. These values are inferred by maximising the likelihood given by Equation (1). As the initial values for the parameters are assigned randomly, it is necessary to perform multiple initialisations until a consistent result is reached. We found that seven initialisations were required to reach a stable result for our dataset.

After visual inspection of the kinematics of each identified component, we concluded that GMM identifies a halo and three disc-like structures, which, as we discuss in Section 3.2, can be associated with high-, intermediate-, and low-α components of the stellar disc. The age distribution of each of the three discs is shown in Fig. 2, while the average kinematical properties, along with MW’s disc properties, are reported in Table 2. In addition, Fig. 3 depicts the distribution of rotational velocities as a function of age and metallicity. As expected, the average rotational velocity decreases with increasing age of the component. Surprisingly, the coherent rotation of the high-α disc increases in the last 1.34 Gyr of evolution of the galaxy, the period within which a weak bar develops in Romeo. This together with the change in the vertical distribution of this profile during the bar development in Romeo (see bottom panel of Fig. 5) might be driven by the bar.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Top: stellar surface density profiles of the Romeo and Vinter-gatan (Agertz et al. 2021) simulated galaxies, together with the MW profile taken from Lian et al. (2024) and the best-fitting morphology of McMillan (2011). The former MW profile has been normalised to enclose a stellar mass of 4×1010 M⊙. The dot-dashed orange line depicts a falling exponential profile with scale-length of 2.6 kpc. Bottom: star formation rate as a function of time in Romeo calculated with timesteps ∆t = 10 Myr (orange), ∆t = 100 Myr (purple) and ∆t = 500 Myr (blue).
Table 2

Mean rotational velocity, v¯φMathematical equation: ${{\bar v}_\varphi }$, and velocity dispersions in the cylindrical R and z coordinates for the high-α, bridge-, and low-α sequences.

Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

1D age distribution of all disc star particles along with the disc components identified by GMM. The vertical grey lines mark key evolutionary stages of the disc identified by GMM (see Sect. 3.2) from the onset of the thick or high-α disc (10.9 Gyr ago) to the time when the bridge or intermediate-α disc begins to dominate over thick disc (7.7 Gyr ago) and the onset of the α−bimodality (3.6 Gyr ago). The vertical grey bands mark the timing of the bar episode.
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

2D distribution of the azimuthal or rotational velocity of all stellar particles in Romeo, R < 40 kpc and |z| < 10 kpc, as a function of age (left panel) and metallicity [Fe/H] (right panel), with the median trend shown in red. The green, blue, and purple lines depict the median trend of the stellar population in the high-α, bridge and low-α discs, respectively. In the left panel, the dashed vertical lines indicate the inferred period at which the MW’s stellar disc settles (see Sect. 3.2), while the dashed black line in the right panel shows the inferred trend for the MW as in Belokurov & Kravtsov (2022).

3.2 Disc components

Interestingly, while the GMM analysis is based solely on ages and kinematics, the identified components also exhibit clear differences in their chemical compositions (Fig. 4) and spatial distributions (Fig. 5), highlighting their distinct evolutionary histories. Because of the clear location on the [Mg/Fe]-[Fe/H] plane of the three disc components, we named them the high-α, bridge, and low-α discs. Furthermore, the three discs match the main disc’s formation phases in Romeo,

  • The high-α disc phase marks the onset of coherent rotation and the settling down in the inner galaxy of a radially compact disc with high velocity dispersion. As shown in Figs. 2 and 3, the spin-up phase timescale is similar to what is expected in the MW, but occurs at younger ages, between 8 and 11 Gyr. This is consistent with previous analyses of disc formation in Romeo (Yu et al. 2021, 2023; McCluskey et al. 2024). Notably, the thick disc formation coincides with an increase in the bursty SFR.

  • As the SFR drops and transitions from bursty to steady at around 7–8 Gyr ago (see bottom panel of Fig. 1), consistent with the estimate in Parul et al. 2023, a thinner intermediate-α and more radially extended disc is able to develop (see the bottom central panel of Fig. 5 compared to the left panel). This disc could correspond to the bridge region identified in the MW (Ciucă et al. 2021).

  • Finally, about 3.6 Gyr ago (at z ∼ 0.3), a mild increase of the overall SFR, coupled with a dilution of the gas in the disc as we shall see in Sect. 3.4, enables the onset of the disc bimodality and the growth of the low-α sequence. This timing agrees with that inferred by Parul et al. (2025) using a different methodology.

The panels in Fig. 4 show the number density of a sample of giants from Gaia DR3 cross-matched with APOGEE DR17, along with the number density and age distribution of disc stars, as identified by GMM, in the [Mg/Fe]–[Fe/H] plane at z = 0.0998. The Gaia and APOGEE data preparation and selection criteria is described in Appendix B. It can be seen from this figure that the distribution of star particles in Romeo, as in the case of the MW, is bimodal in [Mg/Fe]. However, in the MW, the α-sequences are broader and more clearly separated along the [Mg/Fe] axis. In particular, in Romeo, the difference in the peaks in the bimodal 1D [Mg/Fe] distribution at different galactocentric distance is systematically lower than 0.1 dex, compared to at least 0.2 dex in the MW. This appears to be the case for all the FIRE-2 galaxies in Parul et al. (2025) (see their Fig. A1), which likely points to limitations in the subgrid star formation model used in the simulation, such as the use of supernova metal yield prescriptions that do not account for metallicity dependence. Additionally, this systematic lower difference may reduce the number of simulated galaxies identified as having an α-bimodality in Parul et al. (2025). Finally, another striking difference compared to the MW is the relative prominence of the bridge region in Romeo with respect to the analogous counterpart in the MW.

Figure 5 illustrates the spatial evolution of each of the identified disc substructures over time. At z = 0, the low-α disc shapes the outer regions of the galaxy, extending beyond ∼15 kpc, while the bridge dominates the inner regions up to roughly 10 kpc. During the bar epoch, the high-α disc appears to undergo a gradual contraction and its vertical profile deviates from a sech2 profile, consistent with the MW’s qualitative trend (Queiroz et al. 2020). Although it is tempting to relate these spatial changes to the unexpected increase in average rotation of the thick disc at later times and the onset of the bar, further analysis, which is beyond the scope of this paper, would be required to reach this conclusion.

Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Top panel: 2D number density of a sample of giants in Gaia DR3 cross-matched with APOGEE DR17 that spanned a Galactocentric range of 4–16 kpc (see text for details). The red dashed lines delimit the thick-disc, bridge-, low-α regions used to derive the kinematic properties in Table 2, while the excluded region for this derivation is shaded red. Middle panel: 2D number density of star disc particles, as identified by GMM, in the [Mg/Fe] vs [Fe/H] plane. Green, blue and purple mark the contours containing 90% of the high-α, bridge and low-α disc stars, respectively. Bottom panel: mean age distribution of disc stars in the same plane, with the black contour enclosing the region containing 90% of the star particles in the disc. The last two panels show disc star particles at z ∼ 0.1.

3.3 Radial trends

To study the radial trends of the chemical composition of the disc’s components, the stellar disc was restricted vertically to |z| < 1 kpc and binned into concentric annular cylinders. The bin sizes were chosen to facilitate comparison with the radial structure observed in the MW at similar multiples of the disc scale-length, rd (see Fig. 8 in Imig et al. 2023). For the MW, we adopted a scale-length of 2.6 kpc, while for Romeo, a scale-length of 4.5 kpc was determined by fitting the entire disc surface density profile. The top panel of Fig. 6 shows 2D distributions [Mg/Fe] versus [Fe/H] of star disc particles at z = 0.0998 within each ring. However, we also refer to Fig. C.1 for the 2D number density of each GMM disc component.

Overall, the radial structure of the chemical bimodality in Romeo resembles that of the MW (Hayden et al. 2015; Queiroz et al. 2020; Imig et al. 2023), with stars in the low-α disc becoming progressively more metal-poor at larger galactocentric distances. We also find that the ratio of stars with high-α to those with low-α increases with higher |z|. This supports the conclusions of Orkney et al. (2026) that this vertical structure is likely a general result of disc formation in MW-mass galaxies. There are, nonetheless, two very notable differences compared to the MW. Firstly, the high and intermediate-α discs in Romeo remains as dominant as (or even more than) the low-α component out to 3.6×rd (∼16 kpc), unlike the case of the MW. With its more extended stellar disc, Romeo still hosts a notable high-α plus bridge component out to almost 4.7×rd (∼20 kpc), whereas in our Galaxy, the high-α population drops sharply beyond 9 kpc or 3.6×rd. This may suggest that after the spin-up phase and formation of the thick disc, the MW experienced a drop in the SFR, leading to an underpopulated bridge.

Secondly, the change in the location of the [Mg/Fe] versus [Fe/H] plane of the low-α sequence as a function of R is more pronounced in Romeo than in the MW. This second difference could be attributed to subdominant radial stellar migration during the development of the low-α in Romeo (see Appendix C) or it could point to overly large radial fluxes, which could be alleviated by the presence of cosmic ray physics (Trapp et al. 2022) that are not implemented in the current galactic model. Both processes would dilute metallicity differences between radial bins.

The bottom panel of Fig. 6 depicts the age distribution of the star disc particles. By comparing the top and bottom panels of this figure, it can be inferred that there is a deficit of stellar particles with ages around 3–4 Gyr (as can also be seen in Figs. 2 and 4). This deficit driven by a dropping SFR (see bottom panel of Fig. 1) is one of the necessary conditions for forming two separate α-sequences in the absence of radial migration. Following this dip in the number density of stars, the low-α population begins to dominate over the intermediate-α sequence. In our solar neighbourhood, a dearth of stars was found at ages around 6 Gyr (Nissen et al. 2020); however, Fernández-Alvar et al. (2025) instead found a shortage at 3 Gyr. It should be noted that this shortage of stars in the MW is not associated with the onset of the low-α disc, which is estimated to have occurred earlier, around 8–10 Gyr ago (Fernández-Alvar et al. 2025).

Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Time evolution of radial stellar surface density (top) and vertical density (bottom) profiles from z = 1.5 to z = 0 for the disc components identified using GMM. The thickest solid grey line represents the total disc surface density profile at z = 0, while the dashed lines show exponential profiles with scale-lengths rd, obtained by fitting the z = 0 surface density profile of each substructure. The vertical profiles are shown in a cylindrical bin within rd ± 1 kpc (see e.g. Park et al. 2021), where rd is calculated at each look-back time.
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Top: 2D number density of disc star particles within concentric rings with |z| < 1 kpc. The number density distribution is normalised within each ring. Bottom: mean age distribution of disc stars. The red lines depict the chemical evolutionary track of stars formed from in situ gas in each of the rings. In the case of the outer rings, the sequence is discontinuous, as there are periods when no stars form from in situ gas, especially during the spin-up and high-α disc and the bridge-intermediate-α phases. The black contours enclose the region containing 90% of the star particles in the disc.
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Top panel: gas chemical tracks in the inner, central, and outer disc regions. Each track is constructed by, at each snapshot, selecting the gas within the corresponding region and computing its (mass-weighted) average chemical properties. Grey dashed lines indicate schematically the three general trends exhibited by the chemical tracks. Bottom panels: dynamically evolved stellar chemical tracks for the inner (left), central (middle), and outer (right) regions. These are built by selecting, at the final snapshot, the stars located in each region and tracing back their formation times and birth chemistry. Stars are plotted within the track of the region in which they reside at the end, regardless of their place of birth. Each point represents 103 stars, coloured by the time required to form this number of stars; shorter-bluer (longer-redder) formation times correspond to periods of higher (lower) SFR, thereby tracing the SFH of each region.

3.4 Gas and stellar chemical tracks

Figure 7 shows the chemical tracks for the gas (left panel) and for the stars (right panels) in three concentric annular rings. In particular, in an inner, central, and outer ring spanning 0 < R (kpc) < 5, 5 < R (kpc) < 16, and 16 < R (kpc) < 26, respectively. Each cylindrical annuli is restricted vertically to |z| < 2 kpc. While stellar tracks of each annular ring are built by selecting the stars that at the end of the simulation are within that bin, regardless of their birth location, gas tracks are constructed by, at each snapshot, accounting for the gas within a given annular cylinder. The dynamically evolved stellar tracks closely resemble that of the gas. Therefore, stars largely preserve their birth galactocentric distance and stellar migration in Romeo is, on average, negligible or subdominant, as can also be seen in Appendix C. As stars retain the chemical conditions of the cold gas-phase at their birthplace and time of formations, in the absence of substantial radial migration, gas chemical evolution in a given disc region dictates the chemical evolution of the stars.

We can identify three phases in the chemical evolution of the gas and dynamically evolved stellar tracks, based on changes in the slope of these tracks. Interestingly, these three phases are associated with the three disc components identified by GMM and the three phases of the SFH. In particular, we find that:

  • The stars in the thick disc are characterised by high-α abundances with relatively constant values of [Mg/Fe] with increasing [Fe/H] values. These stars form during an early epoch with a rising and bursty SF, when core-collapse supernovae (SNe) exceed SNe type Ia.

  • Next, the star formation begins to decrease and enters a smoother phase. During this period, stars form in the bridge region with intermediate-α values. Their [Mg/Fe] abundance ratios decrease progressively with higher values of [Fe/H], as SNe type Ia begin to outnumber core-collapse SNe, thus dominating the mass yields of Fe mass over the production of Mg. Before the end of bridge formation, around 5 Gyr ago, chemical tracks show a steepening of the slope and even a decrease in [Fe/H] in the outer disc region. As we discuss in Section 4.2, this dilution of the metal content in the disc is caused by the hot, coronal gas.

  • Finally, a mild increase in the overall SFR, together with the metallicity dilution, enables the emergence of a sequence of disc stars with low-α abundance ratios.

It should be noted that in the absence of stellar migration as in Romeo, gas dilution in the disc is a necessary but not sufficient condition for generating the separation and formation of two distinct α sequences. The other necessary condition is either a decrease in the SFR or a rapid and sudden decrease in the [Mg/Fe] abundance ratio of the gas, so that a deficit occurs in the number of stars in the region between the two sequences. In the case of Romeo, the first scenario applies, as can be seen in the lower panels of Fig. 7 by the slight increase in the time elapsed to form 103 stars (decrease in SFR) in the region separating the high and low-α sequences and its subsequent moderate decrease (increase in SFR).

Finally, Fig. 8 shows the chemical tracks of stars in the central bin. In here, tracks are constructed in such a way that the stars are grouped according to whether they formed from the in situ gas or from gas falling into the bin from different directions. It is clear from this plot that radial outflows, primarily driven by stellar winds from massive stars and SNe explosions, carry more metal-rich material outward, producing (on average) more metal-rich tracks. On the other hand, radial inflows bring in lower metallicity gas. Despite some mixing between the infalling and the in situ gas2, stars of the same age at a given radius show significant dispersion in metallicity, reflecting the chemically diverse origins of their progenitor gas. This implies that radial stellar migration is not necessary to explain the chemical spread in a stellar disc and that the commonly adopted assumption of instantaneous mixing in chemical analytical models cannot explain the observed dispersions in metallicity and [Mg/Fe] observed in the Romeo simulated galaxy.

Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Chemical evolutionary tracks of stars, separated according to whether they formed from in situ gas or from infalling gas arriving from different directions.
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Temperature versus number density distribution of the gas in Romeo at tlb = 5.16 Gyr. The horizontal lines mark the cuts used to define the cold, warm, and hot phases, which separate the gas into the cold disc, the cooling bridge connecting the disc with the corona, and the hot corona.

4 Physical drivers of disc evolution

4.1 Gas flow rates

We calculated the gas inflow and outflow rates across the boundaries of the three annular rings defined in Sect. 3.4 and as described in Appendix D. In this way, we evaluate the role of vertical gas accretion and radial gas flows in the development of the three phases of the chemical evolution of the stellar disc. Fluxes are separated by gas temperature, which helps identify their origin. In particular, we distinguish three gas phases: cold (T (K) < 5 × 103), warm (5 × 103 < T (K) < 105), and hot (T (K) > 105) gas. As shown in Fig. 9, these cuts separate the gas into the hot, low-density corona, whose gas gradually cools and shifts into the transition region at intermediate densities, from where it continues to lose energy and slowly accretes onto the cold, dense disc. This transition region also contains originally cold gas from gas-rich mergers that was stripped by ram pressure and heated as its bulk kinetic energy is dissipated into thermal energy.

The above temperature boundaries differ slightly from those adopted by Barbani et al. (2023). The lower cold-warm threshold allows us, on the one hand, to identify as disc gas only the cold component confined to the disc mid-plane; on the other hand, to account for differences between the radiative feedback model used in Barbani et al. (2023) and that adopted in Romeo. In particular, because FIRE-2 does not impose a characteristic temperature for photoionised gas, the slightly lower threshold adopted here ensures that the intermediate-temperature cooling-bridge phase fully encompasses gas at densities ∼10−2 cm−3, whose origin is predominantly the hot corona (see Fig. 9). Additionally, we have verified that adopting either our warm-hot boundary or that of Barbani et al. (2023) does not affect the results of our analysis.

Gas can enter the disc through accretion from the galactic corona, via gas-rich mergers, or through galactic fountains, where stellar feedback ejects gas from the disc into the corona before it eventually falls back (Fraternali & Binney 2008; Marinacci et al. 2019; Barbani et al. 2023). The typical fountain cycle is ∼100 Myr (Fraternali & Binney 2008) and so, outflow and inflow rates are expected to track each other with a delay of this order; however, on average, these rates balance each other out. Such a temporal shift cannot be seen in Fig. 10, as the time span between available snapshots is much larger. Coronal accretion contributes primarily to the warm phase, mergers supply both cold and warm gas, and fountains predominantly recycle cold material.

The upper panels of Fig. 10 showing the vertical net (inflow minus outflow) rates, indicate that, along the disc formation, the vast majority of the gas vertically accreted onto the disc is warm, as can be seen by the overall positive net rate of this gas phase, especially in the central and outer disc regions. This offset is not as substantial as in Barbani et al. (2023), likely reflecting differences in the gas physics model and the absence in Romeo of positive feedback from the efficient mixing between gas ejected from the disc and the hot metal-poor corona, which would otherwise accelerate cooling and boost star formation (Barbani et al. 2023). The accreted warm gas is mainly hot coronal gas that cools down, thus producing an overall positive offset of warm inflows relative to warm outflows. We also arrived to the conclusion that hot coronal gas is the main source of accreted gas by first constructing SUBFIND halo catalogues (Springel et al. 2001) at two look-back times selected just after and before the onset of the thick disc, specifically 9.5 Gyr and 11.6 Gyr ago. We then identified the gas cells that accreted vertically in the disc and located the group or subhalos to which each of these cells belonged in the catalogues. The largest percentage of cells belonged to the ‘outer fuzz,’ which we identified as the intergalactic medium that slowly accretes onto the halo corona and from there into the disc.

Two other trends can be identified from the accretion rates. First, an irregular and accentuated profile of the net vertical rate of warm and cold gas, with high peaks of gas inflow and steep descents where outflows dominate. This initial period, which lasted until approximately 7.7 Gyr ago (corresponding to the beginning of the formation of the bridge or intermediate-α disc), gives rise to a bursty star formation, and is associated with multiple gas-rich mergers that results in the formation of a compact thick disc exhibiting high velocity dispersion like in Brook et al. (2004). This initial chaotic epoch of hierarchical clustering can be seen in the left panel of Fig. 11 that displays the galaxy right before the thick disc starts to settle. At this early time of the galaxy’s collapse, the gas number density traces the primitive filamentary structure of the dark matter Cosmic Web (Joeveer & Einasto 1978; Einasto et al. 1980) and gas is accreted from the intergalactic medium through a filament, with cold gas clouds flowing through it.

After this chaotic epoch, the vertical inflow and outflow rates of the cold phase nearly balance, consistent with recycling-dominated flows, except for a peak around 4–5 Gyr ago, which coincides with a peak in the warm phase. This peak is attributable to minor mergers with gas-rich satellites, as can be seen in the right panel of Fig. 11. Specifically, this panel shows the gas temperature at a look-back time close to the pericentric passages of two minor mergers and clearly reveals cold gas filaments connecting the outer edges of the disc in Romeo with the satellites. The accretion of this gas produces the slight increase in the overall SFR of Romeo, which enables the development of the low-α sequence.

Additionally, the bottom panels of Fig. 10 shows the radial net (inflow minus outflow) rates. From here, we would like to highlight that (in contrast to the vertical accretion) radial accretion is dominated by cold flows in the inner and central disc regions. The inflow of cold, metal-poor gas is mainly a dynamical consequence of low-angular momentum gas accreting vertically onto the rotating disc. Since this gas has a lower specific angular momentum than the disc material, part of it must flow radially inward to conserve angular momentum (Mayor & Vigroux 1981). Thus, although late gas accretion mainly occurs in the outer disc, some of this cold, metal-poor gas provides new material to fuel the star formation of the low-α sequence in the inner regions of the disc.

Thumbnail: Fig. 10 Refer to the following caption and surrounding text. Fig. 10

Vertical (top panels) and radial (bottom panels) net (inflow minus outflow) gas flow rates through an inner (0 < R (kpc) < 5), central (5 < R (kpc) < 16), and outer (16 < R (kpc) < 26) annular rings for cold (T (K) < 5 × 103), warm (5 × 103 < T (K) < 105, middle), and hot (T (K) > 105) gas. Positive net fluxes correspond to inflow (i.e. gas entering the corresponding ring), while negative values indicate gas outflow. The vertical lines in the panels indicate the lookback times (3.6, 7.7 and 10.9 Gyr ago) that separate the formation and evolution of the three disc phases identified in Romeo.
Thumbnail: Fig. 11 Refer to the following caption and surrounding text. Fig. 11

Gas temperature at different lookback times. These are number density-weighted profiles along the perpendicular direction, enhancing the visibility of underdense regions by reducing the dominance of high-density areas. The reference frames are aligned with the principal axes of Romeo at the last simulated snapshot. The yellow circle marks the virial radius, R200,c, at each time.

4.2 Galactic corona

Figure 12 illustrates how the metallicity of the hot corona evolves over time in connection with the accretion history of Romeo. In the top panel, one can observe a metallicity dilution in the outer layers of the corona around 8 Gyr ago. The onset of this dilution in the corona appears to be linked to the infall of a minor merger crossing the virial radius of Romeo (see bottom panel). In fact, in the snapshot just after this infall (tlb = 7.9 Gyr), we separated the coronal gas to whether it belonged right before the infall at tlb = 9.5 Gyr to the ‘outer fuzz’3 other halos or Romeo. The dilution is only found in the gas belonging to other halos. Although this component accounts for only ∼25% of the total coronal gas mass (with the minor merger contributing roughly 15% and the ‘outer fuzz’ contributing about 70%), it is sufficient to drive the observed dilution. At tlb = 11.6 Gyr (z = 3), the minor merger contains slightly less than 1% of Romeo’s stellar mass, about 6% of its gas mass and roughly 6% of its total mass.

As can be seen from the top panel of the same figure, the dilution propagates from the outer to the inner layers of the galactic corona. In the inner regions, the metallicity decline is seen at around 5 Gyr ago, coinciding with the halt (and even reversal) in the [Fe/H] increase in the gas of the disc, as shown by the chemical tracks in the upper panel of Fig. 7. A comparable decline for gas accreted within 10 kpc (similar to 0.1R200,c) was also reported by Parul et al. (2025) (see their Fig. 9), thereby indicating consistency in the different analyses. These results support the interpretation that accreted metal-poor gas from a gas-rich minor merger drives the decline, first in the corona and subsequently in the disc of Romeo.

The bottom panel of this figure also shows the pericentric passages of the gas-rich minor mergers, around 4–5 Gyr ago, responsible for the increase in the inflow rates of the cold and warm gas phases in the central and outer disc regions (see also right panel of Fig. 11). The accretion of gas from these mergers results in a mild increase in the overall SFR and the onset of the low-α sequence. In addition, we can observe that following the pericentric passages, the disc expands (i.e. the accretion of gas from these close approaches facilitates the inside-out growth of the disc). This is consistent with the trends shown in Parul et al. (2025), where the FIRE-2 simulated galaxies with strongest bimodalities exhibit more extended discs.

Thumbnail: Fig. 12 Refer to the following caption and surrounding text. Fig. 12

Top: time evolution of the mass-weighted [Fe/H] of the gaseous corona as a function of lookback time, with linestyles indicating the average metallicity within successive radial cuts (gas within |z| < 2 kpc was excluded to avoid contamination from the disc): r ≤ 0.1R200,c, r ≤ 0.2R200,c, r ≤ 0.5R200,c, and the entire corona. Bottom: density of forming stars as a function of formation time (x-axis) and galactocentric radius r (y-axis), with red indicating higher and blue lower densities. This highlights both the inside-out growth of the disc and the orbits of gas-rich merger events. The black line traces the virial radius R200,c as a function of look-back time.

5 Discussion

5.1 Comparison with previous work

We observe a clear [Mg/Fe]-[Fe/H] dichotomy in the Romeo simulated galaxy, an example of a system with no major mergers and with no bar at present. It should be noted that the high- and low-α sequences are, nonetheless, less separated in [Mg/Fe] than observed in the MW. A similar weak separation is seen in other FIRE-2 galaxies analysed in Parul et al. (2025), which could point to limitations in the sub-grid star formation model, such as the metal yield prescriptions adopted.

We note that while the absence of certain physical processes in Romeo, such as the lack of active galactic nuclei (AGN) feedback, could affect the generality of our findings, our results indicate that neither stellar radial migration nor major mergers are required for a galaxy to develop a disc’s bimodal chemical structure. Instead, we argue that α-bimodalities seem to be a natural consequence of inside-out disc growth in a hierarchical formation scenario. This is consistent with trends shown in Parul et al. (2025), where 11 FIRE-2 simulated galaxies are analysed (including Romeo) and those with a strongest bimodality have the most extended stellar discs.

In the scenario proposed here, the high-α and low-α discs form in two distinct star formation modes, consistent with previous studies (Clarke et al. 2019; Khoperskov et al. 2021). In a similar way to those works, stars on the high-α sequence form, on average, in regions of higher SFR density. Although early star formation is spatially more patchy (whereas at later times it is smoother along the disc), we do not find spatially extended clumps of high SFR density comparable to those reported by Clarke et al. (2019). While Khoperskov et al. (2021) showed that bimodality can arise naturally from inside-out disc growth in isolation, in Romeo minor mergers are required: first to dilute the galactic corona and subsequently to supply fresh gas that fuels the low-α sequence.

5.2 Conclusions

We find that rather than a strict chemical dichotomy of chemically thick versus a thin disc, Romeo exhibits a trichotomy with three main disc assembly phases: high-α, bridge-, and low-α discs, each with distinct formation times, chemical tracks, and kinematic signatures across galactocentric radii. The development of these three disc phases is described below:

  • Similarly to the (semi)isolated galaxy of Brook et al. (2004), the compact, high-α thick disc develops during the spin-up phase, characterised by a turbulent period of hierarchical clustering and a high, bursty star formation rate, resulting in a disc with high velocity dispersion. The spin-up phase, comparable to that of the MW (Belokurov & Kravtsov 2022), is short and lasts ∼2 Gyr. However, in Romeo, it occurs somewhat later and at higher metallicities, between 8 and 11 Gyr ago, compared to 10–13 Gyr ago in the MW. We would like to note that the disc in Romeo is the first to settle among the MW-mass galaxies in the FIRE-2 suite (McCluskey et al. 2024). This aligns with results from other hydrodynamical simulations that tend to form discs later than observed in the MW (Iza et al. 2022; Semenov et al. 2024; McCluskey et al. 2025);

  • Subsequently, a thinner, more extended disc grows as star formation starts decreasing and transitions into a more steady phase, in agreement with early findings by Larson (1976). This disc may correspond to the bridge region identified in the MW (Ciucă et al. 2021), though the prominence of this feature in the MW is significantly weaker than in Romeo. Recent analyses of the MW SFH indicate a quenching phase around 9–10 Gyr ago (Fernández-Alvar et al. 2025), coinciding with the epoch of bridge formation (Ciucă et al. 2021). This would therefore explain the low populated bridge in the MW, which might be a direct consequence of quenching due to the GSE merger event;

  • Finally, around 5 Gyr ago, gas accretes onto the disc due to pericentric passages of gas-rich minor mergers, causing a mild increase in the SFR and favouring the inside-out growth of the stellar disc. Together with the halting of the evolution of the metals due to dilution driven by the galactic corona, this fuels the development of a low-α disc. Although late gas accretion mainly occurs in the outer disc, as expected in an inside-out growth scenario, some of this gas radially flows inwards providing new material to fuel star formation in the inner disc.

In summary, our conclusions are as follows:

  • Stellar discs can develop an α-bimodality as a natural consequence of inside-out disc growth driven by the hierarchical assembly of gas-rich minor mergers. This is not necessary in conflict with observations, which suggest that minor mergers contribute little to the gas accretion budget in the Local Universe (Di Teodoro & Fraternali 2014);

  • As a testable prediction of this scenario, we expect that α-bimodalities are common among MW-mass disc galaxies and do not correlate strongly with either the presence of a bar or the galaxy’s position within the Cosmic Web. Moreover, since the high- and low-α sequences form through slow secular processes, stars in both components are expected to be close to dynamical equilibrium, and the mean of the metallicity distribution function of stars in the low-α sequences at a given radii is constant or only mildly increasing with age. Although a proper statistical analysis, beyond the scope of this paper, is required to confirm these hypotheses;

  • The chemical evolution of the stellar disc cannot be understood without taking into account its interaction with the surrounding galactic corona, which constitutes the main gas reservoir for the disc and is itself replenished mainly through accretion from the intergalactic medium;

  • In Romeo, at the end of the high-α disc formation, a minor gas-rich merger causes a metallicity dilution of the coronal gas. This dilution propagates from the outer to the inner coronal layers and then to the disc, providing one of the two conditions (in the absence of substantial radial migration) necessary for the formation of two distinct α-sequences: the halting (or even the decrease) of metallicity;

  • At the same time, the SFR begins to decrease, generating the second necessary condition for the development of a bimodality: a dip in the number density of stars, followed by gas accretion at the pericentric passages of minor mergers that fuels the onset of the low-α sequence;

  • Finally, we find that instantaneous mixing does not hold in Romeo and the spreads in the stellar [Fe/H] and [Mg/Fe] abundance ratios at a given galactocentric distance reflect the diverse origins of infalling gas, highlighting the importance of accounting for both vertical and radial gas flows.

Looking ahead, determining the ubiquity of α-bimodalities (see e.g. Pinna et al. 2024) and their correlation with galaxy’s position in the Cosmic Web is crucial, as the large-scale environment might regulate the gas content, inflow timing, and merger history. As a result, these effects will influence the disc chemistry and the emergence of multiple chemical sequences (e.g. Khoperskov et al. 2023).

Acknowledgements

We send a special thanks to Claudio Dalla Vecchia for providing the necessary scripts to construct SUBFIND halo catalogues along with useful discussion. We also thank Jenna Samuel and Andrew Wetzel for valuable discussions. MB is funded by the European Union (ERA Fellowship, Dia-LoGues, 101180670). This work was supported by the Estonian Research Council grants (PSG938, PRG1006, PRG3034) the Estonian Ministry of Education and Research (grant TK202) and the European Union’s Horizon Europe research and innovation programmes (EXCOSM, grant No. 101159513, and EXOHOST, grant No. 101079231). GB and EFA acknowledge support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) under grant “EN LA FRONTERA DE LA ARQUEOLOGÍA GALÁCTICA: EVOLUCIÓN DE LA MATERIA LUMINOSA Y OSCURA DE LA VÍA LÁCTEA Y LAS GALAXIAS ENANAS DEL GRUPO LOCAL EN LA ERA DE GAIA. (FOGALERA)”, the European Regional Development Fund (ERDF) with reference PID2023-150319NB-C21 and PID2023-150319NB-C22. EFA also acknowledges support from HORIZON TMA MSCA Postdoctoral Fellowships Project TEMPOS, number 101066193, call HORIZON-MSCA-2021-PF-01, by the European Research Executive Agency. HR acknowledges funding from UK Research and Innovation (UKRI) under the UK government’s Horizon Europe funding guarantee (grant No. 10051045). SCB acknowledges the support of the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (MCIN/AEI/10.13039/501100011033) under grant nos. PID2021-128131NB-I00 and CNS2022-135482 and the European Regional Development Fund (ERDF) ‘A way of making Europe’ and the ‘NextGenerationEU/PRTR’. The results and figures presented in this work were made possible thanks to the following software libraries: Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Jupyter (Kluyver et al. 2016), Scikit-learn (Pedregosa et al. 2011), Gaia-tools (https://github.com/HEP-KBFI/gaia-tools) and GizmoAnalysis (http://ascl.net/2002.015). This last repository was first used in Wetzel et al. (2016).

References

  1. Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826 [NASA ADS] [CrossRef] [Google Scholar]
  2. Andrae, R., Fouesneau, M., Sordo, R., et al. 2023, A&A, 674, A27 [CrossRef] [EDP Sciences] [Google Scholar]
  3. Ansar, S., Pearson, S., Sanderson, R. E., et al. 2025, ApJ, 978, 37 [Google Scholar]
  4. Barbani, F., Pascale, R., Marinacci, F., et al. 2023, MNRAS, 524, 4091 [Google Scholar]
  5. Belokurov, V., & Kravtsov, A. 2022, MNRAS, 514, 689 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527 [CrossRef] [EDP Sciences] [Google Scholar]
  7. Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
  9. Brook, C. B., Kawata, D., Gibson, B. K., & Freeman, K. C. 2004, ApJ, 612, 894 [NASA ADS] [CrossRef] [Google Scholar]
  10. Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765 [Google Scholar]
  11. Ciucă, I., Kawata, D., Miglio, A., Davies, G. R., & Grand, R. J. J. 2021, MNRAS, 503, 2814 [Google Scholar]
  12. Clarke, A. J., Debattista, V. P., Nidever, D. L., et al. 2019, MNRAS, 484, 3476 [NASA ADS] [CrossRef] [Google Scholar]
  13. Di Teodoro, E. M., & Fraternali, F. 2014, A&A, 567, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Dubay, L. O., Johnson, J. A., Johnson, J. W., & Roberts, J. D. 2025, arXiv eprints [arXiv:2508.00988] [Google Scholar]
  15. Einasto, J., Joeveer, M., & Saar, E. 1980, Nature, 283, 47 [Google Scholar]
  16. Fernández-Alvar, E., Kordopatis, G., Hill, V., et al. 2024, A&A, 685, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Fernández-Alvar, E., Ruiz-Lara, T., Gallart, C., et al. 2025, A&A, 704, A258 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Fraternali, F., & Binney, J. J. 2008, MNRAS, 386, 935 [CrossRef] [Google Scholar]
  19. Fuhrmann, K. 1998, A&A, 338, 161 [NASA ADS] [Google Scholar]
  20. Gallart, C., Surot, F., Cassisi, S., et al. 2024, A&A, 687, A168 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Garrison-Kimmel, S., Hopkins, P. F., Wetzel, A., et al. 2019, MNRAS, 487, 1380 [NASA ADS] [CrossRef] [Google Scholar]
  22. Grand, R. J. J., Bustamante, S., Gómez, F. A., et al. 2018, MNRAS, 474, 3629 [NASA ADS] [CrossRef] [Google Scholar]
  23. Guiglion, G., Nepal, S., Chiappini, C., et al. 2024, A&A, 682, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
  25. Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132 [Google Scholar]
  26. Haywood, M., Di Matteo, P., Lehnert, M. D., Katz, D., & Gómez, A. 2013, A&A, 560, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Hopkins, P. F. 2015, MNRAS, 450, 53 [Google Scholar]
  28. Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
  30. Imig, J., Price, C., Holtzman, J. A., et al. 2023, ApJ, 954, 124 [CrossRef] [Google Scholar]
  31. Iza, F. G., Scannapieco, C., Nuza, S. E., et al. 2022, MNRAS, 517, 832 [NASA ADS] [CrossRef] [Google Scholar]
  32. Joeveer, M., & Einasto, J. 1978, in IAU Symposium, 79, Large Scale Structures in the Universe, eds. M. S. Longair, & J. Einasto, 241 [Google Scholar]
  33. Karukes, E. V., Benito, M., Iocco, F., Trotta, R., & Geringer-Sameth, A. 2020, J. Cosmology Astropart. Phys., 2020, 033 [Google Scholar]
  34. Khanna, S., Yu, J., Drimmel, R., et al. 2025, A&A, 701, A270 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Khoperskov, S., Haywood, M., Snaith, O., et al. 2021, MNRAS, 501, 5176 [NASA ADS] [CrossRef] [Google Scholar]
  36. Khoperskov, S., Minchev, I., Libeskind, N., et al. 2023, A&A, 677, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides, & B. Schmidt (IOS Press), 87 [Google Scholar]
  38. Kobayashi, C., Bhattacharya, S., Arnaboldi, M., & Gerhard, O. 2023, ApJ, 956, L14 [Google Scholar]
  39. Larson, R. B. 1976, MNRAS, 176, 31 [NASA ADS] [CrossRef] [Google Scholar]
  40. Lian, J., Zasowski, G., Chen, B., et al. 2024, Nat. Astron., 8, 1302 [Google Scholar]
  41. Mackereth, J. T., Crain, R. A., Schiavon, R. P., et al. 2018, MNRAS, 477, 5072 [NASA ADS] [CrossRef] [Google Scholar]
  42. Marinacci, F., Sales, L. V., Vogelsberger, M., Torrey, P., & Springel, V. 2019, MNRAS, 489, 4233 [NASA ADS] [CrossRef] [Google Scholar]
  43. Mayor, M., & Vigroux, L. 1981, A&A, 98, 1 [NASA ADS] [Google Scholar]
  44. McCluskey, F., Wetzel, A., Loebman, S. R., et al. 2024, MNRAS, 527, 6926 [Google Scholar]
  45. McCluskey, F., Wetzel, A., Loebman, S., & Moreno, J. 2025, arXiv e-prints [arXiv:2506.11840] [Google Scholar]
  46. McMillan, P. J. 2011, MNRAS, 414, 2446 [Google Scholar]
  47. Miller, M. J., & Bregman, J. N. 2015, ApJ, 800, 14 [NASA ADS] [CrossRef] [Google Scholar]
  48. Minchev, I., Chiappini, C., & Martig, M. 2013, A&A, 558, A9 [CrossRef] [EDP Sciences] [Google Scholar]
  49. Nidever, D. L., Gilbert, K., Tollerud, E., et al. 2024, in IAU Symposium, 377, Early Disk-Galaxy Formation from JWST to the Milky Way, eds. F. Tabatabaei, B. Barbuy, & Y.-S. Ting, 115 [Google Scholar]
  50. Nikakhtar, F., Sanderson, R. E., Wetzel, A., et al. 2021, ApJ, 921, 106 [Google Scholar]
  51. Nissen, P. E., Christensen-Dalsgaard, J., Mosumgaard, J. R., et al. 2020, A&A, 640, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Nuza, S. E., Scannapieco, C., Chiappini, C., et al. 2019, MNRAS, 482, 3089 [NASA ADS] [Google Scholar]
  53. Orkney, M. D. A., Laporte, C. F. P., Grand, R. J. J., & Springel, V. 2026, MNRAS, 545, staf1551 [Google Scholar]
  54. Põder, S., Benito, M., Pata, J., et al. 2023, A&A, 676, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Park, M. J., Yi, S. K., Peirani, S., et al. 2021, ApJS, 254, 2 [Google Scholar]
  56. Parul, H., Bailin, J., Wetzel, A., et al. 2023, MNRAS, 520, 1672 [Google Scholar]
  57. Parul, H., Bailin, J., Loebman, S. R., et al. 2025, MNRAS, 537, 1571 [Google Scholar]
  58. Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., 12, 2825 [Google Scholar]
  59. Pinna, F., Grand, R. J. J., Martig, M., & Fragkoudi, F. 2024, A&A, 691, A61 [Google Scholar]
  60. Queiroz, A. B. A., Anders, F., Chiappini, C., et al. 2020, A&A, 638, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Schönrich, R., & Binney, J. 2009, MNRAS, 396, 203 [Google Scholar]
  62. Sellwood, J. A., & Binney, J. J. 2002, MNRAS, 336, 785 [Google Scholar]
  63. Semenov, V. A., Conroy, C., Chandra, V., Hernquist, L., & Nelson, D. 2024, ApJ, 962, 84 [NASA ADS] [CrossRef] [Google Scholar]
  64. Snaith, O. N., Haywood, M., Di Matteo, P., et al. 2014, ApJ, 781, L31 [CrossRef] [Google Scholar]
  65. Spitoni, E., Aguirre Børsen-Koch, V., Verma, K., & Stokholm, A. 2022, A&A, 663, A174 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726 [Google Scholar]
  67. Trapp, C. W., Kereš, D., Chan, T. K., et al. 2022, MNRAS, 509, 4149 [Google Scholar]
  68. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  69. Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJ, 827, L23 [NASA ADS] [CrossRef] [Google Scholar]
  70. Wetzel, A., Hayward, C. C., Sanderson, R. E., et al. 2023, ApJS, 265, 44 [NASA ADS] [CrossRef] [Google Scholar]
  71. Xiang, M., & Rix, H.-W. 2022, Nature, 603, 599 [NASA ADS] [CrossRef] [Google Scholar]
  72. Yu, S., Bullock, J. S., Klein, C., et al. 2021, MNRAS, 505, 889 [NASA ADS] [CrossRef] [Google Scholar]
  73. Yu, S., Bullock, J. S., Gurvich, A. B., et al. 2023, MNRAS, 523, 6220 [NASA ADS] [CrossRef] [Google Scholar]
1

DM mass within a sphere in which the average DM density equals 200 times the critical density of the universe at a given z.

2

Gas tracks when separated by their origin show greater dispersion in metallicity compared to stellar tracks.

3

The ‘outer fuzz’ corresponds to the set of particles not associated with any friends-of-friends halo constituting a diffuse, unbound component surrounding bound structures (Springel et al. 2001). This component is therefore identified as the intergalactic medium.

5

Varying ∆ within reasonable limits does not affect the results.

Appendix A GMM disc components

In Fig. A.1, we present the BIC values as a function of the number of Gaussian components used to model the joint probability density of 3D cylindrical velocities and ages of star particles within R < 40 kpc and |z| < 10 kpc in Romeo. The decrease in BIC values slows noticeably beyond four components. Moreover, increasing the number of components further introduces nonphysical clusters, such as artificial separations in radial velocity into distinct positive and negative groups. Fig. C.1 shows the 2D number density of high-α, bridge-, and low-α stars at different annular rings with |z| < 1 kpc.

Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

BIC as a function of the number of Gaussian components in the GMM analysis of 3D cylindrical velocities and stellar ages within R < 40 kpc and |z| < 10 kpc in the Romeo simulated galaxy.

Appendix B Gaia-APOGEE data preparation and selection criteria

Gaia DR3-provided astrometry and radial velocities, coupled with stellar distance estimates from GSP-Phot (Andrae et al. 2023), were used to transform the sample to the Galactocentric frame and cylindrical coordinates. This was done by adopting the Sun’s orbital parameters from Põder et al. (2023) and using the transformation procedure and the Gaia-tools repository4 described therein. Finally, counter-rotating stars were removed and the sample was radially restricted to a Galactocentric range of R ∈ [4, 16] kpc and vertically within ±2 kpc.

As for the APOGEE DR17 data, following Fernández-Alvar et al. (2024), we selected only stars with ASPCAPFLAG bits 14–42 equal to 0, STARFLAG bits 1–26 equal to 0, and EXTRATARG bits 2 and 4 equal to 0. We also excluded those objects with the apogee1_target1, apogee1_target2, apogee2_target1, and apogee2_target2 classifying them as clusters, streams, dwarf galaxies, sky, telluric, binaries, radial velocity variables, the bar, extended objects, and the TriAnd, GASS, and A13 disc structures. In addition, we select only stars with 3500 < Teff [K] < 6500 and 1 < log g [dex] < 3.5.

Appendix C Radial migration

In Fig. C.1, we show the difference between cylindrical radius of low-α disc stars at formation and at the snapshot prior to bar formation. The majority of the stars lie close to the zero horizontal line. We observe that, although there is an apparent average inward motion for galactocentric radii greater than 25 kpc, the median line overestimates this effect for two reasons. First, there is a statistical bias because the conditioning on the formation radius skews the median toward lower final radii as a consequence of the intrinsic dispersion of the formation radius. Second, even in the absence of radial migration, we expect some stars to move slightly inward and others slightly outward. However, in the outer regions, since the galaxy (or disc) is not infinite, there are no stars that can arrive from beyond the edge to balance the distribution. This lack of symmetry results in a biased distribution, further reinforcing the apparent inward trend.

Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Dierence between the cylindical radius at formation and at the snapshot prior to bar formation versus formation radii for low-α disc stars. The green solid line shows the median trend with the corresponding 1σ scatter shown by the dashed green lines.

Appendix D Calculation of gas flow rates

Following Marinacci et al. (2019), for each snapshot, the flux through a surface was computed as M˙=∑​vnmΔ,Mathematical equation: $\dot M = {{\mathop \sum \nolimits^ {v_n}m} \over \Delta },$(D.1)

where m and vn are the gas cell mass and its velocity component perpendicular to the surface, and the sum runs over gas cells within a slice of thickness ∆ = 0.3 kpc5. In particular, for each annular ring defined by Rmin < R < Rmax and |z| < 2 kpc, vertical fluxes were measured through the top and bottom surfaces by selecting cells within Rmin < R < Rmax and z ∈ [±2 − ∆/2, ±2 + ∆/2]. Cells with vz z > 0 (vz z < 0) were counted as outflows (inflows). Additionally, radial fluxes were measured through the outer cylindrical surface by selecting gas cells with R ∈ [Rmax − ∆/2, Rmax + ∆/2] and |z| < 2 kpc, with cells with vR R > 0 (vR R < 0) classified as outflows (inflows).

At each snapshot, the concentric annular cylinders are aligned with the coordinate system defined by the galaxy’s principal axes, which are calculated as the eigenvectors of the moment-of-inertia tensor. This tensor was determined using the youngest 25% of star particles that, within 10 kpc, comprise 90% of the total stellar mass (Wetzel et al. 2023). While an alternative approach could be to align the coordinate system with the total angular momentum of the stellar disc, as done in Nuza et al. (2019), the current choice of alignment ensures that primarily disc particles are selected. Furthermore, for a disc-like configuration, the principal axis of the moment-of-inertia tensor generally aligns with the total angular momentum vector.

Thumbnail: Fig. A.2 Refer to the following caption and surrounding text. Fig. A.2

Number density for the high-α (top), bridge (middle) and low-α (bottom) discs. The 2D histograms in each row are scaled by the largest amplitude among the six subplots in that row. In this way, the galactocentric distances at which each disc component is most prominent can be seen.

All Tables

Table 1

Global properties of the Romeo simulated galaxy at z ∼ 0.1, just before the start of its bar episode, and at z = 0, compared with those of the MW.

In the text
Table 2

Mean rotational velocity, v¯φMathematical equation: ${{\bar v}_\varphi }$, and velocity dispersions in the cylindrical R and z coordinates for the high-α, bridge-, and low-α sequences.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Top: stellar surface density profiles of the Romeo and Vinter-gatan (Agertz et al. 2021) simulated galaxies, together with the MW profile taken from Lian et al. (2024) and the best-fitting morphology of McMillan (2011). The former MW profile has been normalised to enclose a stellar mass of 4×1010 M⊙. The dot-dashed orange line depicts a falling exponential profile with scale-length of 2.6 kpc. Bottom: star formation rate as a function of time in Romeo calculated with timesteps ∆t = 10 Myr (orange), ∆t = 100 Myr (purple) and ∆t = 500 Myr (blue).
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

1D age distribution of all disc star particles along with the disc components identified by GMM. The vertical grey lines mark key evolutionary stages of the disc identified by GMM (see Sect. 3.2) from the onset of the thick or high-α disc (10.9 Gyr ago) to the time when the bridge or intermediate-α disc begins to dominate over thick disc (7.7 Gyr ago) and the onset of the α−bimodality (3.6 Gyr ago). The vertical grey bands mark the timing of the bar episode.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

2D distribution of the azimuthal or rotational velocity of all stellar particles in Romeo, R < 40 kpc and |z| < 10 kpc, as a function of age (left panel) and metallicity [Fe/H] (right panel), with the median trend shown in red. The green, blue, and purple lines depict the median trend of the stellar population in the high-α, bridge and low-α discs, respectively. In the left panel, the dashed vertical lines indicate the inferred period at which the MW’s stellar disc settles (see Sect. 3.2), while the dashed black line in the right panel shows the inferred trend for the MW as in Belokurov & Kravtsov (2022).
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Top panel: 2D number density of a sample of giants in Gaia DR3 cross-matched with APOGEE DR17 that spanned a Galactocentric range of 4–16 kpc (see text for details). The red dashed lines delimit the thick-disc, bridge-, low-α regions used to derive the kinematic properties in Table 2, while the excluded region for this derivation is shaded red. Middle panel: 2D number density of star disc particles, as identified by GMM, in the [Mg/Fe] vs [Fe/H] plane. Green, blue and purple mark the contours containing 90% of the high-α, bridge and low-α disc stars, respectively. Bottom panel: mean age distribution of disc stars in the same plane, with the black contour enclosing the region containing 90% of the star particles in the disc. The last two panels show disc star particles at z ∼ 0.1.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Time evolution of radial stellar surface density (top) and vertical density (bottom) profiles from z = 1.5 to z = 0 for the disc components identified using GMM. The thickest solid grey line represents the total disc surface density profile at z = 0, while the dashed lines show exponential profiles with scale-lengths rd, obtained by fitting the z = 0 surface density profile of each substructure. The vertical profiles are shown in a cylindrical bin within rd ± 1 kpc (see e.g. Park et al. 2021), where rd is calculated at each look-back time.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Top: 2D number density of disc star particles within concentric rings with |z| < 1 kpc. The number density distribution is normalised within each ring. Bottom: mean age distribution of disc stars. The red lines depict the chemical evolutionary track of stars formed from in situ gas in each of the rings. In the case of the outer rings, the sequence is discontinuous, as there are periods when no stars form from in situ gas, especially during the spin-up and high-α disc and the bridge-intermediate-α phases. The black contours enclose the region containing 90% of the star particles in the disc.
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Top panel: gas chemical tracks in the inner, central, and outer disc regions. Each track is constructed by, at each snapshot, selecting the gas within the corresponding region and computing its (mass-weighted) average chemical properties. Grey dashed lines indicate schematically the three general trends exhibited by the chemical tracks. Bottom panels: dynamically evolved stellar chemical tracks for the inner (left), central (middle), and outer (right) regions. These are built by selecting, at the final snapshot, the stars located in each region and tracing back their formation times and birth chemistry. Stars are plotted within the track of the region in which they reside at the end, regardless of their place of birth. Each point represents 103 stars, coloured by the time required to form this number of stars; shorter-bluer (longer-redder) formation times correspond to periods of higher (lower) SFR, thereby tracing the SFH of each region.
In the text
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Chemical evolutionary tracks of stars, separated according to whether they formed from in situ gas or from infalling gas arriving from different directions.
In the text
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Temperature versus number density distribution of the gas in Romeo at tlb = 5.16 Gyr. The horizontal lines mark the cuts used to define the cold, warm, and hot phases, which separate the gas into the cold disc, the cooling bridge connecting the disc with the corona, and the hot corona.
In the text
Thumbnail: Fig. 10 Refer to the following caption and surrounding text. Fig. 10

Vertical (top panels) and radial (bottom panels) net (inflow minus outflow) gas flow rates through an inner (0 < R (kpc) < 5), central (5 < R (kpc) < 16), and outer (16 < R (kpc) < 26) annular rings for cold (T (K) < 5 × 103), warm (5 × 103 < T (K) < 105, middle), and hot (T (K) > 105) gas. Positive net fluxes correspond to inflow (i.e. gas entering the corresponding ring), while negative values indicate gas outflow. The vertical lines in the panels indicate the lookback times (3.6, 7.7 and 10.9 Gyr ago) that separate the formation and evolution of the three disc phases identified in Romeo.
In the text
Thumbnail: Fig. 11 Refer to the following caption and surrounding text. Fig. 11

Gas temperature at different lookback times. These are number density-weighted profiles along the perpendicular direction, enhancing the visibility of underdense regions by reducing the dominance of high-density areas. The reference frames are aligned with the principal axes of Romeo at the last simulated snapshot. The yellow circle marks the virial radius, R200,c, at each time.
In the text
Thumbnail: Fig. 12 Refer to the following caption and surrounding text. Fig. 12

Top: time evolution of the mass-weighted [Fe/H] of the gaseous corona as a function of lookback time, with linestyles indicating the average metallicity within successive radial cuts (gas within |z| < 2 kpc was excluded to avoid contamination from the disc): r ≤ 0.1R200,c, r ≤ 0.2R200,c, r ≤ 0.5R200,c, and the entire corona. Bottom: density of forming stars as a function of formation time (x-axis) and galactocentric radius r (y-axis), with red indicating higher and blue lower densities. This highlights both the inside-out growth of the disc and the orbits of gas-rich merger events. The black line traces the virial radius R200,c as a function of look-back time.
In the text
Thumbnail: Fig. A.1 Refer to the following caption and surrounding text. Fig. A.1

BIC as a function of the number of Gaussian components in the GMM analysis of 3D cylindrical velocities and stellar ages within R < 40 kpc and |z| < 10 kpc in the Romeo simulated galaxy.
In the text
Thumbnail: Fig. C.1 Refer to the following caption and surrounding text. Fig. C.1

Dierence between the cylindical radius at formation and at the snapshot prior to bar formation versus formation radii for low-α disc stars. The green solid line shows the median trend with the corresponding 1σ scatter shown by the dashed green lines.
In the text
Thumbnail: Fig. A.2 Refer to the following caption and surrounding text. Fig. A.2

Number density for the high-α (top), bridge (middle) and low-α (bottom) discs. The 2D histograms in each row are scaled by the largest amplitude among the six subplots in that row. In this way, the galactocentric distances at which each disc component is most prominent can be seen.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.