VINTERGATAN-GM: The cosmological imprints of early mergers on Milky-Way-mass galaxies

We present a new suite of cosmological zoom-in hydrodynamical ($\approx 20\, \mathrm{pc}$ spatial resolution) simulations of Milky-Way mass galaxies to study how a varying mass ratio for a Gaia-Sausage-Enceladus (GSE) progenitor impacts the $z=0$ chemodynamics of halo stars. Using the genetic modification approach, we create five cosmological histories for a Milky-Way-mass dark matter halo ($M_{200} \approx 10^{12} \, M_\mathrm{\odot}$), incrementally increasing the stellar mass ratio of a $z\approx2$ merger from 1:25 to 1:2, while fixing the galaxy's final dynamical, stellar mass and large-scale environment. We find markedly different morphologies at $z=0$ following this change in early history, with a growing merger resulting in increasingly compact and bulge-dominated galaxies. Despite this structural diversity, all galaxies show a radially-biased population of inner halo stars like the Milky-Way's GSE which, surprisingly, has a similar magnitude, age, $\rm [Fe/H]$ and $\rm [\alpha/Fe]$ distribution whether the $z\approx2$ merger is more minor or major. This arises because a smaller ex-situ population at $z\approx2$ is compensated by a larger population formed in an earlier merger-driven starburst whose contribution to the GES can grow dynamically over time, with both populations strongly overlapping in the $\rm [Fe/H]-\rm [\alpha/Fe]$ plane. Our study demonstrates that multiple high-redshift histories can lead to similar $z=0$ chemodynamical features in the halo, highlighting the need for additional constraints to distinguish them, and the importance of considering the full spectrum of progenitors when interpreting $z=0$ data to reconstruct our Galaxy's past.


INTRODUCTION
The orbits and chemical abundances of stars within a galaxy encode information about its dynamical and enrichment history, providing us with a window into the main events of its cosmological formation history. The advent of the Gaia space telescope has transformed our ability to perform such analysis in the Milky Way, thanks to a dramatic improvement in the quality and volume of astrometric data sets and reconstructed stellar orbital parameters (Gaia Collaboration et al. 2016a,b, 2018. Combined with chemical abundances and radial velocities acquired by large, spectroscopic surveys (e.g. ; ; 3; Majewski et al. 2017;Martell et al. 2017;Conroy et al. 2019, respectively), this now allows us to isolate coherent chemodynamical structures in the solar neighbourhood and link them to ancient events several billion years back in the Milky Way's history (see Helmi 2020 for a review).
Several such clustered structures have now been identified in the high-dimensional space of orbital parameters and stellar abundances ★ E-mail: martin.rey@physics.ox.ac.uk (e.g. Belokurov et al. 2018;Helmi et al. 2018;Myeong et al. 2019;Kruĳssen et al. 2020;Myeong et al. 2022), with the most striking feature being the excess of stars on radial orbits in the local stellar halo around the Sun (the Gaia-Sausage-Enceladus; Belokurov et al. 2018;Helmi et al. 2018; see also Nissen & Schuster 2010;Koppelman et al. 2018;Haywood et al. 2018). This feature is most commonly interpreted as the remnant trace of our Galaxy's last significant merger, with the metallicities, ages and eccentricities of its stars pointing to a dwarf galaxy colliding with the proto-Milky Way around ≈ 2.
However, the exact mass-scale of the GSE progenitor remains to be pinpointed, with stellar mass estimates extending over an order of magnitude ( ★ ∼ 10 8 − 10 9 M ; e.g. Bonaca et al. 2020;Feuillet et al. 2020;Kruĳssen et al. 2020;Mackereth & Bovy 2020;Naidu et al. 2020Naidu et al. , 2021Limberg et al. 2022). This makes it difficult to quantify its mass ratio with the proto-Milky Way and thus its impact on the early Galaxy. It also remains unclear whether the GSE feature observed at = 0 is a pure population that can be robustly linked to one single event, or whether it contains a superposition of multiple population with distinct origins (e.g. Grand Khoperskov et al. 2022a). Furthermore, other chemokinematic debris in the disc and halo could be associated with the merger event (e.g. the 'Splash'; Bonaca et al. 2017;Haywood et al. 2018;Di Matteo et al. 2019;Belokurov et al. 2020;'Arjuna';e.g. Naidu et al. 2020e.g. Naidu et al. , 2021 but could also be of entirely distinct origin (e.g. Amarante et al. 2020;Donlon et al. 2020;Pagnini et al. 2022).
These uncertainties reflect the difficulties of inferring a galaxy's billions of years of dynamical and chemical evolution, from a single data snapshot. Cosmological simulations of galaxy formation provide an ideal framework for such inference, naturally providing an environment within which mergers, mass growth and subsequent star formation and chemical enrichment are self-consistently seeded from the early Universe. However, resolving the internal dynamical structure of galaxies requires large numbers of particles to adequately sample stellar phase-space orbits and avoid spurious heating ( 10 6 ; e.g. Sellwood 2013; Ludlow et al. 2019Ludlow et al. , 2021. Furthermore, resolving the multiphase, dense structure of the interstellar medium (ISM) from which stars form is key to accurately capture their birth kinematics and thus the subsequent dynamical structure and evolution of a disc (e.g. House et al. 2011;Bird et al. 2013).
Recent progress in numerical methods and computing power now allow modern cosmological zoom simulations to meet these requirements within individual Milky-Way-mass galaxies, enabling us to model a handful of objects sampling varying environments and formation scenarios (e.g. Sawala et al. 2016;Grand et al. 2017Grand et al. , 2021Buck et al. 2020;Font et al. 2020;Applebaum et al. 2021;Agertz et al. 2021;Bird et al. 2021;Khoperskov et al. 2022a;Wetzel et al. 2022). With such tools, we can now quantify the frequency of chemodynamical patterns at = 0, link them to specific events in each galaxy's history, and inform the reconstruction of our Galaxy's formation. (e.g. Bignone et al. 2019;Mackereth et al. 2019;Fattahi et al. 2019;Elias et al. 2020;Dillamore et al. 2022;Khoperskov et al. 2022b,c;Pagnini et al. 2022).
However, a causal interpretation of such suites of individual galaxies still remains challenging. The formation scenario of each given galaxy originates from stochastic early-Universe perturbations, making two simulated galaxies' merger histories entirely unrelated to one another. It then becomes difficult to assess how specific = 0 chemodynamical signatures would respond to a change in the early cosmological merger history, and whether such signatures are unique to this formation scenario or can be produced through multiple routes. This fundamentally limits the robustness with which we can reconstruct our Milky Way's past, and our understanding of the degeneracies associated with this inference.
In this work, we address this challenge by combining two approaches: (i) high-resolution (≈ 20 pc), cosmological zoomed simulations using a physical model that can successfully reproduce the chemodynamical structure of Milky-Way-like galaxies Renaud et al. 2021a,b); and (ii) the genetic modification approach (Roth et al. 2016;Rey & Pontzen 2018;Stopyra et al. 2021). Genetic modifications allow us to create different versions of a chosen cosmological galaxy, introducing targeted changes to its formation scenario such as the mass ratio of a specific merger (e.g. Pontzen et al. 2017;Sanchez et al. 2021), the overall formation time (e.g. Rey et al. 2019b;Davies et al. 2021) or the angular momentum accretion (e.g. Cadiou et al. 2022). Each modified version differs minimally from the original scenario, for example conserving the total dynamical mass at = 0 and the large-scale environment around the galaxy. This enables controlled, comparative studies in a fully cosmological context, isolating how a specific aspect of the formation scenario affects the final observables of a galaxy.
Specifically, in this work, we target a Milky-Way-inspired merger history, studying a ≈ 10 12 M dark matter halo which experiences a merger at ≈ 2 on a radially biased orbit similar to that inferred for a potential GSE progenitor. Using the genetic modification approach, we then make the merger mass ratio incrementally smaller and larger to create a suite of five related formation scenarios, all with similar cosmological environment and all converging to the same total dynamical mass to within 10 per cent (Rey & Starkenburg 2022). Each of these cosmological scenarios is then evolved to = 0 using simulations that resolve the galaxy's ISM multiphase structure and include the detailed star formation, stellar feedback and chemical enrichment model used for the project . The five genetically modified galaxies used in this work further form a subset of a larger suite evolved with this model, which will be described in a forthcoming work (Agertz et al. in preperation).
We present the simulation suite in Section 2, and show in Section 3 that modifying the mass ratio of the GSE-like event at ≈ 2 leads to markedly distinct galactic morphologies at = 0, at fixed dynamical and stellar mass. Despite this structural diversity and large variations in mass ratios, we obtain similar GSE-like phase-space features at = 0, whose stars have similar median ages and [Fe/H] distributions (Section 4). We discuss the consequences of our findings on inferring merger mass ratios from Galactic data in Section 5 and summarize in Section 6.

GENETICALLY-MODIFIED MILKY-WAY GALAXIES
We present and analyse a suite of five genetically modified, Milky-Way-mass galaxies, systematically varying the significance of an early ≈ 2 merger similar to the inferred properties of the GSE progenitor. A thorough description of how we construct genetically modified, cosmological zoomed initial conditions to vary the significance of this event is available in Rey & Starkenburg (2022), while the physical model used to evolve them to = 0 is described in Agertz et al. (2021). We summarize these aspects in Section 2.1, and describe the controlled changes to the merger scenarios of each galaxy introduced by our modifications in Section 2.2.

Numerical setup and galaxy formation physics
We construct cosmological, zoomed initial conditions using thesoftware (Stopyra et al. 2021) and cosmological parameters Ω = 0.3139, ℎ = 0.6727, 8 = 0.8440, = 0.9645 (Planck Collaboration et al. 2016). Starting from a dark matter-only cosmological volume with a box size 50 ℎ −1 Mpc ≈ 73 Mpc and mass resolution DM = 1.2 × 10 8 M , we select a target reference halo with Milky-Way virial mass ( 200 ≈ 10 12 M ) and no neighbours more massive within 5 200 , where 200 is the radius enclosing 200 times the critical density of the Universe. We then refine the mass resolution in the Lagrangian region down to DM = 2.0 × 10 5 M , and apply the procedure described in Pontzen et al. (2021) to strongly damp the bulk velocity of the Lagrangian region and limit advection errors during integration. The initial conditions are evolved using linear theory to = 99 (Zel'dovich 1970), before we start following the evolution of dark matter, stars, and gas to = 0 with the adaptive mesh refinement code (Teyssier 2002). We follow the dynamics of collisionless particles (dark matter and stars) using a multiscale particle-mesh solver (Guillet & Teyssier 2011), while fluid dynamics are solved with an HLLC Riemann solver (Toro et al. 1994) assuming an ideal gas equation of state with adiabatic index = 5/3. Our Lagrangian refinement strategy allows us to reach a spatial resolution of 20 pc throughout the galaxy's ISM . We complement this hydrodynamical setup with an extensive galaxy formation model described in detail in Agertz et al. (2021), which we briefly summarize now.
We follow the equilibrium cooling of a metal-enriched plasma (Courty & Alimi 2004;Rosdahl et al. 2013), and model the spatially uniform, time-dependent heating and photoionization from a cosmic ultraviolet background using an updated version of Haardt & Madau (1996) as implemented in the public version. Gas with H ≥ 0.01 cm −3 self-shields from this heating source (Aubert & Teyssier 2010;Rosdahl & Blaizot 2012), allowing it to condense to densities ≥ 100 p cm −3 and temperatures ≤ 100 K at which we model star formation through a Poisson process following a Schmidt law (Schmidt 1959;Rasera & Teyssier 2006;Agertz et al. 2013). Newborn stellar particles are sampled with 10 4 M masses and are treated as single stellar populations with a Chabrier (2003) initial mass function.
We track the age-dependent injection of momentum, energy, and metals from stellar winds in O, B, and asymptotic giant branch (AGB) stars, and explosions of Type II and Type Ia supernovae (SNe) according to the budget defined in Agertz et al. (2021) (see also Agertz et al. 2013Agertz et al. , 2020Agertz & Kravtsov 2015). Feedback from SNe is injected as thermal energy if the cooling radius is resolved by at least 6 grid cells, and as momentum otherwise (Kim & Ostriker 2015;Martizzi et al. 2015;Agertz et al. 2020Agertz et al. , 2021. We track the evolution of two metals, iron and oxygen, using progenitor-mass-dependent stellar yields for SNe (Woosley & Heger 2007), a delay-time distribution for SNeIa ('field' in Maoz et al. 2012), and a continuous slow, release from AGB stars (Agertz & Kravtsov 2015).
We identify dark matter haloes and subhaloes using the AHF structure finder (Gill et al. 2004;Knollmann & Knebe 2009), retaining only structures with more than 100 particles. We construct merger trees by matching haloes and subhaloes across each simulation snapshot using the (Pontzen et al. 2013) and (Pontzen & Tremmel 2018) libraries. Halo centres are identified using the shrinking-sphere algorithm (Power et al. 2003), and we shift the coordinate frame to ensure that velocities within the central kpc vanish. We define the total stellar mass of each galaxy, ★ , by summing the stellar mass within 200 , and we interpolate a single stellar population model (Girardi et al. 2010) over a grid of ages and metallicities to obtain the luminosities of individual stellar particles. The projected half-light radii, 1/2, is then derived along a random line of sight. We define the iron and abundance ratios following Agertz et al. (2021), equation (3) (see also Escala et al. 2018) converting to solar ratios using Asplund et al. (2009), table 1, and using oxygen as an approximation for the total abundance as it dominates the mass fraction. We compute total metallicities from these two elements following Kim et al. (2014), equation (4).

Genetic modifications and merger scenarios
The genetic modification technique makes targeted changes to a galaxy's cosmological initial conditions to modify its later non-linear merger history in a controlled way. Modifications and initial conditions used in this work are extensively described in Rey & Starkenburg (2022) in a dark-matter-only context (see their 'Milky-Way-like' family) -here, we re-evolve this same family of initial conditions with hydrodynamical galaxy-formation simulations (Section 2.1) and now present their resulting merger scenarios.
Our study targets a dark matter halo with a total dynamical mass 200 ≈ 1×10 12 M and which experiences an early, large interaction with a radially biased orbit (infall at ≈ 2 with a 200 ratio of 1:6 and a radial-to-tangential velocity ratio between the progenitors / = and stellar (bottom) mass assemblies over cosmic time across the suite of genetically modified galaxies. The significance of a merger at ≈ 2 (B, red dot) in the reference case (blue line) is incrementally decreased (navy and purple lines) and increased (cyan and green) using genetic modifications. By construction, the final dynamical mass of all galaxies is fixed, inducing small compensating shifts in an early major event (A, orange) and a late minor merger (C, yellow, see Table 1 for all merger properties). All galaxies converge to similar stellar masses at = 0 (bottom) and evolve in similar large-scale environment ( Figure 2). 16 at this time; see 'Target Progenitor' in Table 1). We choose this history to broadly resemble the inferred properties of the progenitors of the proto-Milky Way (Belokurov et al. 2018;Helmi et al. 2018).
Using genetic modifications, we then aim to increase and decrease the significance of this event, bracketing the range of reported progenitor masses and mass ratios. We thus identify the Lagrangian patch of this early merger in the reference object and define linear modifications to increase or decrease its mean overdensity and control the merger mass ratio (see Rey & Starkenburg 2022, table 1 for the detailed modifications). Figure 1 shows the assembly of dynamical (top) and stellar (bottom) mass in the main progenitor of each merger scenario. Our modifications successfully conserve the early ( ≥ 5) and late ( ≤ 0.5) accretion histories, converging by design to similar dynamical masses at = 0 to within 10 per cent (Rey & Starkenburg 2022). Further, Figure 2 shows that the large-scale environment around the galaxy is largely unchanged between scenarios. This is a natural consequence of the genetic modification approach which aims to make minimal changes to variables untargeted by the modifications while retaining consistency with the ΛCDM cosmology (see also Pontzen   geometry is maintained in each case. Due to correlations and compensations to conserve the total mass, other events ('A' and 'C', orange and yellow, respectively) see their infall time and merger ratios slightly modified (see Table 1 for a complete quantification). Mass growth histories around ≈ 2 however diverge, as expected following our explicit targeting of a merger event at this time. Before quantifying these changes, we briefly assess the likelihood of each genetically modified accretion history to verify any potential rarity in a ΛCDM universe. The relative likelihood of each new modified initial condition to the reference, Δ 2 , remain small compared to the available number of degrees of freedom -Δ 2 = −3.6, −2.77, −0.45, +0.3 for increasing merger mass ratio compared to ≈ 10 6 modes in the zoom region -ensuring their compatibility with the ΛCDM power spectrum. We further compute the median and 68-95% fractional mass growth histories across a population of 28,475 dark matter haloes extracted from the IllustrisTNG simulation (Nelson et al. 2019, grey contours in Figure 1, see Rey & Starkenburg 2022 for further details on the computation). All our mass growth histories are within the 1 contour of the overall population around ≈ 2, demonstrating that the scenarios studied in this work are all plausible cosmological realizations (see also Section 5.2 for further discussion on compatibility with ΛCDM merger rates).
To quantify differences in merger histories, we extract merger events in the reference scenario that (i) bring at least ★ ≥ 10 8 M of accreted material at infall and with (ii) merger mass ratios more significant than 1:30 in 200 . These cuts ensure that we eliminate both very high-redshift events with high mass ratios but low significance to the overall content, and late, low-redshift mergers with small mass ratios. This flags three key progenitors in the reference merger history (labelled 'A', 'B'. and 'C' in Figures 1 and 2), which, as we will see in Section 4 all play a role in defining the = 0 chemokinematic structure within our galaxies.
We cross-match these merger events across our genetically modified simulations, and report in Table 1 their infall redshifts, defined at the time at which the infaller's centre is last outside the main progenitor's virial sphere, and their masses and mass ratios at this time. We further compute the ratio between the radial and tangential velocities of the two halo centres at this time to assess the angular momentum of the encounter, and identify the times of first pericentric passages, pericentre , and end of the interaction, coalescence , using the high-cadence simulation movies.
Our genetic modifications explicitly target merger 'B' (second row in Table 1; red in Figures 1 and 2), which in the reference case is the last major event until = 0. We make it incrementally more significant in each scenario in both mass ratio, dynamical and stellar mass (Table 1), while retaining a strongly radial approach in all cases ( / ≥ 5). The range of merging stellar masses scanned by our suite (2 × 10 8 M ≤ ★ ≤ 2 × 10 9 M ) accurately brackets the proposed masses for a potential GSE progenitor (e.g. Belokurov et al. 2018;Helmi et al. 2018;Feuillet et al. 2020;Kruĳssen et al. 2020;Mackereth & Bovy 2020;Naidu et al. 2020Naidu et al. , 2021Limberg et al. 2022).
Performing these targeted changes to merger 'B' however leads to other alterations to the merger history, due to the correlations inherent to a cosmological context (Roth et al. 2016;Rey & Pontzen 2018). For example, reducing the mass of an event (e.g. 'B') while maintaining the same total mass at = 0 is compensated by growing other events (here mainly 'A', Table 1). Further genetic modifications could force the mass ratio of the earlier event 'A' to match across all scenarios and control for this effect. We leave such finer control to future work, and stress that having multiple significant, high-redshift progenitors is a generic prediction of ΛCDM (see Section 5.2) and is thus inherently reflected by the cosmological nature of the genetic modification approach.
Another effect of our modifications is to slightly alter the infall times of each merger as their significance and mass varies (Table 1). This is most visible in Figure 2, where merger 'A' (orange) is already overlapping with the main progenitor (white) in the rightmost panel, but is further away in the leftmost case. Such shifts in infall time arise from correlations between the linear density, velocity, and potential fields in the initial conditions -as we increase or decrease the local overdensity to modify the merger mass ratio, it smooths or sharpens the local potential gradient towards the main progenitor, in turn modifying the velocity field and introducing shifts in merger timings (see Pontzen et al. 2017;Rey et al. 2019a for further examples). Again, infall times could be fixed to their reference values using additional modifications targeting the velocity or angular momentum structure of the Lagrangian patch (e.g. Cadiou et al. 2021;Pontzen et al. 2021). However, the differences in timings (≈ 300 Myr) remain small compared to uncertainties in dating such early merger events in our Milky Way (e.g. Bonaca et al. 2020;Feuillet et al. 2021), and we thus leave a more detailed setup simultaneously controlling infall times and merger ratios to a future study.
By design of our modifications, all host dark matter haloes converge to the same dynamical mass at = 0, but we note that their central galaxies also match in their final stellar mass ( ★ = 1.8×10 10 M within 20 per cent of each other; Figure 1). All final stellar masses are compatible at 1 with empirical model predictions for this halo mass (Moster et al. 2018;Behroozi et al. 2019), and we will show in a forthcoming work that their growth over time is also compatible with such models. We interpret this convergence in stellar mass as a byproduct of the extremely similar dynamical mass assemblies and large-scale environments of each genetically modified galaxy.
To summarize, our five genetically modified scenarios provide us with a controlled study that systematically varies the significance of an early, radially infalling merger at ≈ 2, while fixing the largescale cosmological environment, and the total dynamical and stellar mass budget of a Milky-Way-mass galaxy.

RESPONSE OF THE GALAXIES' STRUCTURE TO GROWING AN EARLY MERGER
We start by quantifying the response of the = 0 stellar structure of each galaxy as we genetically modify their early assembly history. Figure 3 shows UVI-mock images of the stellar light at = 0 (top row) and luminosity-weighted line-of-sight velocity maps (bottom row) as we incrementally increase the significance of the ≈ 2 merger (left to right). Images do not account for dust attenuation, span a surface brightness interval from 20 to 30 mag arcsec −2 and are viewed along the simulation's z-axis (i.e. a random line of sight physically, but a consistent orientation across panels).
Starting from the reference case (central panels), the galaxy showcases an irregular morphology at = 0, with a rotationally supported, inner galactic disc which extends into a bluer corotating outer structure misaligned with the inner stars. As we make the early history of the galaxy quieter (from central to left-hand panels), the disc orders into a single kinematic component, and spatially grows. By contrast, a more significant merger increases the compactness of the galaxy, forming increasingly brighter central bulges and reducing the overall rotational support (fourth and fifth panels).
The galaxy with the largest merger (right-most panels) exhibits an inner core, counterrotating compared to an outer, low surface brightness (≥ 27 mag arcsec −2 ) star-forming disc. Such structure is reminiscent of the decoupled internal kinematics observed around local elliptical galaxies (e.g. Efstathiou et al. 1982;Bender 1988;Krajnović et al. 2011Krajnović et al. , 2015Johnston et al. 2018;Prichard et al. 2019), for which such low-surface brightness discs have been revealed by deep imaging (e.g. Duc et al. 2015). This galaxy stands out in our suite compared to the more ordered, rotationally supported objects (first four columns), and is difficult to relate to studies of our Galaxy. We present it here as a useful complement to establish trends with growing mass ratio, and will provide a dedicated study of its internal kinematics in a follow-up work.
Our results highlight a broad trend: quieter histories favour larger stellar discs (left), while early major mergers produce more compact, bulge-dominated morphologies at = 0 (right and Section 4). This trend aligns with expectations that mergers drive, in part, morphological transformations of galaxies (e.g. Negroponte & White 1983;Di Matteo et al. 2007;Hopkins et al. 2009;Naab et al. 2014;Zolotov et al. 2015;Pontzen et al. 2017;Martin et al. 2018;Davies et al. 2021). We cleanly isolate this effect here, as our transformations in galactic structure occur at fixed dynamical masses, fixed galaxy stellar masses, and similar large-scale cosmological environments (Figures 1 and 2).
However, our different histories at fixed mass showcase further diversity in galactic morphological and kinematic structure than this broad trend. The effective radii (labelled in each panel) are modified by a factor five from left to right, and do not respond linearly to the growth of the early ≈ 2 merger (e.g. second, third, and fourth columns decrease and increase in turn). We verified that this  (top) and kinematics (bottom) to modifying our galaxies' early histories. As we grow the ≈ 2 merger (left to right), the galaxy transitions from an extended, rotationally supported disc (left-most), to a more compact, bulge-dominated structure (right-most). All changes occur at fixed galaxy stellar and dynamical mass ( Figure 1) and in similar large-scale environment ( Figure 2), showcasing the diversity in galaxy morphology and kinematics introduced by different cosmological histories. Despite this diversity, all scenarios forming a well-ordered galactic disc (first four columns) showcase a GSE-like phase-space structure ( Figure 4, 5), that has distinct cosmological origin but overlapping chemodynamics and ages in every case ( Figure 6, 7).
is also the case for 3D stellar half-mass radii, which are less sensitive to mass-to-light assumptions. Moreover, each galaxy's disc significantly varies in orientation with respect to the same frame of reference in Figure 3. This is to be expected as other physical variables that affect disc formation are evolving across our genetically modified scenarios. In particular, the spin-orbit coupling, the impact parameter, and the resulting tidal fields in the ≈ 2 interaction vary between scenarios and can play a key role in setting the bulk of the angular momentum of the future galaxy (e.g. Hernquist 1993;Springel & Hernquist 2005;Hopkins et al. 2009;Renaud et al. 2009Renaud et al. , 2021a. Furthermore, later aspects of the galaxy's history ( ≤ 1) relevant to disc formation are also modified, notably the constructiveness of filamentary gas accretion (e.g. Dekel et al. 2020;Kretschmer et al. 2020), the structure of the inner CGM and its angular momentum content (e.g. Stern et al. 2021;Hafen et al. 2022;Gurvich et al. 2023), or the orbital parameters and spin of late minor mergers (e.g. Jackson et al. 2022).
None the less, despite this structural diversity in the inner galaxy, we will see in Section 4 that the chemodynamics of halo stars in each galaxy remain largely similar. We hypothesize that the misaligned discs in the reference scenario highlight a 'turning point' in cosmological angular momentum accretion, with these two discs getting constructively aligned when diminishing the early merger to produce larger discs, and destructively counteraligned when growing the early merger, leading to more compact morphologies (e.g. Kretschmer et al. 2020). We leave a confirmation of this scenario to a future study, and now focus on extracting the signatures of each merger scenario in phase-space.

Phase-space structure and radially biased halo stars
To construct the chemodynamical structure of each galaxy, we start by extracting the phase-space distribution of all stars within the inner 50 kpc. Figure 4 shows the fractional distributions of their specific vertical angular momentum, z , and orbital energies, E tot . We define the vertical disc direction from the parallel to the stellar angular momentum in the inner 5 kpc, and compute E tot by summing the kinetic energy of each stellar particle with the local gravitational potential computed by the simulations' Poisson solver. This accounts for asphericity, time-dependence, and long-range cosmological fluctuations in the potential -we thus re-normalize it to have vanishing E tot at the virial radius and ease comparison with isolated dynamical studies. The total number of stars normalizing each histogram is shown at the top of each panel of Figure 4, and changes by less than 20 per cent across scenarios. This follows the minimal changes between each galaxy's total stellar mass (Section 2.2) and ensures that comparisons between fractional phase-space distributions are minimally affected by different normalizations.
Starting from the 'smallest = 2 merger' scenario (left-hand panel), we can first identify different galactic components in phasespace: (i) the bulge, with bound orbits (low E tot ) and pressuresupported kinematics (low z ); (ii) the stellar disc, as a sequence spreading across a range of E tot and extending towards rotationally supported, highz orbits; and (iii) the fainter stellar halo surrounding the galaxy, towards higher E tot and small z . The clustered, coherent structures in the stellar halo's phase space map onto disrupted debris on long dynamical time-scales and surviving dwarf galaxy satellites around each galaxy.
Each of these generic components can be mapped across the different formation scenarios, recovering the visual and kinematic structural trend highlighted in Figure 3. The disc sequence in phase space is the most extended in z for the quietest merger history (left-hand panel), significantly smaller for intermediate scenarios (second, third, and fourth columns) and lacking for the largest merger (right-hand panel). Furthermore, we can identify in phase-space the specifics kinematic structures noted previously. The misaligned stellar discs in the reference scenario (central panel) appear as a double sequence, with an inner, z -tail linking to the bulge, and a population overlapping in z but at higher E tot . Similarly, the decoupled kinematics for . Fractional distributions of specific energy and vertical angular momentum for all stars within 50 kpc at = 0, as we increase the significance of the early merger (left to right). Broad trends and specific galactic structures identified in Figure 3 are recovered in phase space (annotated), in particular the increasingly smaller stellar discs with increasing merger ratio (e.g. lack of a highz sequence in the right-most panel). A population of inner halo stars (low z , E tot comparable to disc stars) is present in all merger scenarios, but becomes more prominent compared to the highz disc population with increasing merger mass ratio (fourth panel). This population of inner halo stars is strongly radially biased, resembling the GSE population in the Milky Way ( Figure 5). Clustered structures in the outer halo (e.g. left-hand panel) map to dwarf galaxy satellites and recently disrupted debris. well-ordered disc galaxies (first four columns), the velocity distribution is bimodal, with disc stars (positive , orange boxes) linking into a radially biased population ( / 1, purple box) alike the Milky-Way's GSE structure. We detect this latter population in all cases, in similar absolute numbers (purple), despite the diversity in final galactic structure and past merger mass ratio at ≈ 2. Further to these kinematic similarities, radially biased GSE-like stars exhibit similar age and [Fe/H] distributions in all cases (Figures 6 and 7). the largest merger (right-hand panel) are visible as a minimal prograde disc sequence for the inner stars (small E tot ), and retrograde, higher E tot orbits for the outer stars.
Beyond these already noted features, the phase-space distributions exhibit a significant population of inner halo stars, with pressuresupported orbits (low z ) that share the same range of E tot as disc stars. This population is visible in all galaxies but becomes particularly prominent for the larger merger mass ratios (e.g. fourth column), as expected if a single early event dominates the assembly of the inner halo. In fact, Figure 4 exhibits a systematic trend: increasing the mass ratio of an early merger makes the inner halo population more prominent relative to disc stars (darker colours compared to the disc stars from left to right).
To quantify this trend in the stellar halo in more detail, we extract the kinematics of the stellar population, excluding the innermost, bulge-dominated stars (i.e. stars with 2 ≤ ≤ 50 kpc) and show their azimuthal and radial velocity fractional distributions in Figure 5. All merger scenarios with a well-ordered disc component (first four columns; the right-most galaxy showcases counterrotating inner and outer kinematics, hence its misalignement in this plot) exhibit a bimodal kinematic distribution, with (i) a population of stars with positive azimuthal velocities rotating with the disc (highlighted by the orange boxes), which links into (ii) a population of stars with small azimuthal velocities over a large radial extent (purple boxes). This latter population of radially biased stars is highly reminiscent of the GSE population in the Milky Way (e.g. Belokurov et al. 2018;Helmi et al. 2018), which is key to claims that our Galaxy underwent a large merger around ≈ 2.
However, we recover this orbital structure in each of our galaxies with ordered rotation (first fourth columns), whether our merger at ≈ 2 is enhanced to become a major event (stellar mass ratio up to 1:4; fourth panel) or decreased to a minor event (down to a merger ratio less than 1:20; left-most). We further extract the number of stars on radial orbits by integrating the distribution within the purple boxes and quote it in each panel. Despite their diversity in early merger history and final galactic structure, all four galaxies have nearly identical absolute number of stars on radial orbits, changing by at most thirty per cent across scenarios. This is comparable to changes in overall stellar mass (numbers in Figure 4) and to variations when observing the galaxy at different nearby timestamps (Section 4.3).
We can further quantify the relative contribution of inner halo stars compared to disc stars defined within the orange boxes in Figure 5 (orange number). The fraction of radial-to-disc stars follows a systematic trend with ≈ 2 merger mass ratio (left to right), going from 0.34 to 0.37 to 0.43 to 0.59 for the four disc galaxies, and to 2.9 in the last scenario which lacks well-ordered rotation. We checked that this comparative trend holds when selecting stars with ≥ 1 kpc and ≥ 3 kpc, or focusing only on halo-like orbits with | z | ≤ 0.5 kpc km s −1 .
Our results thus demonstrate that GSE-like inner halo populations can be assembled through both minor and major events at ≈ 2, and might be relatively common across the Milky-Way analogue popula- tion (see Section 5.2 for further discussion). Further, the overall mass of a GSE-like population does not directly link to the mass ratio of a single, early merger. Rather, a larger merger systematically increases the relative contrast between this population and the galactic disc in phase-space, at fixed galaxy stellar mass and dynamical mass. We thus caution that inferring the merger mass ratio and mass-scale of the GSE progenitor requires quantifying relative weights between subpopulations (e.g. disc to inner halo) which, in turn, calls for a thorough understanding of the completeness and selection functions of both the observed and simulated data (see Section 5.1 for further discussion). This weak connection between GSE-like structures at = 0 and the mass ratio of a ≈ 2 merger questions the cosmological origin of these stars, and we now rewind each progenitor's history to establish their exact nature.

Ages and metallicity distributions
We now focus on testing the cosmological origin of the GSElike stars identified kinematically in Figure 5 (purple box with | | ≤ 50 km s −1 , | | ≤ 240 km s −1 , and 2 ≤ ≤ 50 kpc). We start by plotting their star formation histories in Figure 6 (merger significance increasing from left to right). This does not account for stellar mass-loss across cosmic time, but we checked that our results are unchanged if plotting the mass-weighted distribution of stellar ages.
In all cases, ages show a broad distribution, with most stars being older than 9 Gyr ( ≥ 1.7), as expected if they track ancient stellar populations associated to early merger events. However, despite the order-of-magnitude variation in the merger mass ratio at ≈ 2, all GSE-like populations have similar median ages -11.2, 11.3, 11.5, 11.3, 11.4 Gyr for each panel, respectively -which would be hard to distinguish from data once convolved with realistic uncertainties (≈ 1 Gyr; e.g. Feuillet et al. 2016;Bonaca et al. 2020;Xiang & Rix 2022).
A tail of later forming stars is also visible in all cases ( birth ≥ 6 Gyr), reflecting a contamination from younger stars with overlapping kinematics at = 0. We will confirm in Section 4.2.2 and Appendix A that these are young, [Fe/H]-rich, [ /Fe]-poor stars consistent with being formed within the disc and in satellite dwarf galaxies that infall later into the main progenitor. Such contamination could likely be minimized through optimized chemical and dynamical selections tailored to each galaxy, but we rather choose broader cuts to provide consistent and complete comparisons between scenarios.
Furthermore, we show in Figure  58 from left to right panels, respectively. We thus conclude that the order-of-magnitude scan in merger mass ratio at ≈ 2 not only produces kinematic similarities ( Figure 5), but also similar distributions of ages and metallicities for GSE-like stars.
To understand the nature of these similarities between merger scenarios, we note that the star formation histories of the radially biased halo population at = 0 also exhibit singled-out enhancements compared to the background distribution. We mark at the top of Figure 6 the length of the interaction with mergers 'A' and 'B' (see Table 1 for infall and coalescence 1 ). For all regular disc galaxies (first four panels), peaks of star formation in the GSE-like population coincide with interactions, and we show in Appendix B that such peaks further map to galaxy-wide starbursts in the overall stellar population. Such interaction-driven starbursts have been recently detected in the star formation history of the Milky Way, where peaks of star formation coincide with probable pericentre passages of the Sagittarius dwarf galaxy (e.g. Ruiz-Lara et al. 2020). Starbursts triggered by early, gas-rich mergers thus play a significant role in assembling the GSElike population at = 0, indicating that a significant fraction of its stars are of in-situ origin (see also Grand et al. 2020;Orkney et al. 2022;Khoperskov et al. 2022a) and could elevate the galaxy-wide star formation rate ≈ 10 Gyr ago (e.g. Alzate et al. 2021).
Furthermore, the long tails towards low ([Fe/H] ≤ −1.0) and high ([Fe/H] ≥ −0.5) metallicities, and the spread in ages, imply that multiple populations of different origins are populating the = 0 GSE-like population. In fact, despite having the most massive merger at ≈ 2, the last scenario (right-most) lacks a correspondence between the interaction of 'B' and a peak of star formation. We show in Appendix B that this merger does in fact trigger a large, galaxy- wide, star formation enhancement, but that stars formed during this burst do not populate GSE-like orbits at = 0. The nature of radially biased halo stars and the respective importance of each interaction is thus evolving across the different merger scenarios, which we quantify now.

Multiple co-exisiting in-situ and ex-situ populations
To establish the respective origin of stars and the role of each merger event in assembling the = 0 GSE-like feature, we track each GSElike star through the cosmological merger trees and assign them to eight mutually exclusive categories: In-situ pre-A: Stars that are within the main progenitor before the first infall of merger 'A' (Table 1). We define a star's membership to a dark matter halo from the halo finder, and verified that defining it from spheres of 200 and 2 1/2, around the halo centre does not affect the relative trends reported in this work. Accreted A: Stars that are within 'A' at its first infall, which form a clean sample of ex-situ stars from this event.
Interaction A: Stars that are within the main progenitor between the first infall of merger 'A' and its coalescence . Stars in this subpopulation primarily form during the starburst driven by this interaction (recall the peak in star-formation in Figure 6) and can originate from gas within the main host, from gas within 'A' or from the gas tail produced by the interaction. We avoid distinguishing between in-situ and exsitu origins for this population, as this distinction becomes ever more difficult as the two merging bodies coalescence (although see e.g. Cooper et al. 2015).

In-situ A-B:
Stars that are within the main progenitor between the end of merger 'A' and the first infall of 'B'. These are primarily in-situ stars, although a small, sub-dominant population of ex-situ stars deposited by very minor mergers between 'A' and 'B' also contributes. Accreted B: Stars that are within 'B' at its first infall. Interaction B: Stars within the main progenitor between the first infall of merger 'B' and its coalescence . Accreted C: Stars within 'C' at its first infall.
Other Stars which have not been assigned to the previous categories. As we will see, these are primarily young, [Fe/H]-rich, low-[ /Fe] stars formed in the disc, or deposited later in the halo from other less significant mergers than 'A', 'B' and 'C' (see also Appendix A). In all merger scenarios, a varying mixture of the subpopulations identified above contributes to the metallicity distributions, each peaking at distinct [Fe/H] but always with broad and overlapping spread. The 'Interaction A' and 'Interaction B' (blue and green, middle and bottom row, respectively) subpopulations are particularly prominent, confirming that merger-driven starbursts produce a significant fraction of the total = 0 GSE-like population.
As we increase the importance of merger 'B' (left to right), its accreted population (red, middle row) rises. It only dominates, however, in the most extreme scenario (right-most), for which the galaxy transitions into a bulge-dominated morphology (Figure 3). In fact, when the final galaxy showcases organized rotation at = 0 (first four columns), the population of stars formed during the earlier interaction 'A' (blue, middle row) dominates the overall budget of radial stars for small mass ratios (1:24 and 1:15; first two columns), and has nearly perfectly overlapping [Fe/H] distribution with the accreted population (red) for medium mass ratios (1:8 and 1:4; third and fourth columns). Our results thus reaffirm that multiple components associated to several early accretion events can mix and contribute to a galaxy-wide, metal-poor radial population (see also Grand et al. 2020;Donlon et al. 2022;Donlon & Newberg 2022;Orkney et al. 2022).
Further, within each individual merger scenario, distinguishing such subpopulations chemically could be highly ambiguous. One might expect a population of stars formed in the central regions during the interaction to be more metal-rich than the accreted component associated to the merger body. We recover this trend for each galaxy, where the accreted component of a given interaction peaks toward lower [Fe/H] (orange versus blue; red versus green). However, telling apart an accreted component at ≈ 2 from an earlier in-situ population is difficult, as they peak around the same metallicity and exhibit extended distributions that overlap with one another (blue versus red; see also Renaud et al. 2021a). We further show in Appendix A that, similarly to [Fe/H], they also overlap in E tot and [ /Fe] and their 2D combinations.
Furthermore, these ambiguities strongly propagate to potential inferences of the merger mass ratio at ≈ 2. For the four disc galaxies (first four columns), diminishing the mass ratio at fixed total stellar and dynamical mass (recall Figure 1)  Our results thus highlight strong degeneracies when inferring the early merger history of a galaxy from its chemodynamics, with multiple histories leaving overlapping signatures at = 0. We stress that there likely remains powerful avenues to distinguish each merger scenario, in particular leveraging additional chemical information and precise stellar ages (e.g. Myeong et al. 2019;Bonaca et al. 2020;Feuillet et al. 2021;Matsuno et al. 2022). None the less, we show that a robust inference of our Galaxy's past require a wider scan through merger histories than previously explored, which accounts for the full spectrum of high-redshift progenitors and the response of the protogalaxy to each of them, combined with a careful account of their observables in the solar neighbourhood. We discuss this latter aspect further in Section 5.1 and focus next on quantifying when the GSE-like structure and its subpopulations are assembled across cosmic time.

Are stars born, deposited, or kicked onto radial orbits?
Our results establish that in-situ populations forming in earlier interactions overlap kinematically and chemically with later accreted populations in the inner = 0 halo. However, it remains unclear whether all subpopulations are born on halo-like radial orbits and retain memory of their formation. Or whether they acquire their radial orbits later during the galaxies' history, for example due to subsequent interactions or secular evolution.
To answer these questions, we repeat the analysis performed in Sections 4.1 and 4.2 along the cosmological merger tree of each galaxy. Briefly, we extract face-on azimuthal and radial velocities at each saved snapshot and apply the same cuts as in Section 4.1 to isolate the GSE-like population (| | ≤ 50 km s −1 , | | ≤ 240 km s −1 , and 2 ≤ ≤ 50 comoving kpc). We then repeat the assignment performed in Section 4.2 to identify subpopulations and show their respective numbers over time in Figure 8. We verified that all trends discussed in this section are qualitatively unchanged when plotting the number of stars which are both in a population at a given time and at = 0. We further observe significant variation from one time-step to the next due to shifts in galaxy centres and orientation affecting our kinematic and spatial cuts. Figure 8 thus shows the full time evolution data and a 3-pixel Gaussian-smoothed curve for readability (thin and thick lines, respectively).
We show on Figure 8 the timings of the 'A' and 'B' interactions ( infall to coalescence marked at the top). As expected, these time periods coincide with the rise of the specific subpopulations associated to these events. For all merger scenarios, the respective contributions of ancient populations is roughly in place at ≈ 1, after the most significant interactions have completed. The contaminating, younger-forming population grows later into the GSE-like feature (grey), consistent with contamination from disc stars and later-accreted populations (Section 4.2 and Appendix A).
Interaction populations (blue and green) exhibit the strongest time evolution, although it appears as a continuous process rather than being linked to localized triggers. In particular, the passage of 'B' does not correlate with a significant increase in older populations (violet, orange, blue, and green) on GSE-like orbits, as would be expected if the merger dynamically kicks a significant population of existing central stars onto radial orbits. Such increase is tentatively visible for the larger merger mass-ratios (centre and two right-most panels), but the time-step-to-time-step variability makes a robust interpretation challenging and changes around this time are comparable to later dynamical evolution. We further verified that the galaxy always has time to dynamically relax between 'A' and 'B', despite what appears as two mergers in quick succession (e.g. right-most panel) -in all scenarios, when 'A' has coalesced, the dynamical time √︁ 2 3 / (< ) at 2 1/2, is shorter than 20 Myr, whereas the first infall of 'B' happens at the fastest 150 Myr later.
We thus conclude that merger 'B' plays a key role in (i) delivering accreted stars onto radial orbits (red), (ii) forming stars on radial orbits by driving a starburst during its interaction (green), (iii) but minimally modifies the orbits of existing GSE-like stars.
Rather, dynamical evolution appears more gradually, as a slow increase in the number of stars on GSE-like orbits over time, which is particularly apparent for in-situ populations within our disc galaxies (blue and green lines in first four panels). These components can see their numbers double between = 1 and = 0 and overtake the more stable, accreted components. Our results thus stress the importance of capturing the full dynamical history of a galaxy when interpreting its = 0 observables, as late evolution can significantly affect the composition and mixture of GSE-like features. This slow evolution of the structure put in place at high-redshift points towards secular dynamical mechanisms within the galaxy, or repeated small perturbations over a Hubble time. We expect a general heating of old-stars orbits over time, either driven by physical perturbations such as satellite flybys (Sellwood & Binney 2002;Quillen et al. 2009) and scattering off gas and stellar clumps (e.g. van Donkelaar et al. 2022;Wu et al. 2022), or driven by numerical limitations due to the finite particle numbers and force resolution (e.g. Ludlow et al. 2019). However, it remains unclear whether such mechanisms would produce preferentially radially biased, halo-like orbits, motivating future cosmological studies quantifying their respective importance.

Disentangling high-redshift merger scenarios from a mock solar neighbourhood
Our results in Section 4 focus on galaxy-wide chemodynamical patterns to provide a complete and consistent comparison between merger scenarios. However, the data sets motivating our analysis, and with which we ultimately wish to compare, are only available within our Galaxy and around the particular location of the Solar neighbourhood. Furthermore, these data sets have non-trivial spatial and kinematic selection functions, which could potentially bias inferences of the merger history compared to galaxy-wide trends.
Transforming the stellar populations of cosmological simulations into resolved stars to account for such selection functions requires careful, dedicated modelling (e.g. Grand et al. 2018;Sanderson et al. 2020;Thomas et al. 2021;Lim et al. 2022), and is outside the scope of this work. In this section, we none the less provide a preliminary assessment of whether the multiple radially-biased halo populations highlighted in Section 4 would remain as overlapped chemodynamically in a mock Galactic survey. Our simulated galaxies are less massive than the Milky Way at = 0 ( ★ ≈ 2 × 10 10 M ) and lack some of its defining late-time features such as a bar ( Figure 3) and a Local Group environment. This limits direct, absolute comparisons with Galactic data, but still allows us to establish comparative trends between scenarios. We thus focus on our two intermediate scenarios 'Smaller = 2 merger' and 'Larger = 2 merger' that are closest to our Galaxy: they both show a regular disc structure necessary to define an equivalent Solar neighbourhood, have optical sizes compatible with the Milky Way's thin disc (≈ 2.5 kpc, e.g. McMillan 2017).
To mimic Galactic surveys of halo stars, we then centre on an arbitrary point in the disc plane at 8 kpc from the galactic centre and select stars with inclinations | | > 15 • and distances < 20 kpc. We then repeat the same kinematic cuts as in Section 4.1 to isolate the GSE-like population and perform the same assignment procedure as in Section 4.2 to link stars with early merger events. Figure 9 shows the resulting distributions of the total GSE-like population (purple) and each of its subcomponent (coloured lines matching Figure 6 A promising avenue for future work to distinguish these early merger scenarios could be to focus on the low-[Fe/H] tail of halo stars. The absolute number of stars with [Fe/H] ≤ −1.0 ([Fe/H] ≤ −0.5) in the mock solar neighbourhood does not vary significantly between the two scenarios, with 9677 stars and 10343 for the left-and right-hand columns (25,024 and 27,783 stars, respectively). But the different nature and timing of the progenitors could be potentially distinguished using specific peaks in the metallicity distribution functions, combined with additional chemical abundances sensitive to different enrichment time-scales (e.g. Al and Na; Belokurov & Kravtsov 2022;Feuillet et al. 2022, although see e.g. Buck et al. 2021 for uncertainties related to chemical enrichment). We plan in future work to post-process our simulations to produce such additional abundances and study in more details the patterns of metal-poor stars across our merger scenarios (see also Sestito et al. 2021). For now, we caution that linking = 0 metallicities in a GSE-like feature to the mass ratio of a single GSE-progenitor With this selection function, as with galaxy-wide samples, the magnitude of the GSE-like feature is conserved for both a small and large merger mass ratio at ≈ 2 (1:10 and 1:4 in left-and right-hand columns, respectively) and overall chemodynamical distributions strongly overlap (purple). In both cases, the overall population is a mixture of several subpopulations (colours match Figure 7), which do not separate clearly in chemical or phase-space to unequivocally identify each scenario. event is a delicate task. Robust inferences require broader exploration of cosmological merger histories to understand their potential degeneracies, and a careful account of their observables' in the Solar Neighbourhood to disentangle and weight the respective importance of mixed subpopulations.

How common are two high-redshift merger events for Milky-Way-mass dark matter haloes?
Our results highlight that multiple merger events at early times play a key role in shaping the final population of radially biased, halo stars at = 0. In particular, Section 4.2 highlights the important contribution of stars formed during an earlier interaction at ≈ 3 (merger 'A'), which is difficult to separate from ex-situ stars brought with the radial event at ≈ 2 explicitly targeted by our genetic modifications (merger 'B'). In Section 2.2 and Figure 1, we show that the mass growths of major progenitors across all merger scenarios are within the 2 population from ≈ 6. However, mass assembly can be achieved through both mergers and smooth accretion, and we now wish to assess the likelihood of having two significant mergers at such early times in a ΛCDM universe.
To answer this question formally, one should estimate the joint likelihood of two merger events of different mass ratios one after the other, across the population of Milky-Way-mass dark matter halo in a ΛCDM universe. However, such likelihood is not readily available due to the large parameter space that needs to be sampled across all possible ΛCDM merger histories. Instead, we extract the summary statistics introduced by Fakhouri et al. (2010) from the Millenium simulation suite Boylan-Kolchin et al. 2009) to describe the average population merger rate of dark matter haloes: where is the dark matter merger mass ratio and = 0.0144, 0 = 9.72 × 10 −3 , = 0.133, = −1.995, = 0.263, = 0.0993 are taken to the best-fitting parameters determined by Fakhouri et al. (2010).
Since we wish to verify the likelihood of early, major events, we integrate the merger rate between = 5 and = 1 along the median mass accretion history determined in Section 2.2, first for mass ratios greater than 1:5. We obtain (0.2 ≤ ≤ 1) = 2.01, which is consistent with our merger scenarios in which Milky-Waymass progenitors experience two significant events between = 5 and = 1 (Table 1). This confirms that GSE-like feature are likely to be common across the Milky-Way analogue population (see also Bignone et al. 2019;Mackereth et al. 2019;Fattahi et al. 2019;Elias et al. 2020;Dillamore et al. 2022), and that our findings that GSE-like stars have several origins due to multiple high-redshift progenitors is likely generic (see also Grand et al. 2020;Orkney et al. 2022).
The average merger rate in Equation 1 strongly depends on the merger mass ratio , rapidly suppressing increasingly major mergers (e.g. Fakhouri & Ma 2008;Fakhouri et al. 2010). Re-integrating the merger rate with ≥ 1:4, 1:3 and 1:2 gives = 1.67, = 1.28, and = 0.76, respectively. Our genetically modified scenario with the most massive merger undergoes two events with dark matter mass ratios greater than 1:3 between = 5 and = 1, and is thus a less likely realization according to this metric. Interestingly, such rarer formation scenario could provide a natural mechanism to populate the small population of low-mass, bulge-dominated, central galaxies surrounded by low-surface brightness, star-forming rings (e.g. Cappellari et al. 2011;Duc et al. 2015).
Despite this, we do not consider this last scenario (or any others) to be unlikely in a ΛCDM universe. Our genetic modifications generate minimal shifts in the likelihood of modified initial conditions compared to their reference, ensuring that all of them are highly consistent with the ΛCDM power spectrum (Section 2.2 and Figure 1). Furthermore, calculations of average merger rates suffer from systematic uncertainties due to issues in identifying haloes, trimming their merger trees and defining an infall mass ratio (e.g. Fakhouri & Ma 2008). For example, using an updated merger tree algorithm and different snapshot cadence, Poole et al. (2017) report a much flatter dependence on the merger mass ratio for the same fitting function, particularly affecting the major merger rate. Recomputing the integral using their best-fitting parameters (their table 2), we find = 1.90 for ≥ 1:3, i.e. a 50 per cent increase making it compatible at face value with our most extreme scenario.
Future studies more robustly assessing the likelihood of a given merger history would thus be highly beneficial. This could be achieved by leveraging the extremely large halo catalogues of modern simulations (e.g. Ishiyama et al. 2021;Maksimova et al. 2021) to sample additional parameter space explicitly and construct joint likelihoods, or using more advanced statistical models than fitting functions to summarize the complex information content of halo merger trees (e.g. Robles et al. 2022).

CONCLUSION
We present a suite of genetically modified, cosmological zoomed simulations of Milky-Way-mass galaxies, systematically studying how varying the mass ratio of a ≈ 2 merger impacts the = 0 chemodynamics of halo stars.
All simulations are evolved to = 0 with the extensive hydrodynamical, galaxy formation model of the project , that can capture the multiphase structure of the galaxy's ISM (≈ 20 pc) and a Hubble time of disc formation and dynamical evolution. Each genetically modified galaxy has, by construction, the same dynamical and stellar mass at = 0 ( Figure 1) and a similar large-scale environment (Figure 2). Comparing them with one another then provides a controlled, study isolating the signatures of a galaxy's early merger history on its chemodynamical structure at = 0. At fixed dynamical and stellar mass, modifications of the mass ratio at ≈ 2 result in a large diversity of galactic structure and morphology at = 0 ( Figure 3). Smaller mergers favour the growth of an extended, rotationally supported stellar disc, while, by contrast, larger mergers leads to more compact, bulge-dominated morphologies after a Hubble time. In fact, our largest merger mass ratio at = 2 lacks a well-defined galactic disc, and exhibits a counterrotating inner core surrounded by a faint, blue star-forming disc similar to those observed around nearby ellipticals (e.g. Duc et al. 2015).
Despite the large diversity in galaxy morphology and structure, each galaxy with well-ordered stellar rotation exhibits a population of radially biased, inner halo stars overlapping in energy with disc stars (Figure 4 and Figure 5). We detect this kinematic GSE-like population with similar magnitude in all merger scenarios, whether we enhance the ≈ 2 merger to become a major event or diminish it into a minor event, at odds with common interpretations linking GSE-like features to a single radially biased collision.
Furthermore, we show, for all merger mass ratios, stars within these GSE-like features have similar median peaks in ages (Figure 6), [Fe/H] (Figure 7), and overlapping 2D chemodynamical distributions (Appendix A). This demonstrates that the existence and magnitude of a GSE-like population does not directly and causally relate to the mass ratio of a single, early merger at ≈ 2, cautioning studies that aim to reconstruct the properties of an assumed single progenitor from = 0 chemodynamical data.
To identify the source of such degenerate signatures between multiple merger histories, we track back every GSE-like stars through the cosmological merger trees and pinpoint their origin. We find that halo stars on radial orbits at = 0 have diverse origins, originating from a mixture of in-situ and ex-situ subpopulations with overlapping chemodynamical distributions (Figure 7; see also Grand et al. 2020;Orkney et al. 2022).
In particular, we find that a significant fraction of the = 0 feature originates from stars formed during an earlier, high-redshift interaction at ≈ 3, which have similar chemical abundances to the accreted population at ≈ 2. As we decrease the significance of the progenitor at ≈ 2 using genetic modifications, the contribution of its accreted component diminishes accordingly. However, this is compensated by the growth of the earlier, starburst population to conserve the same total stellar and dynamical mass, resulting in similar chemodynamics for the GSE-like feature.
Our study highlights the importance of modelling the full cosmological formation scenario when interpreting Galactic data at = 0. We explicitly demonstrate the importance of capturing the multiple high-redshift interactions, as well as the response of the central galaxy to each of them to obtain a full account of the several contributors to the final metal-poor, radial halo population. Robust inferences of past merger progenitors thus require a wide exploration of cosmological early scenarios, sampling a spectrum of progenitors, mass ratios, and infall geometry, to establish robust signatures distinguishing them. As shown by this study, the genetic modification approach combined with high-resolution zoomed simulations offers great potential to efficiently achieve this objective, and inform us on the cosmological formation scenario of our Milky Way. In particular, further genetic modifications controlling the angular momentum accretion history from the initial conditions (Cadiou et al. , 2022 could allow us to vary the infall times and orbital parameters of GSE progenitors, in addition to their mass ratios studied here. Despite the mechanisms and chemodynamical degeneracies between merger scenarios highlighted in this study, it remains unclear how to link one of our cosmological scenario (if any) to our Milky Way. All of our genetically modified galaxies have lower final stellar masses ( ★ ≈ 2 × 10 10 M ) than our Galaxy (≥ 6 × 10 10 M ), and blindly upscaling stellar masses across the merger tree would bring our GSE progenitors towards the massive end of commonly inferred values ( ★ ≥ 10 9 M ), without changing their merger mass ratios. Our quieter genetically modified histories that prefer smaller dwarfs could then relate more closely to the potential history of our Milky Way, although a key finding from our study is that inferring past merger mass ratios from = 0 chemodynamical data is challenging. Furthermore, the numerical modelling of the stellar mass growth of galaxies, particularly that of dwarf galaxies at ≈ 2, remains highly uncertain (e.g. Somerville & Davé 2015;Naab & Ostriker 2017). This adds further uncertainty and scatter to the dark matter halo masses that should host a given stellar-mass progenitor, and thus their dynamical impact on the protogalaxy. In addition to explorations of cosmological merger scenarios, future studies scanning through galaxy formation models will be key to understand how such uncertainties affect the reconstruction of our Milky Way's history.
Furthermore, our galaxies lack defining late-time features of our Galaxy (e.g. a bar), pointing to different evolution, gas accretion, and star formation activity at late times ( < 1) between our scenarios and the Milky Way. Since early and late histories are highly correlated in a ΛCDM Universe, constraining the later evolution of our Galaxy will in turn reduce the available freedom in possible early, high-redshift history. The structure of the Milky Way's discs (e.g. Belokurov & Kravtsov 2022;Xiang & Rix 2022; Figure 3) and their potential accreted components (e.g. Ruchti et al. 2015;Feuillet et al. 2022), the density profile of the stellar halo (e.g. Deason et al. 2011Deason et al. , 2014Han et al. 2022), or the distribution of globular clusters (e.g. Kruĳssen et al. 2020, although see Pagnini et al. 2022) offer us additional and complementary clues that can be combined with the GSE's properties to constrain the cosmological formation scenario of the Milky Way.
The richness of Galactic data also offers numerous additional promising avenues to find more detailed diagnosis that would allow us to distinguish each merger scenario. Notably, using additional abundance ratios tracking different enrichment time-scales (e.g. Belokurov & Kravtsov 2022;Feuillet et al. 2022;Matsuno et al. 2022), combined with more detailed clustering techniques in the high-dimensional chemodynamical space to isolate subpopulations (e.g. Myeong et al. 2022;Dodd et al. 2023) will become ever-more powerful with the next-generation of spectroscopic surveys and the improvement in data quality and sample sizes in the coming years (e.g. WEAVE, 4MOST).

APPENDIX A: ADDITIONAL CHEMODYNAMICAL CHARACTERIZATION OF THE GSE-LIKE POPULATIONS
In Section 4, we establish that, despite a large diversity in final morphology and past merger mass ratio at ≈ 2, each galaxy has similar GSE-like kinematic features at fixed stellar mass (Section 4.1). We further show in Section 4.2 that GSE-like stars, in all merger scenarios at ≈ 2, have broad age and [Fe/H] distributions that peak around the same median, due to overlapping populations with varying cosmological origins. In this Appendix, we show that distinguishing subpopulations within each merger scenario, or merger scenarios from one another, remains difficult even when leveraging additional chemodynamical data for the GSE-like stars. Figure A1 shows the distributions of E tot and [ /Fe] of each GSE subpopulation identified in Section 4.2. As in Section 4, we recover strong similarities in the overall GSE-like population (purple) between the different merger scenarios, with similar peaks in abundances -[ /Fe] = 0.48, 0.48, 0.48, 0.47, 0.47, from left to right, respectively -and broad coinciding distributions in E tot .
Breaking down the overall population into its subpopulations (colours matching Figure 7), we first confirm that, for all merger scenarios, the contaminating population (grey) has high E tot and low [ /Fe], combining with its higher [Fe/H]. This is consistent with our interpretation that this contamination primarily arises from outer disc and halo stars that formed after the ancient structures of interest.
Focusing on the subpopulations originating at high-redshift (colours), we recover a strong overlap in E tot and [ /Fe], adding to the kinematic and chemical overlap already highlighted in Section 4. In particular, stars formed during early interactions (light blue and green) can populate high-E tot orbits and have similar [ /Fe] ratios than the accreted population of the genetically modified merger at ≈ 2 ('B'; red).
To further visualize this, we show in Figure A2 the positions of subpopulations in the E totz , [ /Fe]-[Fe/H], and − planes (top, middle and bottom, respectively) compared to the overall population of stars within the galaxy (grey background). We estimate respective 50% isodensity contour using a 2D kernel density estimate with Gaussian smoothing according to Scott's rule. The small sample sizes, large spread and non-Gaussian nature of chemodynam-ical distributions can lead to noisy KDE estimates 4 , and generates the large isodensity contours for very early populations that have large spreads in [ /Fe] and small sample sizes. For completeness, we choose to show all subpopulations, even when noisy, but stress that they might not contribute a large fraction of the total GSE-like feature (absolute contributions can be estimated from Figure 7 and A1).
As with 1D distributions, we find substantial similarities between merger scenarios (left to right). Individual subpopulations have distinct peaks and a clear progressing sequence in [ /Fe]-[Fe/H]. However, the main contributors to the metal-poor GSE-like population (blue formed during the earlier interaction 'A', and red brought with the merger body 'B') all significantly overlap in these planes. We thus do not find differences in chemodynamics of the GSE-like population that unequivocally allows us to differentiate merger mass ratio at ≈ 2.

APPENDIX B: COMPARING GALAXY-WIDE AND GSE STAR FORMATION HISTORIES
In Section 4.2, we show that star formation peaks in the GSE-like population coincide in timing with the coalescence of early mergers (Figure 6), when a significant fraction of its stars are formed (Figure 7). In this Appendix, we quantify whether this surplus in star formation reflect a general galaxy-wide enhancement, or is specific to stars on halo-like, radial orbits at = 0.
We show in Figure A3 the star formation histories of all inner stars ( ≤ 50 kpc; black) for each merger scenario (mass ratio at ≈ 2 increasing from left to right), compared to those of GSE-like stars (purple; see also Figure 6). We recover that early interactions ( infall to coalescence marked at the top) drives large, galaxy-wide star formation enhancements consistent with merger induced starbursts. For lower mass ratios at ≈ 2 (first three panels), the elevated star formation coincide with peaks in star formation in the GSE-like population, reflecting the dominant contribution of starburst stars in this population at = 0 ( Figure 7).
For larger mass ratios (right-hand panels), the mapping between galaxy-wide and GSE star formation is less clear, particularly for interaction 'B' (red). This is particularly visible for the largest mass ratio (right-most panel), where the galaxy undergoes a large starburst during interaction 'B' (black), but stars formed at this time do not populate radially biased orbits at = 0 (purple). In fact, the origin of GSE-like at = 0 population is dominated by accreted stars previously formed in the merging body, rather than stars formed during the interaction (Figure 7). This paper has been typeset from a T E X/L A T E X file prepared by the author.    Figure 4 (purple) and their high-redshift origin (colours matching Figure 7). In all merger scenarios (increasing mass ratio from left to right), in-situ and ex-situ contributors substantially overlap chemically and kinematically, and cannot be easily distinguished (see also Figure A2). The contaminating population (grey) lies preferentially at high E tot and low [ /Fe].  Figure 7. Contours show 50% isodensities in these planes, normalized to each subpopulation's total (their absolute contribution to the overall GSE-like feature can be estimated from Figure 7 and A1). The most important subpopulations to the metal-poor GSE feature (blue and red) always overlap in all planes, making it difficult to distinguish increasing merger mass ratios at ≈ 2 in the GSE-like progenitor (left to right). The grey histogram shows the distribution of all inner stars ( < 50 kpc). Largest z = 2 merger Figure A3. Same as Figure 6, but comparing the galactic and GSE-like star formation histories. Peaks of star formation in GSE-like stars broadly map onto galaxy-wide star formation enhancements during early interactions (marked at the top). However, the opposite is not verified -galaxy-wide starbursts do not necessarily form stars that end up on radial halo orbits at = 0. This is particularly visible for the largest mass ratio (right-most panel) where interaction 'B' drives a large galaxy-wide starburst that is missing in the GSE-like population.