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

Martin P. Rey, Oscar Agertz, Tjitske K. Starkenburg, Florent Renaud, Gandhali D. Joshi, Andrew Pontzen, Nicolas F. Martin, Diane K. Feuillet, Justin I. Read

<figure><img src="image_1.png" /></figure>

## To cite this version:

Martin P. Rey, Oscar Agertz, Tjitske K. Starkenburg, Florent Renaud, Gandhali D. Joshi, et al. VINTERGATAN-GM: The cosmological imprints of early mergers on Milky-Way-mass galaxies. Monthly Notices of the Royal Astronomical Society, 2023, 521 (1), pp.995-1012. 10.1093/mn-ras/stad513. hal-03901573

HAL Id: hal-03901573

https://hal.science/hal-03901573

Submitted on 24 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_1.png" /></figure>

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

Martin P. Rey, Oscar Agertz, Tjitske K. Starkenburg, Florent Renaud, Gandhali D. Joshi, Andrew Pontzen, Nicolas F. Martin, Diane K. Feuillet, Justin I. Read

<figure><img src="image_2.png" /></figure>

## To cite this version:

Martin P. Rey, Oscar Agertz, Tjitske K. Starkenburg, Florent Renaud, Gandhali D. Joshi, et al. VINTERGATAN-GM: The cosmological imprints of early mergers on Milky-Way-mass galaxies. Monthly Notices of the Royal Astronomical Society, 2023, 521 (1), pp.995-1012. 10.1093/mn-ras/stad513. hal-03901573

HAL Id: hal-03901573

https://hal.science/hal-03901573v1

Submitted on 24 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>
# VINTERGATAN-GM: The cosmological imprints of early mergers on Milky-Way-mass galaxies

Martin P. Rey<sup>1,1*</sup> Oscar Agertz<sup>1,2</sup> Tjitske K. Starkenburg,<sup>3</sup> Florent Renaud<sup>1,2</sup> Gandhali D. Joshi<sup>1,4</sup> Andrew Pontzen,<sup>4</sup> Nicolas F. Martin<sup>1,5,6</sup> Diane K. Feuillet<sup>2</sup> and Justin I. Read<sup>1,7</sup>

Accepted 2023 February 12. Received 2023 February 10; in original form 2022 December 2

## ABSTRACT

We present a new suite of cosmological zoom-in hydrodynamical (≈ 20 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_{\odot}$), incrementally increasing the stellar mass ratio of a $z \approx 2$ 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, [Fe/H], and [α/Fe] distribution whether the $z \approx 2$ merger is more minor or major. This arises because a smaller ex-situ population at $z \approx 2$ is compensated by a larger population formed in an earlier merger-driven starburst whose contribution to the GES can grow dynamically over time, and with both populations strongly overlapping in the [Fe/H] — [α/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.

**Key words:** methods: numerical – Galaxy: formation – Galaxy: halo – galaxies: formation – galaxies: kinematics and dynamics.

## 1 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 2016a, b, 2018, 2021, 2022). Combined with chemical abundances and radial velocities acquired by large, spectroscopic surveys (e.g. APOGEE; GALAH; H3; 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.g. Belokurov et al. 2018; Helmi et al. 2018; Myeong et al. 2019; Kruijssen 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; Haywood et al. 2018; Koppelman, Helmi & Veljanoski 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 $z \approx 2$.

However, the exact mass-scale of the GSE progenitor remains to be pinpointed, with stellar mass estimates extending over an order of magnitude ($M_* \sim 10^8 - 10^9 M_{\odot}$; e.g. Bonaca et al. 2020; Feuillet et al. 2020; Kruijssen et al. 2020; Mackereth & Bovy 2020; Naidu et al. 2020, 2021; Limberg 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 $z = 0$ is a pure population that can be robustly

<footnote n="1"><sup>1</sup>Sub-department of Astrophysics, University of Oxford, DWB, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="2"><sup>2</sup>Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-221 00 Lund, Sweden</footnote>

<footnote n="3"><sup>3</sup>Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA</footnote>

<footnote n="4"><sup>4</sup>Department of Physics and Astronomy, University College London, London WC1E 6BT, UK</footnote>

<footnote n="5"><sup>5</sup>Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="6"><sup>6</sup>Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="7"><sup>7</sup>Department of Physics, University of Surrey, Guildford GU2 7XH, UK</footnote>

<footnote>*E-mail: martin.rey@physics.ox.ac.uk</footnote>
# VINTERGATAN-GM: The cosmological imprints of early mergers on Milky-Way-mass galaxies

Martin P. Rey<sup>1,★</sup> Oscar Agertz<sup>1,2</sup> Tjitske K. Starkenburg,<sup>3</sup> Florent Renaud<sup>1,2</sup> Gandhali D. Joshi<sup>1,4</sup> Andrew Pontzen,<sup>4</sup> Nicolas F. Martin<sup>1,5,6</sup> Diane K. Feuillet<sup>2</sup> and Justin I. Read<sup>1,7</sup>

Accepted 2023 February 12. Received 2023 February 10; in original form 2022 December 2

## ABSTRACT

We present a new suite of cosmological zoom-in hydrodynamical (≈ 20 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_{\odot}$), incrementally increasing the stellar mass ratio of a $z \approx 2$ 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, [Fe/H], and [α/Fe] distribution whether the $z \approx 2$ merger is more minor or major. This arises because a smaller ex-situ population at $z \approx 2$ is compensated by a larger population formed in an earlier merger-driven starburst whose contribution to the GES can grow dynamically over time, and with both populations strongly overlapping in the [Fe/H] — [α/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.

**Key words:** methods: numerical – Galaxy: formation – Galaxy: halo – galaxies: formation – galaxies: kinematics and dynamics.

## 1 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 2016a, b, 2018, 2021, 2022). Combined with chemical abundances and radial velocities acquired by large, spectroscopic surveys (e.g. APOGEE; GALAH; H3; 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.g. Belokurov et al. 2018; Helmi et al. 2018; Myeong et al. 2019; Kruijssen 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; Haywood et al. 2018; Koppelman, Helmi & Veljanoski 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 $z \approx 2$.

However, the exact mass-scale of the GSE progenitor remains to be pinpointed, with stellar mass estimates extending over an order of magnitude ($M_* \sim 10^8 – 10^9 M_{\odot}$; e.g. Bonaca et al. 2020; Feuillet et al. 2020; Kruijssen et al. 2020; Mackereth & Bovy 2020; Naidu et al. 2020, 2021; Limberg 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 $z = 0$ is a pure population that can be robustly

<footnote n="1"><sup>1</sup>Sub-department of Astrophysics, University of Oxford, DWB, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="2"><sup>2</sup>Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-221 00 Lund, Sweden</footnote>

<footnote n="3"><sup>3</sup>Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA</footnote>

<footnote n="4"><sup>4</sup>Department of Physics and Astronomy, University College London, London WC1E 6BT, UK</footnote>

<footnote n="5"><sup>5</sup>Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="6"><sup>6</sup>Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="7"><sup>7</sup>Department of Physics, University of Surrey, Guildford GU2 7XH, UK</footnote>

<footnote>*E-mail: martin.rey@physics.ox.ac.uk</footnote>
linked to one single event, or whether it contains a superposition of multiple population with distinct origins (e.g. Grand et al. 2020; Donlon II et al. 2022; Donlon & Newberg 2022; Khoperskov et al. 2022a; Myeong et al. 2022; Orkney et al. 2022). 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. 2020, 2021) but could also be of entirely distinct origin (e.g. Amarante et al. 2020; Donlon II 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 ($N \gtrsim 10^6$; e.g. Sellwood 2013; Ludlow et al. 2019, 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. 2017; Buck et al. 2020; Font et al. 2020; Agertz et al. 2021; Applebaum et al. 2021; Bird et al. 2021; Grand et al. 2021; Khoperskov et al. 2022a; Wetzel et al. 2022). With such tools, we can now quantify the frequency of chemodynamical patterns at $z=0$, link them to specific events in each galaxy’s history, and inform the reconstruction of our Galaxy’s formation (e.g. Bignone, Helmi & Tissera 2019; Fattahi et al. 2019; Mackereth 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 $z=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 ($\approx 20$ pc), cosmological zoomed simulations using a physical model that can successfully reproduce the chemodynamical structure of Milky-Way-like galaxies (Agertz et al. 2021; Renaud et al. 2021a, b); and (ii) the genetic modification approach (Roth, Pontzen & Peiris 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, Crain & Pontzen 2021), or the angular momentum accretion (e.g. Cadiou, Pontzen & Peiris 2022). Each modified version differs minimally from the original scenario, for example conserving the total dynamical mass at $z=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 $\approx 10^{12}$ M$_\odot$ dark matter halo which experiences a merger at $z \approx 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 $z=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 VINTERGATAN project (Agertz et al. 2021). 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 preparation).

We present the simulation suite in Section 2, and show in Section 3 that modifying the mass ratio of the GSE-like event at $z \approx 2$ leads to markedly distinct galactic morphologies at $z=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 $z=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.

## 2 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 $z \approx 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 $z=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.

### 2.1 Numerical setup and galaxy formation physics

We construct cosmological, zoomed initial conditions using the GENETIC software (Stopyra et al. 2021) and cosmological parameters $\Omega_m = 0.3139$, $h = 0.6727$, $\sigma_8 = 0.8440$, $n_s = 0.9645$ (Planck Collaboration XIII 2016). Starting from a dark matter-only cosmological volume with a box size $50 h^{-1} \text{Mpc} \approx 73 \text{Mpc}$ and mass resolution $m_{\text{DM}} = 1.2 \times 10^8 \text{ M}_\odot$, we select a target reference halo with Milky-Way virial mass ($M_{200} \approx 10^{12} \text{ M}_\odot$) and no neighbours more massive within $5 r_{200}$, where $r_{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 $m_{\text{DM}} = 2.0 \times 10^5 \text{ M}_\odot$, 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 $z=99$ (Zel’dovich 1970), before we start following the evolution of dark matter, stars, and gas to $z=0$ with the adaptive mesh refinement code RAMSES (Teyssier 2002).
linked to one single event, or whether it contains a superposition of multiple population with distinct origins (e.g. Grand et al. 2020; Donlon II et al. 2022; Donlon & Newberg 2022; Khoperskov et al. 2022a; Myeong et al. 2022; Orkney et al. 2022). 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. 2020, 2021) but could also be of entirely distinct origin (e.g. Amarante et al. 2020; Donlon II 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 ($N \gtrsim 10^6$; e.g. Sellwood 2013; Ludlow et al. 2019, 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. 2017; Buck et al. 2020; Font et al. 2020; Agertz et al. 2021; Applebaum et al. 2021; Bird et al. 2021; Grand et al. 2021; Khoperskov et al. 2022a; Wetzel et al. 2022). With such tools, we can now quantify the frequency of chemodynamical patterns at $z=0$, link them to specific events in each galaxy's history, and inform the reconstruction of our Galaxy's formation (e.g. Bignone, Helmi & Tissera 2019; Fattahi et al. 2019; Mackereth 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 $z=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 ($\approx 20$ pc), cosmological zoomed simulations using a physical model that can successfully reproduce the chemodynamical structure of Milky-Way-like galaxies (Agertz et al. 2021; Renaud et al. 2021a, b); and (ii) the genetic modification approach (Roth, Pontzen & Peiris 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, Crain & Pontzen 2021), or the angular momentum accretion (e.g. Cadiou, Pontzen & Peiris 2022). Each modified version differs minimally from the original scenario, for example conserving the total dynamical mass at $z=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 $\approx 10^{12}$ M$_\odot$ dark matter halo which experiences a merger at $z \approx 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 $z=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 VINTERGATAN project (Agertz et al. 2021). 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 preparation).

We present the simulation suite in Section 2, and show in Section 3 that modifying the mass ratio of the GSE-like event at $z \approx 2$ leads to markedly distinct galactic morphologies at $z=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 $z=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.

## 2 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 $z \approx 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 $z=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.

### 2.1 Numerical setup and galaxy formation physics

We construct cosmological, zoomed initial conditions using the GENETIC software (Stopyra et al. 2021) and cosmological parameters $\Omega_m = 0.3139$, $h = 0.6727$, $\sigma_8 = 0.8440$, $n_s = 0.9645$ (Planck Collaboration XIII 2016). Starting from a dark matter-only cosmological volume with a box size $50 h^{-1} \text{Mpc} \approx 73 \text{Mpc}$ and mass resolution $m_{\text{DM}} = 1.2 \times 10^8 \text{M}_\odot$, we select a target reference halo with Milky-Way virial mass ($M_{200} \approx 10^{12} \text{M}_\odot$) and no neighbours more massive within $5 r_{200}$, where $r_{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 $m_{\text{DM}} = 2.0 \times 10^5 \text{M}_\odot$, 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 $z=99$ (Zel'dovich 1970), before we start following the evolution of dark matter, stars, and gas to $z=0$ with the adaptive mesh refinement code RAMSES (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, Spruce & Speares 1994) assuming an ideal gas equation of state with adiabatic index $\gamma = 5/3$. Our Lagrangian refinement strategy allows us to reach a spatial resolution of 20 pc throughout the galaxy’s ISM (Agertz et al. 2021). 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 RAMSES version. Gas with $n_H \ge 0.01 \text{ cm}^{-3}$ self-shields from this heating source (Aubert & Teyssier 2010; Rosdahl & Blaizot 2012), allowing it to condense to densities $\rho \ge 100 \text{ m}_\text{p} \text{ cm}^{-3}$ and temperatures $T \le 100 \text{ 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 \text{ M}_\odot$ 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. 2013; Agertz & Kravtsov 2015; Agertz et al. 2020). 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, Faucher-Giguère & Quataert 2015; Agertz et al. 2020, 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, Mannucci & Brandt 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, Knebe & Gibson 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 PYNBODY (Pontzen et al. 2013) and TANGOS (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, $M_\star$, by summing the stellar mass within $r_{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, $r_{1/2}$, $v$ is then derived along a random line of sight. We define the iron and $\alpha$ 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 $\alpha$ abundance as it dominates the mass fraction. We compute total metallicities from these two elements following Kim et al. (2014), equation (4).

## 2.2 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 $M_{200} \approx 1 \times 10^{12} \text{ M}_\odot$ and which experiences an early, large interaction with a radially biased orbit (infall at $z \approx 2$ with a $M_{200}$ ratio of 1:6 and a radial-to-tangential velocity ratio between the progenitors $v_r/v_\phi = 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).

Fig. 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 ($z \ge 5$) and late ($z \le 0.5$) accretion histories, converging by design to similar dynamical masses at $z = 0$ to within 10 per cent (Rey & Starkenburg 2022). Further, Fig. 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 $\Lambda$CDM cosmology (see also Pontzen et al. 2017; Rey, Pontzen & Saintonge 2019a for further visuals).

Mass growth histories around $z \approx 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 $\Lambda$CDM universe. The relative likelihood of each new modified initial condition to the reference, $\Delta\chi^2$, remain small compared to the available number of degrees of freedom – $\Delta\chi^2 = -3.6, -2.77, -0.45, +0.3$ for increasing merger mass ratio compared to $\approx 10^6$ modes in the zoom region – ensuring their compatibility with the $\Lambda$CDM power spectrum. We further compute the median and 68–95 per cent 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 Fig. 1, see Rey & Starkenburg 2022 for further details on the computation). All our mass growth histories are within the $1\sigma$ contour of the overall population around $z \approx 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 $\Lambda$CDM merger rates).

To quantify differences in merger histories, we extract merger events in the reference scenario that (i) bring at least $M_\star \ge 10^8 \text{ M}_\odot$ of accreted material at infall and with (ii) merger mass ratios more significant than 1:30 in $M_{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 Figs 1 and 2), which, as we will see in Section 4 all play a role in defining the $z = 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
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, Spruce & Speares 1994) assuming an ideal gas equation of state with adiabatic index $\gamma = 5/3$. Our Lagrangian refinement strategy allows us to reach a spatial resolution of 20 pc throughout the galaxy’s ISM (Agertz et al. 2021). 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 RAMSES version. Gas with $n_H \ge 0.01 \text{ cm}^{-3}$ self-shields from this heating source (Aubert & Teyssier 2010; Rosdahl & Blaizot 2012), allowing it to condense to densities $\rho \ge 100 \text{ m}_\text{p} \text{ cm}^{-3}$ and temperatures $T \le 100 \text{ 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 \text{ M}_\odot$ 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. 2013; Agertz & Kravtsov 2015; Agertz et al. 2020). 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, Faucher-Giguère & Quataert 2015; Agertz et al. 2020, 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, Mannucci & Brandt 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, Knebe & Gibson 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 PYNBODY (Pontzen et al. 2013) and TANGOS (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, $M_\star$, by summing the stellar mass within $r_{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, $r_{1/2}$, $v$ is then derived along a random line of sight. We define the iron and $\alpha$ 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 $\alpha$ abundance as it dominates the mass fraction. We compute total metallicities from these two elements following Kim et al. (2014), equation (4).

## 2.2 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 $M_{200} \approx 1 \times 10^{12} \text{ M}_\odot$ and which experiences an early, large interaction with a radially biased orbit (infall at $z \approx 2$ with a $M_{200}$ ratio of 1:6 and a radial-to-tangential velocity ratio between the progenitors $v_r/v_\phi = 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).

Fig. 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 ($z \ge 5$) and late ($z \le 0.5$) accretion histories, converging by design to similar dynamical masses at $z = 0$ to within 10 per cent (Rey & Starkenburg 2022). Further, Fig. 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 $\Lambda$CDM cosmology (see also Pontzen et al. 2017; Rey, Pontzen & Saintonge 2019a for further visuals).

Mass growth histories around $z \approx 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 $\Lambda$CDM universe. The relative likelihood of each new modified initial condition to the reference, $\Delta\chi^2$, remain small compared to the available number of degrees of freedom – $\Delta\chi^2 = -3.6, -2.77, -0.45, +0.3$ for increasing merger mass ratio compared to $\approx 10^6$ modes in the zoom region – ensuring their compatibility with the $\Lambda$CDM power spectrum. We further compute the median and 68–95 per cent 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 Fig. 1, see Rey & Starkenburg 2022 for further details on the computation). All our mass growth histories are within the $1\sigma$ contour of the overall population around $z \approx 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 $\Lambda$CDM merger rates).

To quantify differences in merger histories, we extract merger events in the reference scenario that (i) bring at least $M_\star \ge 10^8 \text{ M}_\odot$ of accreted material at infall and with (ii) merger mass ratios more significant than 1:30 in $M_{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 Figs 1 and 2), which, as we will see in Section 4 all play a role in defining the $z = 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
<table><caption>**Table 1.** Properties of the three main merger events (A, B, and C, marked in Figs 1 and 2) cross-matched across each genetically modified scenario (individual columns). For each event, we quote the redshift at which the infaller's centre is last outside the main body's $r_{200}$, $z_{\text{infall}}$. At this time, we report the progenitor's virial mass, stellar mass, mass ratios with the central body and radial to tangential velocity ratio at $r_{200}$. We identify the redshift of the interaction's first pericentre passage and final coalescence (see Section 2.2). Self-consistently cross-matching structures between simulations can lead to mass ratios less than unity if the identity of the major progenitor swaps at very early times (merger 'A' in first two columns).</caption><thead><tr><th></th><th>Smallest z = 2 merger</th><th>Smaller z = 2 merger</th><th>Reference</th><th>Larger z = 2 merger</th><th>Largest z = 2 merger</th></tr></thead><tbody><tr><td>Target progenitor B</td><td>$z_{\text{infall}} = 1.94$</td><td>$z_{\text{infall}} = 2.03$</td><td>$z_{\text{infall}} = 1.99$</td><td>$z_{\text{infall}} = 1.99$</td><td>$z_{\text{infall}} = 2.24$</td></tr><tr><td>Infall $M_{200}$</td><td>$2.2 \times 10^{10} M_{\odot}$</td><td>$2.3 \times 10^{10} M_{\odot}$</td><td>$3.8 \times 10^{10} M_{\odot}$</td><td>$7.4 \times 10^{10} M_{\odot}$</td><td>$1.1 \times 10^{11} M_{\odot}$</td></tr><tr><td>Infall $M_{\star}$</td><td>$2.2 \times 10^{8} M_{\odot}$</td><td>$3.1 \times 10^{8} M_{\odot}$</td><td>$6.0 \times 10^{8} M_{\odot}$</td><td>$1.2 \times 10^{9} M_{\odot}$</td><td>$2.3 \times 10^{9} M_{\odot}$</td></tr><tr><td>Ratio $M_{200}$</td><td>1:10</td><td>1:9.8</td><td>1:6.0</td><td>1:2.9</td><td>1:2.1</td></tr><tr><td>Ratio $M_{\star}$</td><td>1:24</td><td>1:15</td><td>1:8.1</td><td>1:4.3</td><td>1:2.1</td></tr><tr><td>Infall $v_r/v_{\theta}$</td><td>9.5</td><td>9.1</td><td>16.</td><td>8.8</td><td>5.5</td></tr><tr><td>$z_{\text{pericentre}}$</td><td>1.70</td><td>1.76</td><td>1.79</td><td>1.80</td><td>1.98</td></tr><tr><td>$z_{\text{coalescence}}$</td><td>1.40</td><td>1.45</td><td>1.52</td><td>1.47</td><td>1.65</td></tr><tr><td>Earlier progenitor A</td><td>$z_{\text{infall}} = 3.11$</td><td>$z_{\text{infall}} = 3.17$</td><td>$z_{\text{infall}} = 3.35$</td><td>$z_{\text{infall}} = 3.35$</td><td>$z_{\text{infall}} = 3.49$</td></tr><tr><td>Infall $M_{200}$</td><td>$9.6 \times 10^{10} M_{\odot}$</td><td>$7.7 \times 10^{10} M_{\odot}$</td><td>$5.8 \times 10^{10} M_{\odot}$</td><td>$5.1 \times 10^{10} M_{\odot}$</td><td>$3.7 \times 10^{10} M_{\odot}$</td></tr><tr><td>Infall $M_{\star}$</td><td>$6.2 \times 10^{8} M_{\odot}$</td><td>$5.2 \times 10^{8} M_{\odot}$</td><td>$4.5 \times 10^{8} M_{\odot}$</td><td>$3.8 \times 10^{8} M_{\odot}$</td><td>$3.1 \times 10^{8} M_{\odot}$</td></tr><tr><td>Ratio $M_{200}$</td><td>1:0.7</td><td>1:0.9</td><td>1:1.3</td><td>1:1.4</td><td>1:1.9</td></tr><tr><td>Ratio $M_{\star}$</td><td>1:1.8</td><td>1:2.1</td><td>1:2.4</td><td>1:1.8</td><td>1:2.4</td></tr><tr><td>Infall $v_r/v_{\theta}$</td><td>9.5</td><td>4.1</td><td>2.3</td><td>2.3</td><td>2.0</td></tr><tr><td>$z_{\text{pericentre}}$</td><td>2.74</td><td>2.72</td><td>2.86</td><td>2.94</td><td>3.00</td></tr><tr><td>$z_{\text{coalescence}}$</td><td>2.17</td><td>2.30</td><td>2.34</td><td>2.34</td><td>2.38</td></tr><tr><td>Later progenitor C</td><td>$z_{\text{infall}} = 1.94$</td><td>$z_{\text{infall}} = 1.86$</td><td>$z_{\text{infall}} = 1.73$</td><td>$z_{\text{infall}} = 1.70$</td><td>$z_{\text{infall}} = 1.62$</td></tr><tr><td>Infall $M_{200}$</td><td>$1.5 \times 10^{10} M_{\odot}$</td><td>$1.8 \times 10^{10} M_{\odot}$</td><td>$1.8 \times 10^{10} M_{\odot}$</td><td>$1.8 \times 10^{10} M_{\odot}$</td><td>$1.7 \times 10^{10} M_{\odot}$</td></tr><tr><td>Infall $M_{\star}$</td><td>$1.7 \times 10^{8} M_{\odot}$</td><td>$2.8 \times 10^{8} M_{\odot}$</td><td>$2.3 \times 10^{8} M_{\odot}$</td><td>$2.8 \times 10^{8} M_{\odot}$</td><td>$2.2 \times 10^{8} M_{\odot}$</td></tr><tr><td>Ratio $M_{200}$</td><td>1:15</td><td>1:14</td><td>1:18</td><td>1:19</td><td>1:22</td></tr><tr><td>Ratio $M_{\star}$</td><td>1:31</td><td>1:20</td><td>1:29</td><td>1:26</td><td>1:50</td></tr><tr><td>Infall $v_r/v_{\theta}$</td><td>0.5</td><td>0.6</td><td>0.5</td><td>0.4</td><td>0.6</td></tr><tr><td>$z_{\text{pericentre}}$</td><td>1.61</td><td>1.59</td><td>1.48</td><td>1.41</td><td>1.35</td></tr><tr><td>$z_{\text{coalescence}}$</td><td>0.09</td><td>0.04</td><td>0.17</td><td>0.22</td><td>0.20</td></tr></tbody></table>

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, $z_{\text{pericentre}}$, and end of the interaction, $z_{\text{coalescence}}$, using the high-cadence simulation movies.

Our genetic modifications explicitly target merger 'B' (second row in Table 1; red in Figs 1 and 2), which in the reference case is the last major event until $z = 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 ($v_r/v_{\theta} \ge 5$). The range of merging stellar masses scanned by our suite ($2 \times 10^8 M_{\odot} \le M_{\star} \le 2 \times 10^9 M_{\odot}$) accurately brackets the proposed masses for a potential GSE progenitor (e.g. Belokurov et al. 2018; Helmi et al. 2018; Feuillet et al. 2020; Kruijssen et al. 2020; Mackereth & Bovy 2020; Naidu et al. 2020, 2021; Limberg 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 $z = 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 $\Lambda$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 Fig. 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, Pontzen & Peiris 2021; Pontzen et al. 2021). However, the differences in timings ($\approx 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 $z = 0$, but we note that their central galaxies also match in their final stellar mass ($M_{\star} = 1.8 \times 10^{10} M_{\odot}$ within 20 per cent of each other; Fig. 1). All final stellar masses are compatible at $1\sigma$ with empirical model predictions for this halo mass (Moster, Naab & White 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 $z \approx 2$, while fixing the large-
<table><caption>**Table 1.** Properties of the three main merger events (A, B, and C, marked in Figs 1 and 2) cross-matched across each genetically modified scenario (individual columns). For each event, we quote the redshift at which the infaller's centre is last outside the main body's $r_{200}$, $z_{\text{infall}}$. At this time, we report the progenitor's virial mass, stellar mass, mass ratios with the central body and radial to tangential velocity ratio at $r_{200}$. We identify the redshift of the interaction's first pericentre passage and final coalescence (see Section 2.2). Self-consistently cross-matching structures between simulations can lead to mass ratios less than unity if the identity of the major progenitor swaps at very early times (merger 'A' in first two columns).</caption><thead><tr><th></th><th>Smallest z = 2 merger</th><th>Smaller z = 2 merger</th><th>Reference</th><th>Larger z = 2 merger</th><th>Largest z = 2 merger</th></tr></thead><tbody><tr><td>Target progenitor B</td><td>$z_{\text{infall}} = 1.94$</td><td>$z_{\text{infall}} = 2.03$</td><td>$z_{\text{infall}} = 1.99$</td><td>$z_{\text{infall}} = 1.99$</td><td>$z_{\text{infall}} = 2.24$</td></tr><tr><td>Infall $M_{200}$</td><td>$2.2 \times 10^{10} M_{\odot}$</td><td>$2.3 \times 10^{10} M_{\odot}$</td><td>$3.8 \times 10^{10} M_{\odot}$</td><td>$7.4 \times 10^{10} M_{\odot}$</td><td>$1.1 \times 10^{11} M_{\odot}$</td></tr><tr><td>Infall $M_{\star}$</td><td>$2.2 \times 10^{8} M_{\odot}$</td><td>$3.1 \times 10^{8} M_{\odot}$</td><td>$6.0 \times 10^{8} M_{\odot}$</td><td>$1.2 \times 10^{9} M_{\odot}$</td><td>$2.3 \times 10^{9} M_{\odot}$</td></tr><tr><td>Ratio $M_{200}$</td><td>1:10</td><td>1:9.8</td><td>1:6.0</td><td>1:2.9</td><td>1:2.1</td></tr><tr><td>Ratio $M_{\star}$</td><td>1:24</td><td>1:15</td><td>1:8.1</td><td>1:4.3</td><td>1:2.1</td></tr><tr><td>Infall $v_r/v_{\theta}$</td><td>9.5</td><td>9.1</td><td>16.</td><td>8.8</td><td>5.5</td></tr><tr><td>$z_{\text{pericentre}}$</td><td>1.70</td><td>1.76</td><td>1.79</td><td>1.80</td><td>1.98</td></tr><tr><td>$z_{\text{coalescence}}$</td><td>1.40</td><td>1.45</td><td>1.52</td><td>1.47</td><td>1.65</td></tr><tr><td>Earlier progenitor A</td><td>$z_{\text{infall}} = 3.11$</td><td>$z_{\text{infall}} = 3.17$</td><td>$z_{\text{infall}} = 3.35$</td><td>$z_{\text{infall}} = 3.35$</td><td>$z_{\text{infall}} = 3.49$</td></tr><tr><td>Infall $M_{200}$</td><td>$9.6 \times 10^{10} M_{\odot}$</td><td>$7.7 \times 10^{10} M_{\odot}$</td><td>$5.8 \times 10^{10} M_{\odot}$</td><td>$5.1 \times 10^{10} M_{\odot}$</td><td>$3.7 \times 10^{10} M_{\odot}$</td></tr><tr><td>Infall $M_{\star}$</td><td>$6.2 \times 10^{8} M_{\odot}$</td><td>$5.2 \times 10^{8} M_{\odot}$</td><td>$4.5 \times 10^{8} M_{\odot}$</td><td>$3.8 \times 10^{8} M_{\odot}$</td><td>$3.1 \times 10^{8} M_{\odot}$</td></tr><tr><td>Ratio $M_{200}$</td><td>1:0.7</td><td>1:0.9</td><td>1:1.3</td><td>1:1.4</td><td>1:1.9</td></tr><tr><td>Ratio $M_{\star}$</td><td>1:1.8</td><td>1:2.1</td><td>1:2.4</td><td>1:1.8</td><td>1:2.4</td></tr><tr><td>Infall $v_r/v_{\theta}$</td><td>9.5</td><td>4.1</td><td>2.3</td><td>2.3</td><td>2.0</td></tr><tr><td>$z_{\text{pericentre}}$</td><td>2.74</td><td>2.72</td><td>2.86</td><td>2.94</td><td>3.00</td></tr><tr><td>$z_{\text{coalescence}}$</td><td>2.17</td><td>2.30</td><td>2.34</td><td>2.34</td><td>2.38</td></tr><tr><td>Later progenitor C</td><td>$z_{\text{infall}} = 1.94$</td><td>$z_{\text{infall}} = 1.86$</td><td>$z_{\text{infall}} = 1.73$</td><td>$z_{\text{infall}} = 1.70$</td><td>$z_{\text{infall}} = 1.62$</td></tr><tr><td>Infall $M_{200}$</td><td>$1.5 \times 10^{10} M_{\odot}$</td><td>$1.8 \times 10^{10} M_{\odot}$</td><td>$1.8 \times 10^{10} M_{\odot}$</td><td>$1.8 \times 10^{10} M_{\odot}$</td><td>$1.7 \times 10^{10} M_{\odot}$</td></tr><tr><td>Infall $M_{\star}$</td><td>$1.7 \times 10^{8} M_{\odot}$</td><td>$2.8 \times 10^{8} M_{\odot}$</td><td>$2.3 \times 10^{8} M_{\odot}$</td><td>$2.8 \times 10^{8} M_{\odot}$</td><td>$2.2 \times 10^{8} M_{\odot}$</td></tr><tr><td>Ratio $M_{200}$</td><td>1:15</td><td>1:14</td><td>1:18</td><td>1:19</td><td>1:22</td></tr><tr><td>Ratio $M_{\star}$</td><td>1:31</td><td>1:20</td><td>1:29</td><td>1:26</td><td>1:50</td></tr><tr><td>Infall $v_r/v_{\theta}$</td><td>0.5</td><td>0.6</td><td>0.5</td><td>0.4</td><td>0.6</td></tr><tr><td>$z_{\text{pericentre}}$</td><td>1.61</td><td>1.59</td><td>1.48</td><td>1.41</td><td>1.35</td></tr><tr><td>$z_{\text{coalescence}}$</td><td>0.09</td><td>0.04</td><td>0.17</td><td>0.22</td><td>0.20</td></tr></tbody></table>

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, $z_{\text{pericentre}}$, and end of the interaction, $z_{\text{coalescence}}$, using the high-cadence simulation movies.

Our genetic modifications explicitly target merger 'B' (second row in Table 1; red in Figs 1 and 2), which in the reference case is the last major event until $z = 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 ($v_r/v_{\theta} \ge 5$). The range of merging stellar masses scanned by our suite ($2 \times 10^8 M_{\odot} \le M_{\star} \le 2 \times 10^9 M_{\odot}$) accurately brackets the proposed masses for a potential GSE progenitor (e.g. Belokurov et al. 2018; Helmi et al. 2018; Feuillet et al. 2020; Kruijssen et al. 2020; Mackereth & Bovy 2020; Naidu et al. 2020, 2021; Limberg 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 $z = 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 $\Lambda$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 Fig. 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, Pontzen & Peiris 2021; Pontzen et al. 2021). However, the differences in timings ($\approx 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 $z = 0$, but we note that their central galaxies also match in their final stellar mass ($M_{\star} = 1.8 \times 10^{10} M_{\odot}$) within 20 per cent of each other; Fig. 1). All final stellar masses are compatible at $1\sigma$ with empirical model predictions for this halo mass (Moster, Naab & White 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 $z \approx 2$, while fixing the large-
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Dynamical (top) and stellar (bottom) mass assemblies over cosmic time across the suite of genetically modified galaxies. The significance of a merger at $z \approx 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 $z=0$ (bottom) and evolve in similar large-scale environment (Fig. 2).</figcaption></figure>

scale cosmological environment, and the total dynamical and stellar mass budget of a Milky-Way-mass galaxy.

## 3 RESPONSE OF THE GALAXIES’ STRUCTURE TO GROWING AN EARLY MERGER

We start by quantifying the response of the $z=0$ stellar structure of each galaxy as we genetically modify their early assembly history. Fig. 3 shows UVI-mock images of the stellar light at $z=0$ (top row) and luminosity-weighted line-of-sight velocity maps (bottom row) as we incrementally increase the significance of the $z \approx 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 $z=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 ($\ge 27$ mag arcsec$^{-2}$) star-forming disc. Such structure is reminiscent of the decoupled internal kinematics observed around local elliptical galaxies (e.g. Efstathiou, Ellis & Carter 1982; Bender 1988; Krajnović et al. 2011, 2015; Johnston et al. 2018; Prichard, Vaughan & Davies 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 $z=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 (Figs 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 $z \approx 2$ merger (e.g. second, third, and fourth columns decrease and increase in turn). We verified that this 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 Fig. 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 $z \approx 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. 2009, 2021a). Furthermore, later aspects of the galaxy’s history ($z \le 1$) relevant to disc formation are also modified, notably the constructiveness of filamentary gas accretion (e.g. Dekel et al. 2020; Kretschmer, Agertz & Teyssier 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
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Dynamical (top) and stellar (bottom) mass assemblies over cosmic time across the suite of genetically modified galaxies. The significance of a merger at $z \approx 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 $z = 0$ (bottom) and evolve in similar large-scale environment (Fig. 2).</figcaption></figure>

scale cosmological environment, and the total dynamical and stellar mass budget of a Milky-Way-mass galaxy.

## 3 RESPONSE OF THE GALAXIES’ STRUCTURE TO GROWING AN EARLY MERGER

We start by quantifying the response of the $z = 0$ stellar structure of each galaxy as we genetically modify their early assembly history. Fig. 3 shows UVI-mock images of the stellar light at $z = 0$ (top row) and luminosity-weighted line-of-sight velocity maps (bottom row) as we incrementally increase the significance of the $z \approx 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 $z = 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 ($\ge$ 27 mag arcsec$^{-2}$) star-forming disc. Such structure is reminiscent of the decoupled internal kinematics observed around local elliptical galaxies (e.g. Efstathiou, Ellis & Carter 1982; Bender 1988; Krajnović et al. 2011, 2015; Johnston et al. 2018; Prichard, Vaughan & Davies 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 $z = 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 (Figs 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 $z \approx 2$ merger (e.g. second, third, and fourth columns decrease and increase in turn). We verified that this 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 Fig. 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 $z \approx 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. 2009, 2021a). Furthermore, later aspects of the galaxy’s history ($z \le 1$) relevant to disc formation are also modified, notably the constructiveness of filamentary gas accretion (e.g. Dekel et al. 2020; Kretschmer, Agertz & Teyssier 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
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Projected dark matter density maps at z = 3.5 showing the large-scale environment of each merger scenario. As we enhance the significance of the z ≈ 2 merger ('B', red) using genetic modifications, its progenitor becomes visibly larger (left to right, circles show r<sub>200</sub>), while the cosmological filamentary 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).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Response of the z = 0 stellar light (top) and kinematics (bottom) to modifying our galaxies’ early histories. As we grow the z ≈ 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 (Fig. 1) and in similar large-scale environment (Fig. 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 (Figs 4, 5), that has distinct cosmological origin but overlapping chemodynamics and ages in every case (Figs 6, 7).</figcaption></figure>

to a future study, and now focus on extracting the signatures of each merger scenario in phase-space.

# 4 CHEMODYNAMICAL SIGNATURES OF HIGH-REDSHIFT MERGER HISTORIES

## 4.1 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. Fig. 4 shows the fractional distributions of their specific vertical angular momentum, $L_z$, and orbital energies, $E_{\text{tot}}$. We define the vertical disc direction from the parallel to the stellar angular momentum in the inner 5 kpc, and compute $E_{\text{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_{\text{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 Fig. 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 z = 2 merger’ scenario (left-hand panel), we can first identify different galactic components in phase-space: (i) the bulge, with bound orbits (low $E_{\text{tot}}$) and pressure-supported kinematics (low $L_z$); (ii) the stellar disc, as a sequence spreading across a range of $E_{\text{tot}}$ and extending towards rotationally supported, high-$L_z$ orbits; and (iii) the fainter stellar halo surrounding the galaxy, towards higher $E_{\text{tot}}$ and small $L_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 Fig. 3. The disc sequence in phase space is the most extended in $L_z$ for the quietest merger history (left-hand panel), significantly smaller for intermediate scenarios
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Projected dark matter density maps at z = 3.5 showing the large-scale environment of each merger scenario. As we enhance the significance of the z ≈ 2 merger ('B', red) using genetic modifications, its progenitor becomes visibly larger (left to right, circles show r<sub>200</sub>), while the cosmological filamentary 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).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Response of the z = 0 stellar light (top) and kinematics (bottom) to modifying our galaxies’ early histories. As we grow the z ≈ 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 (Fig. 1) and in similar large-scale environment (Fig. 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 (Figs 4, 5), that has distinct cosmological origin but overlapping chemodynamics and ages in every case (Figs 6, 7).</figcaption></figure>

to a future study, and now focus on extracting the signatures of each merger scenario in phase-space.

# 4 CHEMODYNAMICAL SIGNATURES OF HIGH-REDSHIFT MERGER HISTORIES

## 4.1 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. Fig. 4 shows the fractional distributions of their specific vertical angular momentum, $L_z$, and orbital energies, $E_{\text{tot}}$. We define the vertical disc direction from the parallel to the stellar angular momentum in the inner 5 kpc, and compute $E_{\text{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_{\text{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 Fig. 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 z = 2 merger’ scenario (left-hand panel), we can first identify different galactic components in phase-space: (i) the bulge, with bound orbits (low $E_{\text{tot}}$) and pressure-supported kinematics (low $L_z$); (ii) the stellar disc, as a sequence spreading across a range of $E_{\text{tot}}$ and extending towards rotationally supported, high-$L_z$ orbits; and (iii) the fainter stellar halo surrounding the galaxy, towards higher $E_{\text{tot}}$ and small $L_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 Fig. 3. The disc sequence in phase space is the most extended in $L_z$ for the quietest merger history (left-hand panel), significantly smaller for intermediate scenarios
<figure><img src="image_6.png" /></figure>

**Figure 4.** Fractional distributions of specific energy and vertical angular momentum for all stars within 50 kpc at $z=0$, as we increase the significance of the early merger (left to right). Broad trends and specific galactic structures identified in Fig. 3 are recovered in phase space (annotated), in particular the increasingly smaller stellar discs with increasing merger ratio (e.g. lack of a high-$L_z$ sequence in the right-most panel). A population of inner halo stars (low $L_z$, $E_{\text{tot}}$ comparable to disc stars) is present in all merger scenarios, but becomes more prominent compared to the high-$L_z$ 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 (Fig. 5). Clustered structures in the outer halo (e.g. left-hand panel) map to dwarf galaxy satellites and recently disrupted debris.

<figure><img src="image_7.png" /></figure>

**Figure 5.** Azimuthal and radial velocity distribution of stars with $2 \le r \le 50$ kpc, as we increase the significance of the early merger (left to right). For well-ordered disc galaxies (first four columns), the velocity distribution is bimodal, with disc stars (positive $v_\phi$, orange boxes) linking into a radially biased population ($v_r/v_\phi \gg 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 $z \approx 2$. Further to these kinematic similarities, radially biased GSE-like stars exhibit similar age and [Fe/H] distributions in all cases (Figs 6 and 7).

(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, $L_z$-tail linking to the bulge, and a population overlapping in $L_z$ but at higher $E_{\text{tot}}$. Similarly, the decoupled kinematics for the largest merger (right-hand panel) are visible as a minimal prograde disc sequence for the inner stars (small $E_{\text{tot}}$), and retrograde, higher $E_{\text{tot}}$ orbits for the outer stars.

Beyond these already noted features, the phase-space distributions exhibit a significant population of inner halo stars, with pressure-supported orbits (low $L_z$) that share the same range of $E_{\text{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, Fig. 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 \le r \le 50$ kpc) and show their azimuthal and radial velocity fractional distributions in Fig. 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 misalignment 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 $z \approx 2$.

However, we recover this orbital structure in each of our galaxies with ordered rotation (first fourth columns), whether our merger at $z \approx 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 Fig. 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><img src="image_7.png" /></figure>

**Figure 4.** Fractional distributions of specific energy and vertical angular momentum for all stars within 50 kpc at $z=0$, as we increase the significance of the early merger (left to right). Broad trends and specific galactic structures identified in Fig. 3 are recovered in phase space (annotated), in particular the increasingly smaller stellar discs with increasing merger ratio (e.g. lack of a high-$L_z$ sequence in the right-most panel). A population of inner halo stars (low $L_z$, $E_{\text{tot}}$ comparable to disc stars) is present in all merger scenarios, but becomes more prominent compared to the high-$L_z$ 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 (Fig. 5). Clustered structures in the outer halo (e.g. left-hand panel) map to dwarf galaxy satellites and recently disrupted debris.

<figure><img src="image_8.png" /></figure>

**Figure 5.** Azimuthal and radial velocity distribution of stars with $2 \le r \le 50$ kpc, as we increase the significance of the early merger (left to right). For well-ordered disc galaxies (first four columns), the velocity distribution is bimodal, with disc stars (positive $v_\phi$, orange boxes) linking into a radially biased population ($v_r/v_\phi \gg 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 $z \approx 2$. Further to these kinematic similarities, radially biased GSE-like stars exhibit similar age and [Fe/H] distributions in all cases (Figs 6 and 7).

(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, $L_z$-tail linking to the bulge, and a population overlapping in $L_z$ but at higher $E_{\text{tot}}$. Similarly, the decoupled kinematics for the largest merger (right-hand panel) are visible as a minimal prograde disc sequence for the inner stars (small $E_{\text{tot}}$), and retrograde, higher $E_{\text{tot}}$ orbits for the outer stars.

Beyond these already noted features, the phase-space distributions exhibit a significant population of inner halo stars, with pressure-supported orbits (low $L_z$) that share the same range of $E_{\text{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, Fig. 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 \le r \le 50$ kpc) and show their azimuthal and radial velocity fractional distributions in Fig. 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 misalignment 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 $z \approx 2$.

However, we recover this orbital structure in each of our galaxies with ordered rotation (first fourth columns), whether our merger at $z \approx 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 Fig. 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
Fig. 5 (orange number). The fraction of radial-to-disc stars follows a systematic trend with $z \approx 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 $r \ge 1$ kpc and $r \ge 3$ kpc, or focusing only on halo-like orbits with $|L_z| \le 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 $z \approx 2$, and might be relatively common across the Milky-Way analogue population (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 $z=0$ and the mass ratio of a $z \approx 2$ merger questions the cosmological origin of these stars, and we now rewind each progenitor's history to establish their exact nature.

## 4.2 The origin of radially biased halo stars

### 4.2.1 Ages and metallicity distributions

We now focus on testing the cosmological origin of the GSE-like stars identified kinematically in Fig. 5 (purple box with $|v_\phi| \le 50$ km s$^{-1}$, $|v_r| \le 240$ km s$^{-1}$, and $2 \le r \le 50$ kpc). We start by plotting their star formation histories in Fig. 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 ($z \ge 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 $z \approx 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 ($\approx 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 ($t_{\text{birth}} \ge 6$ Gyr), reflecting a contamination from younger stars with overlapping kinematics at $z=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 Fig. 7 the [Fe/H] distributions of the GSE-like stars (purple). Again, all scenarios exhibit broad distributions of metallicities reminiscent of that observed in the Milky Way's (e.g. Feuillet et al. 2020; Naidu et al. 2021), that all peak around a similar median – $\langle [\text{Fe/H}] \rangle = -0.61, -0.61, -0.70, -0.59, -0.58$ from left to right panels, respectively. We thus conclude that the order-of-magnitude scan in merger mass ratio at $z \approx 2$ not only produces

kinematic similarities (Fig. 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 $z=0$ also exhibit singled-out enhancements compared to the background distribution. We mark at the top of Fig. 6 the length of the interaction with mergers ‘A’ and ‘B’ (see Table 1 for $z_{\text{infall}}$ and $z_{\text{coalescence}}$<fnref n="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 GSE-like population at $z=0$, indicating that a significant fraction of its stars are of in-situ origin (see also Grand et al. 2020; Khoperskov et al. 2022a; Orkney et al. 2022) and could elevate the galaxy-wide star formation rate $\approx 10$ Gyr ago (e.g. Alzate, Bruzual & Díaz-González 2021).

Furthermore, the long tails towards low ($[\text{Fe/H}] \le -1.0$) and high ($[\text{Fe/H}] \ge -0.5$) metallicities, and the spread in ages, imply that multiple populations of different origins are populating the $z=0$ GSE-like population. In fact, despite having the most massive merger at $z \approx 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 $z=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.

### 4.2.2 Multiple co-existing in-situ and ex-situ populations

To establish the respective origin of stars and the role of each merger event in assembling the $z=0$ GSE-like feature, we track each GSE-like 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 $r_{200}$ and $2r_{1/2}$, $v$ 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 $z_{\text{coalescence}}$. Stars in this subpopulation primarily form during the starburst driven by this interaction (recall the peak in star-formation in Fig. 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 ex-situ 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).

<footnote n="1">¹We omit merger ‘C’ for clarity, as, despite its early infall, its long inspiral which lasts until $z=0$ does not drive a starburst.</footnote>
Fig. 5 (orange number). The fraction of radial-to-disc stars follows a systematic trend with $z \approx 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 $r \ge 1$ kpc and $r \ge 3$ kpc, or focusing only on halo-like orbits with $|L_z| \le 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 $z \approx 2$, and might be relatively common across the Milky-Way analogue population (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 $z=0$ and the mass ratio of a $z \approx 2$ merger questions the cosmological origin of these stars, and we now rewind each progenitor's history to establish their exact nature.

## 4.2 The origin of radially biased halo stars

### 4.2.1 Ages and metallicity distributions

We now focus on testing the cosmological origin of the GSE-like stars identified kinematically in Fig. 5 (purple box with $|v_\phi| \le 50$ km s$^{-1}$, $|v_r| \le 240$ km s$^{-1}$, and $2 \le r \le 50$ kpc). We start by plotting their star formation histories in Fig. 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 ($z \ge 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 $z \approx 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 ($\approx 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 ($t_{\text{birth}} \ge 6$ Gyr), reflecting a contamination from younger stars with overlapping kinematics at $z=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 Fig. 7 the [Fe/H] distributions of the GSE-like stars (purple). Again, all scenarios exhibit broad distributions of metallicities reminiscent of that observed in the Milky Way's (e.g. Feuillet et al. 2020; Naidu et al. 2021), that all peak around a similar median – $\langle[\text{Fe/H}]\rangle = -0.61, -0.61, -0.70, -0.59, -0.58$ from left to right panels, respectively. We thus conclude that the order-of-magnitude scan in merger mass ratio at $z \approx 2$ not only produces

kinematic similarities (Fig. 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 $z=0$ also exhibit singled-out enhancements compared to the background distribution. We mark at the top of Fig. 6 the length of the interaction with mergers ‘A’ and ‘B’ (see Table 1 for $z_{\text{infall}}$ and $z_{\text{coalescence}}$<fnref n="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 GSE-like population at $z=0$, indicating that a significant fraction of its stars are of in-situ origin (see also Grand et al. 2020; Khoperskov et al. 2022a; Orkney et al. 2022) and could elevate the galaxy-wide star formation rate $\approx 10$ Gyr ago (e.g. Alzate, Bruzual & Díaz-González 2021).

Furthermore, the long tails towards low ([Fe/H] $\le -1.0$) and high ([Fe/H] $\ge -0.5$) metallicities, and the spread in ages, imply that multiple populations of different origins are populating the $z=0$ GSE-like population. In fact, despite having the most massive merger at $z \approx 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 $z=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.

### 4.2.2 Multiple co-existing in-situ and ex-situ populations

To establish the respective origin of stars and the role of each merger event in assembling the $z=0$ GSE-like feature, we track each GSE-like 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 $r_{200}$ and $2r_{1/2}$, $v$ 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 $z_{\text{coalescence}}$. Stars in this subpopulation primarily form during the starburst driven by this interaction (recall the peak in star-formation in Fig. 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 ex-situ 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).

<footnote n="1">¹We omit merger ‘C’ for clarity, as, despite its early infall, its long inspiral which lasts until $z=0$ does not drive a starburst.</footnote>
<figure><img src="image_8.png" /></figure>

**Figure 6.** Star formation histories of GSE-like stars identified in Fig. 5 for each merger scenario. Their ages show extended distributions, with similar medians (≈ 11 Gyr) across all scenarios consistent with an ancient population. Interactions with progenitors ‘A’ and ‘B’ (interaction timings marked at the top) map onto isolated bursts of formation of GSE-like stars, indicating that a significant fraction of this population originates from early gas-rich mergers. The lack of such correspondence for merger ‘B’ and its largest merger mass ratio (right-most panel) suggests however that the origin of stars on GSE-like orbits is evolving across merger scenarios (Fig. 7). A younger, metal-rich, [α/Fe]-reduced population contaminates the sample of ancient populations due to our broad kinematic cuts to define GSE-like stars (see Section 4.2.2 and Appendix A for further characterization).

<figure><img src="image_9.png" /></figure>

**Figure 7.** Iron metallicity distribution of the GSE-like stars in each merger scenario (left to right), broken down into subpopulations linked to early merger events (coloured lines, formed before interaction ‘A’, between interactions ‘A’ and ‘B’ and afterwards, in top, middle, and bottom rows, respectively). In all cases, the overall population (purple) shows extended [Fe/H] distributions, with similar medians despite the wide-scan in mass ratio for the GSE-like merger progenitor. The origin of the stars however shifts between merger scenarios. The accreted component of the GSE-like progenitor (‘B’, red, second row) grows with its mass ratio (left to right), but stars formed during the interaction with an earlier merger (‘A’; light blue, second row) always significantly contribute and have overlapping [Fe/H]. A small mass ratio at $z \approx 2$ (left-hand panels) is thus compensated by the growth of this earlier-formed population, resulting in similar chemodynamics and magnitude for the GSE feature at $z = 0$ at fixed galaxy stellar and dynamical mass. Earlier populations (top row) contribute little to the overall population but are significant for the metal-poor tail of the distribution. Later populations (bottom row) represent a significant fraction of the total, but could be easily filtered due to their higher metallicities and lower [α/Fe] (Fig. A2).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Star formation histories of GSE-like stars identified in Fig. 5 for each merger scenario. Their ages show extended distributions, with similar medians (≈ 11 Gyr) across all scenarios consistent with an ancient population. Interactions with progenitors ‘A’ and ‘B’ (interaction timings marked at the top) map onto isolated bursts of formation of GSE-like stars, indicating that a significant fraction of this population originates from early gas-rich mergers. The lack of such correspondence for merger ‘B’ and its largest merger mass ratio (right-most panel) suggests however that the origin of stars on GSE-like orbits is evolving across merger scenarios (Fig. 7). A younger, metal-rich, [α/Fe]-reduced population contaminates the sample of ancient populations due to our broad kinematic cuts to define GSE-like stars (see Section 4.2.2 and Appendix A for further characterization).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Iron metallicity distribution of the GSE-like stars in each merger scenario (left to right), broken down into subpopulations linked to early merger events (coloured lines, formed before interaction ‘A’, between interactions ‘A’ and ‘B’ and afterwards, in top, middle, and bottom rows, respectively). In all cases, the overall population (purple) shows extended [Fe/H] distributions, with similar medians despite the wide-scan in mass ratio for the GSE-like merger progenitor. The origin of the stars however shifts between merger scenarios. The accreted component of the GSE-like progenitor (‘B’, red, second row) grows with its mass ratio (left to right), but stars formed during the interaction with an earlier merger (‘A’; light blue, second row) always significantly contribute and have overlapping [Fe/H]. A small mass ratio at $z \approx 2$ (left-hand panels) is thus compensated by the growth of this earlier-formed population, resulting in similar chemodynamics and magnitude for the GSE feature at $z = 0$ at fixed galaxy stellar and dynamical mass. Earlier populations (top row) contribute little to the overall population but are significant for the metal-poor tail of the distribution. Later populations (bottom row) represent a significant fraction of the total, but could be easily filtered due to their higher metallicities and lower [α/Fe] (Fig. A2).</figcaption></figure>
**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 $z_{\text{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-[$\text{[Fe/H]}$] 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).

Fig. 7 shows the respective contribution of each subpopulation to the total [Fe/H] distribution of GSE-like stars (top). We further show each subpopulation's respective distributions in $E_{\text{tot}}$ and [α/Fe] in 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 $z=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 (Fig. 3). In fact, when the final galaxy showcases organized rotation at $z=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 $z \approx 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_{\text{tot}}$ and [α/Fe] and their 2D combinations.

Furthermore, these ambiguities strongly propagate to potential inferences of the merger mass ratio at $z \approx 2$. For the four disc galaxies (first four columns), diminishing the mass ratio at fixed total stellar and dynamical mass (recall Fig. 1) reduces the accreted contribution from this merger, but is compensated by growing the earlier *in-situ* population. This systematic variation in origin of the $z=0$ GSE-like stars still results in similar overall [Fe/H] distributions, combining

with their overlapping star formation histories (Fig. 6), kinematics (Fig. 5), [α/Fe], and $E_{\text{tot}}$ (Appendix A).

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 $z=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.

## 4.3 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 $z=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 ($|v_\phi| \le 50 \text{ km s}^{-1}$, $|v_r| \le 240 \text{ km s}^{-1}$, and $2 \le r \le 50$ comoving kpc). We then repeat the assignment performed in Section 4.2 to identify subpopulations and show their respective numbers over time in Fig. 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 $z=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. Fig. 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 Fig. 8 the timings of the 'A' and 'B' interactions ($z_{\text{infall}}$ to $z_c$ 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 $z \approx 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
**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 $z_{\text{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-[$\text{[Fe/H]}$] 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).

Fig. 7 shows the respective contribution of each subpopulation to the total [Fe/H] distribution of GSE-like stars (top). We further show each subpopulation's respective distributions in $E_{\text{tot}}$ and [α/Fe] in 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 $z=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 (Fig. 3). In fact, when the final galaxy showcases organized rotation at $z=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 $z \approx 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_{\text{tot}}$ and [α/Fe] and their 2D combinations.

Furthermore, these ambiguities strongly propagate to potential inferences of the merger mass ratio at $z \approx 2$. For the four disc galaxies (first four columns), diminishing the mass ratio at fixed total stellar and dynamical mass (recall Fig. 1) reduces the accreted contribution from this merger, but is compensated by growing the earlier *in-situ* population. This systematic variation in origin of the $z=0$ GSE-like stars still results in similar overall [Fe/H] distributions, combining

with their overlapping star formation histories (Fig. 6), kinematics (Fig. 5), [α/Fe], and $E_{\text{tot}}$ (Appendix A).

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 $z=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.

## 4.3 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 $z=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 ($|v_\phi| \le 50 \text{ km s}^{-1}$, $|v_r| \le 240 \text{ km s}^{-1}$, and $2 \le r \le 50$ comoving kpc). We then repeat the assignment performed in Section 4.2 to identify subpopulations and show their respective numbers over time in Fig. 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 $z=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. Fig. 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 Fig. 8 the timings of the 'A' and 'B' interactions ($z_{\text{infall}}$ to $z_c$ 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 $z \approx 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
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Build-up over time of each subpopulation within the GSE-like feature (colours matching Fig. 7, thick lines are 3-pixel Gaussian smoothings of full-data thin lines). The timing of high-redshift interactions ('A' and 'B' shown at the top) mark the build-up of their associated populations. The composition in ancient populations of the GSE-like feature observed at $z=0$ is already largely in place at $z=1$ (all colours but the younger, contaminating grey population). Larger mass ratios of merger 'B' (centre and right-hand panels) show tentative 'Splash'-like population, visible as an increased number of *in-situ* stars on radial orbits following the interaction (blue), but this remains small compared to time-step-to-time-step noise and later dynamical evolution.</figcaption></figure>

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 $\sqrt{2r^3/GM(<r)}$ at $2r_{1/2,v}$ 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 $z=1$ and $z=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 $z=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. Wu et al. 2022; van Donkelaar, Agertz & Renaud 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.

# 5 DISCUSSION

## 5.1 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 $z=0$ ($M_* \approx 2 \times 10^{10} M_\odot$) and lack some of its defining late-time features such as a bar (Fig. 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 $z=2$ merger' and 'Larger $z=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 ($\approx 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 $|b| > 15^\circ$ and distances $D < 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.

Fig. 9 shows the resulting distributions of the total GSE-like population (purple) and each of its subcomponent (coloured lines matching Fig. 6) in [Fe/H], [$\alpha$/Fe]$^2$, and $E_{tot}$ (top to bottom). With this tailored selection function, we recover the same trends and conclusions as with galaxy-wide samples: (i) despite a factor-three difference in their $z \approx 2$ merger's stellar and dynamical mass, both scenarios have similar total number of GSE-like stars (38,664 and 41 828 for left and right respectively); (ii) the radial population of halo stars, at both low and mid metallicities ($[Fe/H] \le -1.0$ and $-1.0 < [Fe/H] \le -0.5$) is a mixture of multiple *in-situ* and *ex-situ* components originating from several high-redshift interactions; (iii) these populations do not separate clearly in [$\alpha$/Fe] and $E_{tot}$, posing a challenge to the reconstruction of early merger properties from $z=0$ data.

<footnote>$^2$Using the same stellar evolution model as in Agertz et al. (2021), we also recover higher [$\alpha$/Fe] ratios than observed in the Milky Way (see their Section 4 for further discussion).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Build-up over time of each subpopulation within the GSE-like feature (colours matching Fig. 7, thick lines are 3-pixel Gaussian smoothings of full-data thin lines). The timing of high-redshift interactions ('A' and 'B' shown at the top) mark the build-up of their associated populations. The composition in ancient populations of the GSE-like feature observed at $z=0$ is already largely in place at $z=1$ (all colours but the younger, contaminating grey population). Larger mass ratios of merger 'B' (centre and right-hand panels) show tentative 'Splash'-like population, visible as an increased number of *in-situ* stars on radial orbits following the interaction (blue), but this remains small compared to time-step-to-time-step noise and later dynamical evolution.</figcaption></figure>

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 $\sqrt{2r^3/GM(<r)}$ at $2r_{1/2,v}$ 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 $z=1$ and $z=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 $z=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. Wu et al. 2022; van Donkelaar, Agertz & Renaud 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.

# 5 DISCUSSION

## 5.1 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 $z=0$ ($M_* \approx 2 \times 10^{10} M_\odot$) and lack some of its defining late-time features such as a bar (Fig. 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 $z=2$ merger' and 'Larger $z=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 ($\approx 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 $|b| > 15^\circ$ and distances $D < 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.

Fig. 9 shows the resulting distributions of the total GSE-like population (purple) and each of its subcomponent (coloured lines matching Fig. 6) in [Fe/H], [$\alpha$/Fe]$^2$, and $E_{\text{tot}}$ (top to bottom). With this tailored selection function, we recover the same trends and conclusions as with galaxy-wide samples: (i) despite a factor-three difference in their $z \approx 2$ merger's stellar and dynamical mass, both scenarios have similar total number of GSE-like stars (38,664 and 41 828 for left and right respectively); (ii) the radial population of halo stars, at both low and mid metallicities ($[\text{Fe/H}] \le -1.0$ and $-1.0 < [\text{Fe/H}] \le -0.5$) is a mixture of multiple *in-situ* and *ex-situ* components originating from several high-redshift interactions; (iii) these populations do not separate clearly in [$\alpha$/Fe] and $E_{\text{tot}}$, posing a challenge to the reconstruction of early merger properties from $z=0$ data.

<footnote>$^2$Using the same stellar evolution model as in Agertz et al. (2021), we also recover higher [$\alpha$/Fe] ratios than observed in the Milky Way (see their Section 4 for further discussion).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** [Fe/H], [α/Fe], and $E_{tot}$ distributions (top to bottom) of kinematically identified GSE-like stars in a mock Galactic halo survey from the Solar neighbourhood. 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 $z \approx 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 Fig. 7), which do not separate clearly in chemical or phase-space to unequivocally identify each scenario.</figcaption></figure>

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] $\le -1.0$ ([Fe/H] $\le -0.5$) in the mock solar neighbourhood does not vary significantly between the two scenarios, with 9677 stars and 10 343 for the left-hand 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 $z=0$ metallicities in a GSE-like feature to the mass ratio of a single GSE-progenitor 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.

## 5.2 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 $z=0$. In particular, Section 4.2 highlights the important contribution of stars formed during an earlier interaction at $z \approx 3$ (merger 'A'), which is difficult to separate from ex-situ stars brought with the radial event at $z \approx 2$ explicitly targeted by our genetic modifications (merger 'B'). In Section 2.2 and Fig. 1, we show that the mass growths of major progenitors across all merger scenarios are within the $2\sigma$ population from $z \approx 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 $\Lambda$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 $\Lambda$CDM universe. However, such likelihood is not readily available due to the large parameter space that needs to be sampled across all possible $\Lambda$CDM merger histories. Instead, we extract the summary statistics introduced by Fakhouri, Ma & Boylan-Kolchin (2010) from the Millenium simulation suite (Springel et al. 2005; Boylan-Kolchin et al. 2009) to describe the average population merger rate of dark matter haloes:

$$ \frac{dN}{d\xi dz}(M_{200}, \xi, z) = A \left( \frac{M_{200}}{10^{12} M_{\odot}} \right)^{\alpha} \xi^{\beta} e^{\left(\frac{\xi}{\xi_0}\right)^{\gamma}} (1+z)^{\eta}, \quad (1) $$

where $\xi$ is the dark matter merger mass ratio and $A = 0.0144$, $\xi_0 = 9.72 \times 10^{-3}$, $\alpha = 0.133$, $\beta = -1.995$, $\gamma = 0.263$, $\eta = 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 $z=5$ and $z=1$ along the median mass accretion history determined in Section 2.2, first for mass ratios greater than 1:5. We obtain $N(0.2 \le \xi \le 1) = 2.01$, which is consistent with our merger scenarios in which Milky-Way-mass progenitors experience two significant events between $z=5$ and $z=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; Fattahi et al. 2019; Mackereth 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 $\xi$, rapidly suppressing increasingly major mergers (e.g. Fakhouri & Ma 2008; Fakhouri et al. 2010). Re-integrating the merger rate with $\xi \ge 1:4$, 1:3, and 1:2 gives $N=1.67$, $N=1.28$, and $N=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 $z=5$ and $z=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 $\Lambda$CDM universe. Our genetic modifications generate minimal shifts in the likelihood of modified initial con-
<figure><img src="image_12.png" /><figcaption>**Figure 9.** [Fe/H], [α/Fe], and $E_{tot}$ distributions (top to bottom) of kinematically identified GSE-like stars in a mock Galactic halo survey from the Solar neighbourhood. 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 $z \approx 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 Fig. 7), which do not separate clearly in chemical or phase-space to unequivocally identify each scenario.</figcaption></figure>

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] $\le -1.0$ ([Fe/H] $\le -0.5$) in the mock solar neighbourhood does not vary significantly between the two scenarios, with 9677 stars and 10 343 for the left-hand 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 $z=0$ metallicities in a GSE-like feature to the mass ratio of a single GSE-progenitor 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.

## 5.2 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 $z=0$. In particular, Section 4.2 highlights the important contribution of stars formed during an earlier interaction at $z \approx 3$ (merger 'A'), which is difficult to separate from ex-situ stars brought with the radial event at $z \approx 2$ explicitly targeted by our genetic modifications (merger 'B'). In Section 2.2 and Fig. 1, we show that the mass growths of major progenitors across all merger scenarios are within the $2\sigma$ population from $z \approx 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 $\Lambda$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 $\Lambda$CDM universe. However, such likelihood is not readily available due to the large parameter space that needs to be sampled across all possible $\Lambda$CDM merger histories. Instead, we extract the summary statistics introduced by Fakhouri, Ma & Boylan-Kolchin (2010) from the Millenium simulation suite (Springel et al. 2005; Boylan-Kolchin et al. 2009) to describe the average population merger rate of dark matter haloes:

$$ \frac{dN}{d\xi dz}(M_{200}, \xi, z) = A \left( \frac{M_{200}}{10^{12} M_{\odot}} \right)^{\alpha} \xi^{\beta} e^{\left(\frac{\xi}{\xi_0}\right)^{\gamma}} (1+z)^{\eta}, \quad (1) $$

where $\xi$ is the dark matter merger mass ratio and $A = 0.0144$, $\xi_0 = 9.72 \times 10^{-3}$, $\alpha = 0.133$, $\beta = -1.995$, $\gamma = 0.263$, $\eta = 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 $z=5$ and $z=1$ along the median mass accretion history determined in Section 2.2, first for mass ratios greater than 1:5. We obtain $N(0.2 \le \xi \le 1) = 2.01$, which is consistent with our merger scenarios in which Milky-Way-mass progenitors experience two significant events between $z=5$ and $z=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; Fattahi et al. 2019; Mackereth 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 $\xi$, rapidly suppressing increasingly major mergers (e.g. Fakhouri & Ma 2008; Fakhouri et al. 2010). Re-integrating the merger rate with $\xi \ge 1:4$, 1:3, and 1:2 gives $N=1.67$, $N=1.28$, and $N=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 $z=5$ and $z=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 $\Lambda$CDM universe. Our genetic modifications generate minimal shifts in the likelihood of modified initial con-
ditions compared to their reference, ensuring that all of them are highly consistent with the ΛCDM power spectrum (Section 2.2 and Fig. 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 $N = 1.90$ for $\xi \ge 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).

## 6 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 $z \approx 2$ merger impacts the $z = 0$ chemodynamics of halo stars.

We start from a reference galaxy in a Milky-Way-mass dark matter halo ($M_{200} \approx 10^{12} M_{\odot}$) which undergoes a large merger at $z \approx 2$ ($M_{200}$ ratio of 1:6; $v_r/v_\theta = 16$). We then construct four alternative versions of this galaxy, systematically growing and decreasing the significance of this merger using the genetic modification approach (Roth et al. 2016; Rey & Pontzen 2018; Stopyra et al. 2021). The suite scans stellar mass ratios of 1:24, 1:15, 1:8, 1:4, and 1:2 for the interaction, covering the range of inferred values for a potential GSE-progenitor when interpreting Galactic data (e.g. Helmi et al. 2018; Kruijssen et al. 2020; Mackereth & Bovy 2020; Naidu et al. 2020; Feuillet et al. 2021; Naidu et al. 2021; Limberg et al. 2022).

All simulations are evolved to $z = 0$ with the extensive hydrodynamical, galaxy formation model of the VINTERGATAN project (Agertz et al. 2021), that can capture the multiphase structure of the galaxy's ISM ($\approx 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 $z = 0$ (Fig. 1) and a similar large-scale environment (Fig. 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 $z = 0$.

At fixed dynamical and stellar mass, modifications of the mass ratio at $z \approx 2$ result in a large diversity of galactic structure and morphology at $z = 0$ (Fig. 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 $z = 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 (Figs 4 and 5). We detect this kinematic GSE-like population with similar magnitude in all merger scenarios, whether we enhance

the $z \approx 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 (Fig. 6), [Fe/H] (Fig. 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 $z \approx 2$, cautioning studies that aim to reconstruct the properties of an assumed single progenitor from $z = 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 $z = 0$ have diverse origins, originating from a mixture of in-situ and ex-situ subpopulations with overlapping chemodynamical distributions (Fig. 7; see also Grand et al. 2020; Orkney et al. 2022).

In particular, we find that a significant fraction of the $z = 0$ feature originates from stars formed during an earlier, high-redshift interaction at $z \approx 3$, which have similar chemical abundances to the accreted population at $z \approx 2$. As we decrease the significance of the progenitor at $z \approx 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 $z = 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. 2021, 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 ($M_* \approx 2 \times 10^{10} M_{\odot}$) than our Galaxy ($\ge 6 \times 10^{10} M_{\odot}$), and blindly upscaling stellar masses across the merger tree would bring our GSE progenitors towards the massive end of commonly inferred values ($M_* \ge 10^9 M_{\odot}$), 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 $z = 0$ chemodynamical data is challenging. Furthermore, the numerical modelling of the stellar mass growth of galaxies, particularly that of dwarf galaxies at $z \approx 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 proto-galaxy. In addition
ditions compared to their reference, ensuring that all of them are highly consistent with the ΛCDM power spectrum (Section 2.2 and Fig. 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 $N = 1.90$ for $\xi \ge 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).

## 6 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 $z \approx 2$ merger impacts the $z = 0$ chemodynamics of halo stars.

We start from a reference galaxy in a Milky-Way-mass dark matter halo ($M_{200} \approx 10^{12} M_{\odot}$) which undergoes a large merger at $z \approx 2$ ($M_{200}$ ratio of 1:6; $v_r/v_\theta = 16$). We then construct four alternative versions of this galaxy, systematically growing and decreasing the significance of this merger using the genetic modification approach (Roth et al. 2016; Rey & Pontzen 2018; Stopyra et al. 2021). The suite scans stellar mass ratios of 1:24, 1:15, 1:8, 1:4, and 1:2 for the interaction, covering the range of inferred values for a potential GSE-progenitor when interpreting Galactic data (e.g. Helmi et al. 2018; Kruijssen et al. 2020; Mackereth & Bovy 2020; Naidu et al. 2020; Feuillet et al. 2021; Naidu et al. 2021; Limberg et al. 2022).

All simulations are evolved to $z = 0$ with the extensive hydrodynamical, galaxy formation model of the VINTERGATAN project (Agertz et al. 2021), that can capture the multiphase structure of the galaxy's ISM ($\approx 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 $z = 0$ (Fig. 1) and a similar large-scale environment (Fig. 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 $z = 0$.

At fixed dynamical and stellar mass, modifications of the mass ratio at $z \approx 2$ result in a large diversity of galactic structure and morphology at $z = 0$ (Fig. 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 $z = 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 (Figs 4 and 5). We detect this kinematic GSE-like population with similar magnitude in all merger scenarios, whether we enhance

the $z \approx 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 (Fig. 6), [Fe/H] (Fig. 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 $z \approx 2$, cautioning studies that aim to reconstruct the properties of an assumed single progenitor from $z = 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 $z = 0$ have diverse origins, originating from a mixture of in-situ and ex-situ subpopulations with overlapping chemodynamical distributions (Fig. 7; see also Grand et al. 2020; Orkney et al. 2022).

In particular, we find that a significant fraction of the $z = 0$ feature originates from stars formed during an earlier, high-redshift interaction at $z \approx 3$, which have similar chemical abundances to the accreted population at $z \approx 2$. As we decrease the significance of the progenitor at $z \approx 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 $z = 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. 2021, 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 ($M_* \approx 2 \times 10^{10} M_{\odot}$) than our Galaxy ($\ge 6 \times 10^{10} M_{\odot}$), and blindly upscaling stellar masses across the merger tree would bring our GSE progenitors towards the massive end of commonly inferred values ($M_* \ge 10^9 M_{\odot}$), 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 $z = 0$ chemodynamical data is challenging. Furthermore, the numerical modelling of the stellar mass growth of galaxies, particularly that of dwarf galaxies at $z \approx 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 proto-galaxy. 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 ($z < 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; Fig. 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, Belokurov & Evans 2011; Deason et al. 2014; Han et al. 2022), or the distribution of globular clusters (e.g. Kruijssen 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).

## ACKNOWLEDGEMENTS

We would like to thank the referee for a constructive review that improved the quality of the manuscript. MR would like to thank Vasily Belokurov, Alyson Brooks, Amina Helmi, Kathryn Johnston, Harley Katz and Sarah Pearson for insightful discussions and comments during the construction of this work. MR is supported by the Beecroft Fellowship funded by Adrian Beecroft. OA and FR acknowledge support from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (grant 2019–04659) and the Royal Physiographic Society in Lund. TS is supported by a CIERA Postdoctoral Fellowship. AP is supported by the Royal Society. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 818085 GMGalaxies. We acknowledge PRACE for awarding us access to Joliot-Curie at GENCI/CEA, France to perform the simulations presented in this work. Computations presented in this work were in part performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Tetralith supercomputer, part of the National Supercomputer Centre, Linköping University. Development work was carried out on facilities supported by the Research Capital Investment Fund (RCIF) provided by UKRI and partially funded by the UCL Cosmoparticle Initiative.

We thank the developers of GENETIC (Stopyra et al. 2021), JUPYTER (Ragan-Kelley et al. 2014), MATPLOTLIB (Hunter 2007), NUMPY (van der Walt, Colbert & Varoquaux 2011), PYNBODY (Pontzen et al. 2013), SCIPY (Virtanen et al. 2020), and TANGOS (Pontzen & Tremmel 2018) for providing open-source softwares used extensively in this work. The Astrophysics Data Service (ADS) and arXiv preprint repository were used in this work.

**Author contributions:** The main roles of the authors were, using the CRediT (Contribution Roles Taxonomy) system<fnref n="3" />:

MR: Conceptualization; Data curation; Formal analysis; Investigation; Writing – original draft. OA: Funding Acquisition; Investigation; Resources; Software; Writing – review and editing. TS: Conceptualization; Writing – review and editing. FR: Writing – review and editing. GJ: Data Curation; AP: Conceptualization; Writing – review and editing. NM: Writing – review and editing. DF: Writing – review and editing. JR: Writing – review and editing.

## DATA AVAILABILITY

The data underlying this article will be shared upon reasonable request to the corresponding author. The IllustrisTNG data used in this article is available at https://www.tng-project.org, and we will share derived data products upon reasonable requests.

## REFERENCES

Agertz O., Kravtsov A. V., 2015, ApJ, 804, 18

Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

Agertz O. et al., 2020, MNRAS, 491, 1656

Agertz O. et al., 2021, MNRAS, 503, 5826

Alzate J. A., Bruzual G., Dfaz-González D. J., 2021, MNRAS, 501, 302

Amarante J. A. S., e Silva L. B., Debattista V. P., Smith M. C., 2020, ApJ, 891, L30

Applebaum E., Brooks A. M., Christensen C. R., Munshi F., Quinn T. R., Shen S., Tremmel M., 2021, ApJ, 906, 96

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Aubert D., Teyssier R., 2010, ApJ, 724, 244

Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143

Belokurov V., Kravtsov A., 2022, MNRAS

Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611

Belokurov V., Sanders J. L., Fattahi A., Smith M. C., Deason A. J., Evans N. W., Grand R. J. J., 2020, MNRAS, 494, 3880

Bender R., 1988, A&A, 202, L5

Bignone L. A., Helmi A., Tissera P. B., 2019, ApJ, 883, L5

Bird J. C., Kazantzidis S., Weinberg D. H., Guedes J., Callegari S., Mayer L., Madau P., 2013, ApJ, 773, 43

Bird J. C., Loebman S. R., Weinberg D. H., Brooks A. M., Quinn T. R., Christensen C. R., 2021, MNRAS, 503, 1815

Bonaca A., Conroy C., Wetzel A., Hopkins P. F., Kereš D., 2017, ApJ, 845, 101

Bonaca A. et al., 2020, ApJ, 897, L18

Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., Lemson G., 2009, MNRAS, 398, 1150

Buck T., Obreja A., Macciò A. V., Minchev I., Dutton A. A., Ostriker J. P., 2020, MNRAS, 491, 1656

Buck T., Rybizki J., Buder S., Obreja A., Macciò A. V., Pfrommer C., Steinmetz M., Ness M., 2021, MNRAS, 508, 3365

Cadiou C., Pontzen A., Peiris H. V., 2021, MNRAS, 502, 5480

Cadiou C., Pontzen A., Peiris H. V., 2022, MNRAS, 517, 3459

Cappellari M. et al., 2011, MNRAS, 416, 1680

Chabrier G., 2003, PASP, 115, 763

Conroy C., Naidu R. P., Zaritsky D., Bonaca A., Cargile P., Johnson B. D., Caldwell N., 2019, ApJ, 887, 237

Cooper A. P., Parry O. H., Lowing B., Cole S., Frenk C., 2015, MNRAS, 454, 3185

Courty S., Alimi J. M., 2004, A&A, 416, 875

Davies J. J., Crain R. A., Pontzen A., 2021, MNRAS, 501, 236

<footnote n="3">³https://authorservices.wiley.com/author-resources/Journal-Authors/open-a-ccess/credit.html</footnote>
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 ($z < 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; Fig. 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, Belokurov & Evans 2011; Deason et al. 2014; Han et al. 2022), or the distribution of globular clusters (e.g. Kruijssen 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).

## ACKNOWLEDGEMENTS

We would like to thank the referee for a constructive review that improved the quality of the manuscript. MR would like to thank Vasily Belokurov, Alyson Brooks, Amina Helmi, Kathryn Johnston, Harley Katz and Sarah Pearson for insightful discussions and comments during the construction of this work. MR is supported by the Beecroft Fellowship funded by Adrian Beecroft. OA and FR acknowledge support from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (grant 2019–04659) and the Royal Physiographic Society in Lund. TS is supported by a CIERA Postdoctoral Fellowship. AP is supported by the Royal Society. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 818085 GMGalaxies. We acknowledge PRACE for awarding us access to Joliot-Curie at GENCI/CEA, France to perform the simulations presented in this work. Computations presented in this work were in part performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Tetralith supercomputer, part of the National Supercomputer Centre, Linköping University. Development work was carried out on facilities supported by the Research Capital Investment Fund (RCIF) provided by UKRI and partially funded by the UCL Cosmoparticle Initiative.

We thank the developers of GENETIC (Stopyra et al. 2021), JUPYTER (Ragan-Kelley et al. 2014), MATPLOTLIB (Hunter 2007), NUMPY (van der Walt, Colbert & Varoquaux 2011), PYNBODY (Pontzen et al. 2013), SCIPY (Virtanen et al. 2020), and TANGOS (Pontzen & Tremmel 2018) for providing open-source softwares used extensively in this work. The Astrophysics Data Service (ADS) and arXiv preprint repository were used in this work.

**Author contributions:** The main roles of the authors were, using the CRediT (Contribution Roles Taxonomy) system<fnref n="3" />:

MR: Conceptualization; Data curation; Formal analysis; Investigation; Writing – original draft. OA: Funding Acquisition; Investigation; Resources; Software; Writing – review and editing. TS: Conceptualization; Writing – review and editing. FR: Writing – review and editing. GJ: Data Curation; AP: Conceptualization; Writing – review and editing. NM: Writing – review and editing. DF: Writing – review and editing. JR: Writing – review and editing.

## DATA AVAILABILITY

The data underlying this article will be shared upon reasonable request to the corresponding author. The IllustrisTNG data used in this article is available at https://www.tng-project.org, and we will share derived data products upon reasonable requests.

## REFERENCES

Agertz O., Kravtsov A. V., 2015, ApJ, 804, 18

Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

Agertz O. et al., 2020, MNRAS, 491, 1656

Agertz O. et al., 2021, MNRAS, 503, 5826

Alzate J. A., Bruzual G., Dfaz-González D. J., 2021, MNRAS, 501, 302

Amarante J. A. S., e Silva L. B., Debattista V. P., Smith M. C., 2020, ApJ, 891, L30

Applebaum E., Brooks A. M., Christensen C. R., Munshi F., Quinn T. R., Shen S., Tremmel M., 2021, ApJ, 906, 96

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Aubert D., Teyssier R., 2010, ApJ, 724, 244

Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143

Belokurov V., Kravtsov A., 2022, MNRAS

Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611

Belokurov V., Sanders J. L., Fattahi A., Smith M. C., Deason A. J., Evans N. W., Grand R. J. J., 2020, MNRAS, 494, 3880

Bender R., 1988, A&A, 202, L5

Bignone L. A., Helmi A., Tissera P. B., 2019, ApJ, 883, L5

Bird J. C., Kazantzidis S., Weinberg D. H., Guedes J., Callegari S., Mayer L., Madau P., 2013, ApJ, 773, 43

Bird J. C., Loebman S. R., Weinberg D. H., Brooks A. M., Quinn T. R., Christensen C. R., 2021, MNRAS, 503, 1815

Bonaca A., Conroy C., Wetzel A., Hopkins P. F., Kereš D., 2017, ApJ, 845, 101

Bonaca A. et al., 2020, ApJ, 897, L18

Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., Lemson G., 2009, MNRAS, 398, 1150

Buck T., Obreja A., Macciò A. V., Minchev I., Dutton A. A., Ostriker J. P., 2020, MNRAS, 491, 1656

Buck T., Rybizki J., Buder S., Obreja A., Macciò A. V., Pfrommer C., Steinmetz M., Ness M., 2021, MNRAS, 508, 3365

Cadiou C., Pontzen A., Peiris H. V., 2021, MNRAS, 502, 5480

Cadiou C., Pontzen A., Peiris H. V., 2022, MNRAS, 517, 3459

Cappellari M. et al., 2011, MNRAS, 416, 1680

Chabrier G., 2003, PASP, 115, 763

Conroy C., Naidu R. P., Zaritsky D., Bonaca A., Cargile P., Johnson B. D., Caldwell N., 2019, ApJ, 887, 237

Cooper A. P., Parry O. H., Lowing B., Cole S., Frenk C., 2015, MNRAS, 454, 3185

Courty S., Alimi J. M., 2004, A&A, 416, 875

Davies J. J., Crain R. A., Pontzen A., 2021, MNRAS, 501, 236

<footnote n="3">³https://authorservices.wiley.com/author-resources/Journal-Authors/open-a-ccess/credit.html</footnote>
Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS, 416, 2903

Deason A. J., Belokurov V., Koposov S. E., Rockosi C. M., 2014, ApJ, 787, 30

Dekel A., Ginzburg O., Jiang F., Freundlich J., Lapiner S., Ceverino D., Primack J., 2020, MNRAS, 493, 4126

Di Matteo P., Combes F., Melchior A.-L., Semelin B., 2007, A&A, 468, 61

Di Matteo P., Haywood M., Lehnert M. D., Katz D., Khoperskov S., Snaith O. N., Gómez A., Robichon N., 2019, A&A, 632, A4

Dillamore A. M., Belokurov V., Font A. S., McCarthy I. G., 2022, MNRAS, 513, 1867

Dodd E., Callingham T. M., Helmi A., Matsuno T., Ruiz-Lara T., Balbinot E., Lovdal S., 2023, A&A, 670, L2

Donlon T., II, Newberg H. J., 2022, preprint(arXiv:2211.12576)

Donlon T., II, Newberg H. J., Sanderson R., Widrow L. M., 2020, ApJ, 902, 119

Donlon T., II, Newberg H. J., Kim B., Lepine S., 2022, ApJ, 932, L16

Duc P.-A. et al., 2015, MNRAS, 446, 120

Efstathiou G., Ellis R. S., Carter D., 1982, MNRAS, 201, 975

Elias L. M., Sales L. V., Helmi A., Hernquist L., 2020, MNRAS, 495, 29

Escala I. et al., 2018, MNRAS, 474, 2194

Fakhouri O., Ma C.-P., 2008, MNRAS, 386, 577

Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267

Fattahi A. et al., 2019, MNRAS, 484, 4471

Feuillet D. K., Bovy J., Holtzman J., Girardi L., MacDonald N., Majewski S. R., Nidever D. L., 2016, ApJ, 817, 40

Feuillet D. K., Feltzing S., Sahlholdt C., Casagrande L., 2020, MNRAS, 497, 109

Feuillet D. K., Sahlholdt C. L., Feltzing S., Casagrande L., 2021, MNRAS, 508, 1489

Feuillet D. K., Feltzing S., Sahlholdt C., Bensby T., 2022, ApJ, 934, 21

Font A. S. et al., 2020, MNRAS, 498, 1765

Gaia Collaboration 2016a, A&A, 595, A1

Gaia Collaboration 2016b, A&A, 595, A2

Gaia Collaboration 2018, A&A, 616, A1

Gaia Collaboration 2021, A&A, 649, A1

Gaia Collaboration 2022, A&A, 657, A7

Gill S. P. D., Knebe A., Gibson B. K., 2004, MNRAS, 351, 399

Girardi L. et al., 2010, ApJ, 724, 1030

Grand R. J. J. et al., 2017, MNRAS, 467, 179

Grand R. J. J. et al., 2018, MNRAS, 481, 1726

Grand R. J. J. et al., 2020, MNRAS, 497, 1603

Grand R. J. J. et al., 2021, MNRAS, 507, 4953

Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756

Gurvich A. B. et al., 2023, MNRAS, 519, 2598

Haardt F., Madau P., 1996, ApJ, 461, 20

Hafen Z. et al., 2022, MNRAS, 514, 5056

Han J. J. et al., 2022, AJ, 164, 249

Haywood M., Di Matteo P., Lehnert M. D., Snaith O., Khoperskov S., Gómez A., 2018, ApJ, 863, 113

Helmi A., 2020, ARA&A, 58, 205

Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85

Hernquist L., 1993, ApJ, 409, 548

Hopkins P. F., Cox T. J., Younger J. D., Hernquist L., 2009, ApJ, 691, 1168

House E. L. et al., 2011, MNRAS, 415, 2652

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Ishiyama T. et al., 2021, MNRAS, 506, 4210

Jackson R. A., Kaviraj S., Martin G., Devriendt J. E. G., Noakes-Kettel E. A., Silk J., Ogle P., Dubois Y., 2022, MNRAS, 511, 607

Johnston E. J., Hau G. K. T., Coccato L., Herrera C., 2018, MNRAS, 480, 3215

Khoperskov S. et al., 2022a, preprint(arXiv:2206.04521)

Khoperskov S. et al., 2022b, preprint(arXiv:2206.04522)

Khoperskov S. et al., 2022c, preprint(arXiv:2206.05491)

Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99

Kim J.-h. et al., 2014, ApJS, 210, 14

Knollmann S. R., Knebe A., 2009, ApJS, 182, 608

Koppelman H., Helmi A., Veljanoski J., 2018, ApJ, 860, L11

Krajnović D. et al., 2011, MNRAS, 414, 2923

Krajnović D. et al., 2015, MNRAS, 452, 2

Kretschmer M., Agertz O., Teyssier R., 2020, MNRAS, 497, 4346

Kruijssen J. M. D. et al., 2020, MNRAS, 498, 2472

Lim S. H., Raman K. A., Buckley M. R., Shih D., 2022, preprint(arXiv:2211.11765)

Limberg G., Souza S. O., Pérez-Villegas A., Rossi S., Perottoni H. D., Santucci R. M., 2022, ApJ, 935, 109

Ludlow A. D., Schaye J., Schaller M., Richings J., 2019, MNRAS, 488, L123

Ludlow A. D., Fall S. M., Schaye J., Obreschkow D., 2021, MNRAS, 508, 5114

Mackereth J. T., Bovy J., 2020, MNRAS, 492, 3631

Mackereth J. T. et al., 2019, MNRAS, 482, 3426

Majewski S. R. et al., 2017, AJ, 154, 94

Maksimova N. A., Garrison L. H., Eisenstein D. J., Hadzhiyska B., Bose S., Satterthwaite T. P., 2021, MNRAS, 508, 4017

Maoz D., Mannucci F., Brandt T. D., 2012, MNRAS, 426, 3282

Martell S. L. et al., 2017, MNRAS, 465, 3203

Martin G., Kaviraj S., Devriendt J. E. G., Dubois Y., Pichon C., 2018, MNRAS, 480, 2266

Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504

Matsuno T., Koppelman H. H., Helmi A., Aoki W., Ishigaki M. N., Suda T., Yuan Z., Hattori K., 2022, A&A, 661, A103

McMillan P. J., 2017, MNRAS, 465, 76

Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822

Myeong G. C., Vasiliev E., Iorio G., Evans N. W., Belokurov V., 2019, MNRAS, 488, 1235

Myeong G. C., Belokurov V., Aguado D. S., Evans N. W., Caldwell N., Bradley J., 2022, ApJ, 938, 21

Naab T., Ostriker J. P., 2017, ARA&A, 55, 59

Naab T. et al., 2014, MNRAS, 444, 3357

Naidu R. P., Conroy C., Bonaca A., Johnson B. D., Ting Y.-S., Caldwell N., Zaritsky D., Cargile P. A., 2020, ApJ, 901, 48

Naidu R. P. et al., 2021, ApJ, 923, 92

Negroponte J., White S. D. M., 1983, MNRAS, 205, 1009

Nelson D. et al., 2019, Comput. Astron. Cosmol., 6, 2

Nissen P. E., Schuster W. J., 2010, A&A, 511, L10

Orkney M. D. A., Laporte C. F. P., Grand R. J. J., Gómez F. A., van de Voort F., Marinacci F., Pakmor R., Springel V., 2022, MNRAS, 517, L138

Pagnini G., Di Matteo P., Khoperskov S., Mastrobuono-Battisti A., Haywood M., Renaud F., Combes F., 2022, preprint(arXiv:2210.04245)

Planck Collaboration XIII, 2016, A&A, 594, A13

Pontzen A., Tremmel M., 2018, ApJS, 237, 23

Pontzen A., Roškar R., Stinson G., Woods R., 2013, Astrophysics Source Code Library, record ascl:1305.002

Pontzen A., Tremmel M., Roth N., Peiris H. V., Saintonge A., Volonteri M., Quinn T., Governato F., 2017, MNRAS, 465, 547

Pontzen A., Rey M. P., Cadiou C., Agertz O., Teyssier R., Read J., Orkney M. D. A., 2021, MNRAS, 501, 1755

Poole G. B., Mutch S. J., Croton D. J., Wyithe S., 2017, MNRAS, 472, 3659

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Prichard L. J., Vaughan S. P., Davies R. L., 2019, MNRAS, 488, 1679

Quillen A. C., Minchev I., Bland-Hawthorn J., Haywood M., 2009, MNRAS, 397, 1599

Ragan-Kelley M., Perez F., Granger B., Kluyver T., Ivanov P., Frederic J., Bussonnier M., 2014, The Jupyter/Python Architecture: A Unified View of Computational Research, from Interactive Exploration to Communication and Publication, Fall Meeting 2014, abstract id. H44D-07. American Geophysical Union, Washington

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Renaud F., Boily C. M., Naab T., Theis C., 2009, ApJ, 706, 67

Renaud F., Agertz O., Read J. I., Ryde N., Andersson E. P., Bensby T., Rey M. P., Feuillet D. K., 2021a, MNRAS, 503, 5846

Renaud F., Agertz O., Andersson E. P., Read J. I., Ryde N., Bensby T., Rey M. P., Feuillet D. K., 2021b, MNRAS, 503, 5868

Rey M. P., Pontzen A., 2018, MNRAS, 474, 45

Rey M. P., Starkenburg T. K., 2022, MNRAS, 510, 4208
Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS, 416, 2903

Deason A. J., Belokurov V., Koposov S. E., Rockosi C. M., 2014, ApJ, 787, 30

Dekel A., Ginzburg O., Jiang F., Freundlich J., Lapiner S., Ceverino D., Primack J., 2020, MNRAS, 493, 4126

Di Matteo P., Combes F., Melchior A.-L., Semelin B., 2007, A&A, 468, 61

Di Matteo P., Haywood M., Lehnert M. D., Katz D., Khoperskov S., Snaith O. N., Gómez A., Robichon N., 2019, A&A, 632, A4

Dillamore A. M., Belokurov V., Font A. S., McCarthy I. G., 2022, MNRAS, 513, 1867

Dodd E., Callingham T. M., Helmi A., Matsuno T., Ruiz-Lara T., Balbinot E., Lovdal S., 2023, A&A, 670, L2

Donlon T., II, Newberg H. J., 2022, preprint(arXiv:2211.12576)

Donlon T., II, Newberg H. J., Sanderson R., Widrow L. M., 2020, ApJ, 902, 119

Donlon T., II, Newberg H. J., Kim B., Lepine S., 2022, ApJ, 932, L16

Duc P.-A. et al., 2015, MNRAS, 446, 120

Efstathiou G., Ellis R. S., Carter D., 1982, MNRAS, 201, 975

Elias L. M., Sales L. V., Helmi A., Hernquist L., 2020, MNRAS, 495, 29

Escala I. et al., 2018, MNRAS, 474, 2194

Fakhouri O., Ma C.-P., 2008, MNRAS, 386, 577

Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267

Fattahi A. et al., 2019, MNRAS, 484, 4471

Feuillet D. K., Bovy J., Holtzman J., Girardi L., MacDonald N., Majewski S. R., Nidever D. L., 2016, ApJ, 817, 40

Feuillet D. K., Feltzing S., Sahlholdt C., Casagrande L., 2020, MNRAS, 497, 109

Feuillet D. K., Sahlholdt C. L., Feltzing S., Casagrande L., 2021, MNRAS, 508, 1489

Feuillet D. K., Feltzing S., Sahlholdt C., Bensby T., 2022, ApJ, 934, 21

Font A. S. et al., 2020, MNRAS, 498, 1765

Gaia Collaboration 2016a, A&A, 595, A1

Gaia Collaboration 2016b, A&A, 595, A2

Gaia Collaboration 2018, A&A, 616, A1

Gaia Collaboration 2021, A&A, 649, A1

Gaia Collaboration 2022, A&A, 657, A7

Gill S. P. D., Knebe A., Gibson B. K., 2004, MNRAS, 351, 399

Girardi L. et al., 2010, ApJ, 724, 1030

Grand R. J. J. et al., 2017, MNRAS, 467, 179

Grand R. J. J. et al., 2018, MNRAS, 481, 1726

Grand R. J. J. et al., 2020, MNRAS, 497, 1603

Grand R. J. J. et al., 2021, MNRAS, 507, 4953

Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756

Gurvich A. B. et al., 2023, MNRAS, 519, 2598

Haardt F., Madau P., 1996, ApJ, 461, 20

Hafen Z. et al., 2022, MNRAS, 514, 5056

Han J. J. et al., 2022, AJ, 164, 249

Haywood M., Di Matteo P., Lehnert M. D., Snaith O., Khoperskov S., Gómez A., 2018, ApJ, 863, 113

Helmi A., 2020, ARA&A, 58, 205

Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85

Hernquist L., 1993, ApJ, 409, 548

Hopkins P. F., Cox T. J., Younger J. D., Hernquist L., 2009, ApJ, 691, 1168

House E. L. et al., 2011, MNRAS, 415, 2652

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Ishiyama T. et al., 2021, MNRAS, 506, 4210

Jackson R. A., Kaviraj S., Martin G., Devriendt J. E. G., Noakes-Kettel E. A., Silk J., Ogle P., Dubois Y., 2022, MNRAS, 511, 607

Johnston E. J., Hau G. K. T., Coccato L., Herrera C., 2018, MNRAS, 480, 3215

Khoperskov S. et al., 2022a, preprint(arXiv:2206.04521)

Khoperskov S. et al., 2022b, preprint(arXiv:2206.04522)

Khoperskov S. et al., 2022c, preprint(arXiv:2206.05491)

Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99

Kim J.-h. et al., 2014, ApJS, 210, 14

Knollmann S. R., Knebe A., 2009, ApJS, 182, 608

Koppelman H., Helmi A., Veljanoski J., 2018, ApJ, 860, L11

Krajnović D. et al., 2011, MNRAS, 414, 2923

Krajnović D. et al., 2015, MNRAS, 452, 2

Kretschmer M., Agertz O., Teyssier R., 2020, MNRAS, 497, 4346

Kruijssen J. M. D. et al., 2020, MNRAS, 498, 2472

Lim S. H., Raman K. A., Buckley M. R., Shih D., 2022, preprint(arXiv:2211.11765)

Limberg G., Souza S. O., Pérez-Villegas A., Rossi S., Perottoni H. D., Santucci R. M., 2022, ApJ, 935, 109

Ludlow A. D., Schaye J., Schaller M., Richings J., 2019, MNRAS, 488, L123

Ludlow A. D., Fall S. M., Schaye J., Obreschkow D., 2021, MNRAS, 508, 5114

Mackereth J. T., Bovy J., 2020, MNRAS, 492, 3631

Mackereth J. T. et al., 2019, MNRAS, 482, 3426

Majewski S. R. et al., 2017, AJ, 154, 94

Maksimova N. A., Garrison L. H., Eisenstein D. J., Hadzhiyska B., Bose S., Satterthwaite T. P., 2021, MNRAS, 508, 4017

Maoz D., Mannucci F., Brandt T. D., 2012, MNRAS, 426, 3282

Martell S. L. et al., 2017, MNRAS, 465, 3203

Martin G., Kaviraj S., Devriendt J. E. G., Dubois Y., Pichon C., 2018, MNRAS, 480, 2266

Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504

Matsuno T., Koppelman H. H., Helmi A., Aoki W., Ishigaki M. N., Suda T., Yuan Z., Hattori K., 2022, A&A, 661, A103

McMillan P. J., 2017, MNRAS, 465, 76

Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822

Myeong G. C., Vasiliev E., Iorio G., Evans N. W., Belokurov V., 2019, MNRAS, 488, 1235

Myeong G. C., Belokurov V., Aguado D. S., Evans N. W., Caldwell N., Bradley J., 2022, ApJ, 938, 21

Naab T., Ostriker J. P., 2017, ARA&A, 55, 59

Naab T. et al., 2014, MNRAS, 444, 3357

Naidu R. P., Conroy C., Bonaca A., Johnson B. D., Ting Y.-S., Caldwell N., Zaritsky D., Cargile P. A., 2020, ApJ, 901, 48

Naidu R. P. et al., 2021, ApJ, 923, 92

Negroponte J., White S. D. M., 1983, MNRAS, 205, 1009

Nelson D. et al., 2019, Comput. Astron. Cosmol., 6, 2

Nissen P. E., Schuster W. J., 2010, A&A, 511, L10

Orkney M. D. A., Laporte C. F. P., Grand R. J. J., Gómez F. A., van de Voort F., Marinacci F., Pakmor R., Springel V., 2022, MNRAS, 517, L138

Pagnini G., Di Matteo P., Khoperskov S., Mastrobuono-Battisti A., Haywood M., Renaud F., Combes F., 2022, preprint(arXiv:2210.04245)

Planck Collaboration XIII, 2016, A&A, 594, A13

Pontzen A., Tremmel M., 2018, ApJS, 237, 23

Pontzen A., Roškar R., Stinson G., Woods R., 2013, Astrophysics Source Code Library, record ascl:1305.002

Pontzen A., Tremmel M., Roth N., Peiris H. V., Saintonge A., Volonteri M., Quinn T., Governato F., 2017, MNRAS, 465, 547

Pontzen A., Rey M. P., Cadiou C., Agertz O., Teyssier R., Read J., Orkney M. D. A., 2021, MNRAS, 501, 1755

Poole G. B., Mutch S. J., Croton D. J., Wyithe S., 2017, MNRAS, 472, 3659

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Prichard L. J., Vaughan S. P., Davies R. L., 2019, MNRAS, 488, 1679

Quillen A. C., Minchev I., Bland-Hawthorn J., Haywood M., 2009, MNRAS, 397, 1599

Ragan-Kelley M., Perez F., Granger B., Kluyver T., Ivanov P., Frederic J., Bussonnier M., 2014, The Jupyter/Python Architecture: A Unified View of Computational Research, from Interactive Exploration to Communication and Publication, Fall Meeting 2014, abstract id. H44D-07. American Geophysical Union, Washington

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Renaud F., Boily C. M., Naab T., Theis C., 2009, ApJ, 706, 67

Renaud F., Agertz O., Read J. I., Ryde N., Andersson E. P., Bensby T., Rey M. P., Feuillet D. K., 2021a, MNRAS, 503, 5846

Renaud F., Agertz O., Andersson E. P., Read J. I., Ryde N., Bensby T., Rey M. P., Feuillet D. K., 2021b, MNRAS, 503, 5868

Rey M. P., Pontzen A., 2018, MNRAS, 474, 45

Rey M. P., Starkenburg T. K., 2022, MNRAS, 510, 4208