# The stellar halo in Local Group Hestia simulations I. The in-situ component and the effect of mergers

Sergey Khoperskov, Ivan Minchev, Noam Libeskind, Misha Haywood, Paola Di Matteo, Vasily Belokurov, Matthias Steinmetz, Facundo A. Gomez, Robert J.J. Grand, Alexander Knebe, et al.

► To cite this version:

Sergey Khoperskov, Ivan Minchev, Noam Libeskind, Misha Haywood, Paola Di Matteo, et al. The stellar halo in Local Group Hestia simulations I. The in-situ component and the effect of mergers. Astronomy and Astrophysics - A&A, 2023, 677, pp.A89. 10.1051/0004-6361/202244232. hal-03692902

HAL Id: hal-03692902

https://hal.science/hal-03692902

Submitted on 10 Nov 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The stellar halo in Local Group Hestia simulations

## I. The in situ component and the effect of mergers

Sergey Khoperskov¹,*, Ivan Minchev¹, Noam Libeskind¹,², Misha Haywood³, Paola Di Matteo³, Vasily Belokurov⁴,⁵, Matthias Steinmetz¹, Facundo A. Gomez⁶,⁷, Robert J. J. Grand⁸,⁹,¹⁰, Yehuda Hoffman¹¹, Alexander Knebe¹²,¹³,¹⁴, Jenny G. Sorce¹⁵,¹⁶,¹, Martin Spaare¹⁷,¹, Elmo Tempel¹⁸,¹⁹, and Mark Vogelsberger²⁰

¹ Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
e-mail: sergey.khoperskov@gmail.com

² Univ. Lyon, Univ. Claude Bernard Lyon 1, CNRS, IP2I Lyon/IN2P3, IMR 5822, 69622 Villeurbanne, France

³ GEPI, Observatoire de Paris, PSL Research University, CNRS, Place Jules Janssen, 92195 Meudon, France

⁴ Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK

⁵ Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

⁶ Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile

⁷ Departamento de Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena, Chile

⁸ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str 1, 85748 Garching, Germany

⁹ Instituto de Astrofísica de Canarias, Calle Váa Láctea s/n, 38205 La Laguna, Tenerife, Spain

¹⁰ Departamento de Astrofísica, Universidad de La Laguna, Av. del Astrofísico Francisco Sánchez s/n, 38206 La Laguna, Tenerife, Spain

¹¹ Racah Institute of Physics, Hebrew University, Jerusalem 91904, Israel

¹² Departamento de Física Teórica, Módulo 15, Facultad de Ciencias, Universidad Autónoma de Madrid, 28049 Madrid, Spain

¹³ Centro de Investigación Avanzada en Física Fundamental (CIAFF), Facultad de Ciencias, Universidad Autónoma de Madrid, 28049 Madrid, Spain

¹⁴ International Centre for Radio Astronomy Research, University of Western Australia, 35 Stirling Highway, Crawley, Western Australia 6009, Australia

¹⁵ Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, 59000 Lille, France

¹⁶ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France

¹⁷ Institut für Physik und Astronomie, Universität Potsdam, Campus Golm, Haus 28, Karl-Liebknecht Straße 24-25, 14476 Potsdam, Germany

¹⁸ Tartu Observatory, University of Tartu, Observatoorium 1, 61602 Tõravere, Estonia

¹⁹ Estonian Academy of Sciences, Kohtu 6, 10130 Tallinn, Estonia

²⁰ Department of Physics, Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

Received 9 June 2022 / Accepted 1 February 2023

## ABSTRACT

Theory suggests that mergers play an important role in shaping galactic discs and stellar haloes, which was observationally confirmed in the Milky Way (MW) thanks to Gaia data. In this work, aiming to probe the contribution of mergers to the in situ stellar halo formation, we analyse six M 31 and MW analogues from the HESTIA suite of cosmological hydrodynamical zoom-in simulations of the Local Group. We found that all the HESTIA galaxies experience between one to four mergers with stellar mass ratios between 0.2 and 1 relative to the host at the time of the merger. These significant mergers, with a single exception, happened 7–11 Gyr ago. The overall impact of the most massive mergers in HESTIA is clearly seen as a sharp increase in the orbital eccentricity (and a corresponding decrease in the rotational velocity $V_\phi$) of pre-existing disc stars of the main progenitor, thus nicely reproducing the Splash-, Plume-like feature that was discovered in the MW. We do find a correlation between mergers and close pericentric passages of massive satellites and bursts of the star formation in the in situ component. Massive mergers sharply increase the disc velocity dispersion of the in situ stars; however, the latest significant merger often heats up the disc up to the numbers when the contribution of the previous ones is less prominent in the age-velocity dispersion relation. In HESTIA galaxies, the in situ halo is an important component of the inner stellar halo where its fraction is about 30–40%, while in the outer parts it typically does not exceed ≈5% beyond 15 kpc from the galactic centre. The simulations suggest that this component of the stellar haloes continues to grow well after mergers conclude; however, the most significant contribution comes from stars that formed recently before the merger. The orbital analysis of the HESTIA galaxies suggests that wedges in $R_{max}-Z_{max}$ (apocentre – maximum height from the mid-plane) space are mainly populated by the stars born in between significant mergers.

**Key words.** galaxies: evolution – galaxies: halos – galaxies: kinematics and dynamics – galaxies: structure
# 1. Introduction

According to A cold dark matter (CDM) cosmology, galaxies acquire their mass via the continuous merging of smaller satellites (White & Frenk 1991; Springel et al. 2005; Schaye et al. 2015) and accretion of gas from large-scale filaments (Kereš et al. 2005; Dekel et al. 2009). Ancient episodes of gas accretion are hard to detect in local galaxies because of the dissipative nature of the ISM, while stellar merger remnants scattered across large spatial scales around galaxies, still could be seen as shells, streams, filaments and tidal tails (Malin & Carter 1983; Ibata et al. 2001; Belokurov et al. 2006; McConnachie et al. 2009; Martínez-Delgado et al. 2010; Atkinson et al. 2013; Duc et al. 2015). While being minor in terms of the mass budget at $z=0$ (Bell et al. 2008; Deason et al. 2011, 2013, 2019; Mackereth & Bovy 2020), individual accretion events can cause transformation of host galaxies in the past and by contributing to the inner stellar halo (Bell et al. 2008; Zolotov et al. 2009; Purcell et al. 2010; Pillepich et al. 2015; Cooper et al. 2015; Deason et al. 2016; D'Souza & Bell 2018a).

Various models suggest that classical bulges can be formed through mergers (see, e.g. Bournaud et al. 2005; Christensen et al. 2014; Pillepich et al. 2015). In less violent events, a number of different externally driven processes can result in the thickening of galactic discs, for example, heating of a pre-existing thin disc by minor mergers (Quinn et al. 1993; Villalobos & Helmi 2008; Minchev et al. 2015), accretion of disrupted satellites (Chiba & Beers 2000; Abadi et al. 2003), and formation of thick stellar components by gas-rich mergers (Brook et al. 2004, 2005). In theory, disentangling the impact of various types of mergers can be done by comparing the eccentricity distributions for both in situ and accreted stellar populations (see, e.g. Sales et al. 2009; Di Matteo et al. 2011; Qu et al. 2011). However, one could expect that several mechanisms can contribute to the observed orbital composition of stellar populations even in a single galaxy (see, e.g. Wilson et al. 2011; Mackereth et al. 2019), thus making it hard to connect the present-day characteristics with the past.

Generally speaking, redistribution of the in situ formed stars caused by different types of mergers and close passages of massive systems can also result in an inner stellar halo (Brook et al. 2004; Zolotov et al. 2009; Purcell et al. 2010; Font et al. 2011; House et al. 2011; McCarthy et al. 2012; Tissera et al. 2013; Cooper et al. 2015; Pillepich et al. 2015; Monachesi et al. 2016). In the Milky Way (MW) this issue became topical in light of the recent Gaia (ESA) data releases (Gaia Collaboration 2016, 2018). Using the Gaia DR1-TGAS sample of local stars, Bonaca et al. (2017) conclude that locally observed metal-rich ([Fe/H] > −1) stars on halo-like orbits have likely been born in situ. Combining the Gaia DR2 data with chemical abundances from the APOGEE spectroscopic survey Di Matteo et al. (2019) suggest that the vast majority of stars at few kiloparsecs from the Sun with halo-like kinematics have thick disc-like chemical abundances (see also Haywood et al. 2018; Helmi et al. 2018).

Since the age distribution of the heated population has the same age distribution as the Gaia-Sausage-Enceladus (GSE) merger (Gallart et al. 2019), it is now accepted that this event played a key role in the formation of the inner stellar halo of the MW (Haywood et al. 2018; Helmi et al. 2018; Belokurov et al. 2020), as well as the disc (Lu et al. 2022; Ciucă et al. 2023; Xiang & Rix 2022). At the same time, the effect of even earlier mergers on the dynamics of the disc is not yet quantified.

A more detailed analysis of the Gaia DR2 by Belokurov et al. (2020) allowed for the in situ halo (heated

disc) to be defined as the ‘Splash’ stars (earlier dubbed ‘Plume’ Di Matteo et al. 2019) that have little to no angular momentum with a substantial fraction of stars on retrograde orbits. This was first found by Schuster et al. (2012) at the solar vicinity, who analysed the accreted and in situ halo populations from Nissen & Schuster (2010). This latter work does indeed represent the original discovery of GSE in chemo-kinematic spaces and at a few hundred parsecs from the Sun. According to Belokurov et al. (2020), the Splash stars can be found around $[Fe/H] < -0.5$ with $V_\phi < 80$ km s$^{-1}$; however, they have a wide range of eccentricities (0.5–1, see also Mackereth et al. 2019), which makes it difficult to use them as a proxy for the merger parameters. It is important to note, that Amarante et al. (2020) have recently proposed a different model of the Splash-like feature formation, where a scattering of massive clumps at high redshift forms a metal-rich, low angular momentum population, without the need for a major merger.

Thanks to a new generation of cosmological hydrodynamic simulations the study of the accretion histories of simulated galaxies is now possible in great detail (see, e.g. Schaye et al. 2015; Wetzel et al. 2016; Sawala et al. 2016; Pillepich et al. 2018); see also the recent review by Vogelsberger et al. (2020). A few previous studies have investigated the effects of the mergers on the host in cosmological hydrodynamical simulations. Belokurov et al. (2020) analysed the Auriga galaxies, demonstrating a sharp change in the host stars’ kinematics right after the merger. In particular, the MW-like galaxies show a sharp decrease (down to zero) in the rotational velocity of stars, which existed in the disc before the massive merger. Using the same set of models, Grand et al. (2020) found evidence that gas-rich mergers heated the proto-disc of the simulated galaxies, scattering stars onto less circular orbits, such that their rotational velocity and metallicity positively correlate, thus contributing an additional component that connects the galactic thick disc to the inner stellar halo. Grand et al. (2020) also suggest that the counter-rotating fraction of proto-galaxy stars can be used to infer the stellar mass (ratio) of the GSE progenitor. Using a single simulation from EAGLE, Bignone et al. (2019) demonstrate that a GSE-like merger can heat up stars in the early thin disc, forming the present-day thick disc, while the merger debris constitute the stellar halo. Of course, thick discs are much more complicated than what a single event can produce (e.g., the observed strong age gradient in the MW thick disc; Martig et al. 2016) where the outer parts come from the nested flares of mono-age populations (Minchev et al. 2015, 2019; García de la Cruz et al. 2021).

Renaud et al. (2021) present a new zoom-in simulation of a MW-like galaxy formation where the last major GSE-like merger deposited stars with high metallicity and non-zero tangential velocities. This component strongly overlaps with the in situ population in a number of parameters, making its identification very uncertain (see also Jean-Baptiste et al. 2017; Pagnini et al. 2023). Using the NIHAO-UHD simulations, Sestito et al. (2021) found the signature of a low-metallicity in situ population (with both prograde and retrograde rotation) with disc-like kinematics (see also Abadi et al. 2006; Sestito et al. 2019; Di Matteo et al. 2020). Sestito et al. (2021) also found that some of the retrograde stars can be formed in situ or have been deposited at very early times. More recently, Dillamore et al. (2022) have shown that about one-third of galaxies from the ARTEMIS cosmological simulations contain GSE-like features which comprise stellar debris from the most massive accreted satellite. Dillamore et al. (2022) also show that the mergers result in a rapid rotation of the disc, seen as the change of the
spin orientation – which tend to be aligned with the merger plane.

The aim of the paper is to investigate the impact of ancient
mergers on the in situ disc populations in a set of HESTIA con-
strained simulations of the Local Group (LG; Libeskind et al.
2020)<fnref n="1" />. In particular, the HESTIA simulations resemble both
realistic M 31 and MW galaxies in terms of their halo mass,
stellar disc mass, morphology separation, relative velocity, rota-
tion curves, bulge-disc morphology, satellite galaxy stellar mass
function, satellite radial distribution and the presence of a
Magellanic-Cloud-like objects (Libeskind et al. 2020) thus mak-
ing HESTIA simulations the best tool for studying both environ-
mental effects and the assembly history relevant for the MW and
M 31 galaxies.

In a series of works based on a new set of HESTIA high-resolution cosmological simulations of the LG galaxies we investigate the phase-space evolution and the present-day structure of the merger debris (Khoperskov et al. 2023a, hereafter Paper II) and the chemical abundance patterns as a function of stellar ages and kinematics of both accreted and in situ stellar populations (Khoperskov et al. 2023b, hereafter Paper III). In this paper, using six M 31 and MW analogues from the high-resolution hydrodynamical HESTIA simulations, we focus on how much the mergers contribute to the disc heating, its orbital transformation and the emerging of the Splash/Plume-like stellar populations discovered in the MW. The paper is structured as follows. In Sect. 2 we describe the HESTIA cosmological simulations – the initial conditions, physical model, snapshots analysis, and the definition of accreted and in situ stellar populations. In Sect. 3 we describe some parameters of mergers. In Sect. 4 we analyze the impact of the mergers on the in situ populations focussing on the kinematics and orbital composition of stars. In Sect. 5 we discuss the present-day properties and formation paths of the kinematically defined stellar haloes. Finally, in Sect. 6 we summarise our main results.

2. HESTIA simulations

2.1. HESTIA simulations: Initial conditions

In this work, we analyse the three highest resolution HESTIA
simulations of the LG. Each simulation is tailored to reproduce
a number of the LG properties (Libeskind et al. 2020), including
the massive disc galaxies resembling the MW and Andromeda
analogues with the population of smaller satellites at z = 0.
Throughout the paper, we analyse six galaxies where realis-
tic present-day properties of the satellite population provide the
opportunity to study the impact of the M 31- and MW-like accre-
tion history on the in situ formed stellar populations. Next we
provide a brief description of the code, physical model, and post-
processing of the output. For more details, we refer the reader
to HESTIA simulations’ introductory paper (Libeskind et al.
2020).

The HESTIA simulations were performed by using the
AREPO code (Springel 2005; Pakmor et al. 2016), which solves
the ideal magnetohydrodynamics equations on an unstructured
Voronoi mesh with a second order finite volume scheme. Self-
gravity and other source terms were coupled to the ideal
magenta-hydrodynamic equations by operator splitting. Gravi-
tational forces were computed using a hybrid TreePM technique
(Springel 2005) with two Fourier mesh levels, one for the full
box and one centred on the high-resolution region.

The initial conditions for the run are constrained using
Cosmicflows-2 (CF2) peculiar velocities (Tully et al. 2013). The
CF2 catalogue is first grouped (Sorce & Tempel 2018) and
a bias-minimization technique is applied (Sorce 2015). The
technique is based on the Wiener Filter/Constrained Realiza-
tion (Hoffman & Ribak 1991) algorithm, combined with the
reverse Zeldovich approximation, spelled out in Doumler et al.
(2013c,b,a).

The HESTIA simulations use the galaxy formation model
from Grand et al. (2017), which is based on the Illustris model
(Vogelsberger et al. 2013), and implements the most important
physical processes relevant for the formation and evolution of
galaxies. It includes cooling of gas via primordial and metal
cooling (Vogelsberger et al. 2013) and a spatially uniform UV
background (Vogelsberger et al. 2013). The ISM is described by
a subgrid model for a two-phase medium in which cold star-
forming clouds are embedded in a hot volume-filling medium
(Springel & Hernquist 2003). Gas that is denser than n<sub>thres</sub> =
0.13 cm<sup>-3</sup> forms stars following a Schmidt-type star formation
law. Star formation itself is done stochastically and creates star
particles with the target gas mass that represent single stellar
populations. The model includes mass loss and metal return
from asymptotic giant branch stars, core-collapse supernovae,
and Type Ia supernovae that are distributed in the cells around
a star particle. Galactic winds are implemented by creating a
wind with a given velocity and mass loading just outside the star-
forming phase (Grand et al. 2017, 2019; Nelson et al. 2019). The
Auriga model also follows the formation and growth of super-
massive black holes and includes their feedback as active galac-
tic nuclei.

The highest resolutions HESTIA simulations are based on
the re-runs of the low resolution dark matter (DM) only sim-
ulations where two overlapping 3.7 Mpc (2.5 Mpc h⁻<fnref n="1" />) spheres
are drawn around the two main LG members at z = 0 and then
populated with 8192³ effective particles. The mass and spatial
resolution achieved is m<sub>dm</sub> = 1.5 × 10⁵ M<sub>☉</sub>, m<sub>gas</sub> = 2.2 × 10⁴ M<sub>☉</sub>
and ε = 220 pc. HESTIA simulations assume a cosmology con-
sistent with the best fit values (Planck Collaboration XVI 2014):
σ<sub>8</sub> = 0.83 and H<sub>0</sub> = 100 h km s⁻<fnref n="1" /> Mpc⁻<fnref n="1" /> where h = 0.677. We
adopt Ω<sub>Λ</sub> = 0.682 throughout and Ω<sub>M</sub> = 0.270 and Ω<sub>b</sub> = 0.048.

2.2. HESTIA simulations: Post-processing

Halos and sub-haloes are identified at each redshift by using
the publicly available AHF<fnref n="2" /> halo finder (Knollmann & Knebe
2009). Galaxy and halo histories are estimated via merger trees
(MergerTree tool from the AHF package). At each redshift,
accretion events are found by identifying which subhaloes at
a given snapshot are identified as ‘field’ haloes at the previous
snapshot. Haloes may also grow via smooth accretion from the
environment, namely via gravitationally attracting particles in
their vicinity.

For each M 31- and MW-simulated galaxy, as well as each
snapshot, we define a coordinate system (x, y, z) centred on 10%
of the most bound in situ star particles and aligned with the prin-
cipal axes of this in situ stellar component. The disc plane of the
host galaxy thus lies in the xy plane, R<sub>gal</sub> is the cylindrical galac-
tocentric distance, and the rotation is along the z axis. Our study
uses velocity in galactocentric cylindrical coordinates: tangential
velocity (V<sub>φ</sub>), radial velocity (V<sub>r</sub>), and velocity in the z direction
(V<sub>z</sub>). We also make use of the integrals of motion, focussing on

<footnote n="1">¹ https://hestia.aip.de</footnote>

<footnote n="2">² http://popia.ft.uam.es/AHF/index.html</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Face-on (top) and edge-on (bottom) stellar density maps ($M_{\odot} pc^{-2}$) of the M 31 and MW analogues in the HESTIA simulations. Disc components of the simulated galaxies have different extensions; however, all the galaxies show a presence of isolated stellar overdensities – dwarf galaxies orbiting around the host galaxies – and many smooth streams, shells, and tidal features contributing to the smooth halo component made of stars that formed ex situ and that accreted at different epochs. Mock HST images, gas distributions and a visualisation of the HESTIA Local Group are presented in Libeskind et al. (2020).</figcaption></figure>

angular momentum in the z direction ($L_z$) and total orbit energy per unit mass $E$.

To characterise the orbital parameters of stars (eccentricity, maximum cylindrical distance or apocentres $R_{max}$, and maximal vertical excursion from the disc mid-plane $Z_{max}$), we do not analyse the motion of particles across the snapshots, but integrated the orbits of star particles in a smooth potential. This allowed us to obtain the instantaneous values of the orbital parameters for each snapshot independently, which is not possible from the direct output data. To compute the instantaneous orbital parameters of star particles, first we used AGAMA (Vasiliev 2019) to model a smooth gravitational potential. In order to avoid the perturbation of orbits from massive satellites, we interpolated the galaxy potential using only particles associated with the host galaxy. The potential due to dark matter and halo gas is represented by a symmetric expansion in spherical harmonics up to $l = 4$, while the potential of the stars and the gaseous disc was approximated by an azimuthal harmonic expansion up to $m = 4$. We subsequently integrated the orbits of star particles in this potential for 20 Gyr. This timescale was chosen to account for halo particles with small orbital frequencies. The orbits were integrated with AGAMA, using a Runge-Kutta DOP853 integrator with an adaptive time step.

# 3. HESTIA galaxies: Assembly history

## 3.1. In situ and accreted populations

In Fig. 1 we show both face-on and edge-on stellar density maps for all the M 31 and MW galaxies from the three highest-resolution HESTIA simulations (09–17, 17–11 and 37–11 from Libeskind et al. 2020). One can see that all the galaxies exhibit discy stellar components of various sizes surrounded by the smooth, but not featureless, stellar haloes. At the same time, the stellar haloes’ morphology (see bottom panels) reveals a number of tidal streams, shells, and a number of isolated dwarf galaxies, which all together somehow represent the assembly histories of the galaxies.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Formation histories of the HESTIA simulations: formation distance of stars as a function of lookback time, colour-coded by the stellar mass (in $log_{10}$ units). The diagrams include all the in situ stars, merger remnants, and satellites that survived until $z = 0$. The HESTIA galaxies show a diverse assembly history with substantial early accretion and a different number of massive dwarfs in the halo at the present epoch.</figcaption></figure>

The stellar content formation history for the six M 31 and MW HESTIA galaxies is presented in Fig. 2, where we show the amount of stellar mass that formed at different distances from the host galaxy over time (see similar results for the Latte suite of FIRE-2 simulations in Cunningham et al. 2022). One could trace a number of narrow tracks, representing the star formation inside dwarf galaxies that either merged at some point with the host or survived to the present day in the halo. At early times, the
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Orbital decay for merged dwarf galaxies (left) and satellites that survived (right) for each M31 and MW analogue in the HESTIA simulations. The satellites that survived are the dwarf galaxies inside the virial radius of the main progenitor, which are still gravitationally bound to their own sub-halo at $z=0$. Mergers correspond to the galaxies that, at some point in the evolution, become unbound and were disrupted inside the main progenitor. We note that some of these galaxies become unbound well before their mass centre reaches the centre of the main progenitor.</figcaption></figure>

assembly histories of HESTIA galaxies are quite complicated due to a number of mergers barely separated in time with the star formation distributed over a large radial range. We note that most stars (≈85% in the MW 37–11 galaxy and 65–68% in other galaxies) are still formed inside the main progenitor, which can be seen as a broad stripe in the bottom of each panel.

In order to quantify the impact of mergers on the host galaxies, we disentangled ‘in situ’ and ‘accreted’ stellar populations. We consider ‘in situ’ stars as those that were formed gravitationally bound to the main progenitor (either the MW or M31 analogue), while the rest (formed bound to other subhaloes) we mark as ‘accreted’. Throughout the paper, we define a merger event as the accretion of a dwarf galaxy that becomes gravitationally unbound from its own halo and gets bound to the main progenitor, according to the MergerTree tool from AHF. Therefore, all particles associated with a sub-halo right before this event are considered as accreted from a single merger and the last snapshot is used as the time of the merger. According to our definition, the mergers remnants constitute the stellar halo. Meanwhile, all the smaller systems that are still bound at $z=0$, while being inside the virial radius of the main progenitor, represent the population of survived dwarf galaxies of M31 and

MW analogues (same definition as in Fattahi et al. 2020). In other words, the in situ stars form within the main progenitor branch of the galaxy merger tree, independently of the origin of the star-forming gas. Ex situ stars form outside the main progenitor branch and are subsequently accreted onto the host galaxy through mergers or stripping events.

Since we aim to understand the impact of the GSE-like and even more ancient mergers on the MW-like discs, we focus our analysis on the accretion events, which result in the disruption of dwarf galaxies, while the impact of satellites orbiting around – similar to the Sgr dwarf and LMC and SMC systems – is beyond the scope of the present work. In Fig. 3 we show the orbital decay of dwarf galaxies (with the stellar mass $>10^6 M_\odot$ at $z=0$) as a function of time for all six simulated galaxies. We split the populations according to their fate in the potential of the main progenitor: the left column shows satellites that have been accreted (mergers) and have become unbound and disrupted, while on the right we present the ones that survived until $z=0$.

In Fig. 3 (left), we can see a number of mergers (≈10–40) in different galaxies, whose stellar remnants contribute to the diffuse halo component (see Fig. 1). We note that not all of the tracks in Fig. 3 (left) end up in the centre of the main
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Mergers history of the HESTIA galaxies. The in situ stellar mass growth for the six M31 and MW HESTIA galaxies is shown by the black line. Mergers are shown by the coloured circles representing the time of accretion and the stellar mass of a satellite at that time. The circles are colour-coded by the stellar mass ratio relative to the host (μ*, colour bars are on the right). Red horizontal lines show the populations of dwarf galaxies that survive until z = 0. Five of the most significant mergers are marked as M1–M5 where the number increases from the earliest to the most recent ones. The HESTIA galaxies experienced between one to four massive merger events with the stellar mass ratio >10% relative to the host in situ component, while the other mergers are less significant in terms of stellar mass. The number of merger events significantly decreases starting from ≈8 Gyr ago.</figcaption></figure>

progenitor. This means that the satellites become unbound some
time before their core reaches the centre or they are fully dis-
rupted in a shell-like structure without a prominent overdensity
close to the galactic centre at the time of the merger. One can see
that mergers mainly happen at earlier times because apparently
these dwarf galaxies, on average, were closer (or more generally
have smaller total energy) to the host at high redshift on average;
whereas, while dwarf galaxies, especially at their first pericentric
passage close to z = 0, arrive from larger distances (see, also,
Simon 2018; Hammer et al. 2021; Correa Magnus & Vasiliev
2022; Battaglia et al. 2022; Dupuy et al. 2022).

3.2. Mergers statistics

In Fig. 4 we show the stellar mass growth of the main progenitors (in situ stars, black curves), where the final stellar masses of our six HESTIA galaxies cover the range (4–8) × 10¹⁰ M☉, which is comparable to recent estimates of the MW mass (see, e.g. Bland-Hawthorn & Gerhard 2016). All the mergers with Mstars > 10⁶ M☉ are marked by the coloured circles. The colour of the circles corresponds to the stellar mass ratio relative to the main progenitor stellar mass (μ*) at the time of merger. The satellites that survived (bound (sub)haloes, see Fig. 3, right) are shown by red lines near lookback time 0. In this figure we also

mark the most significant five mergers as M1–M5 from the ear-
liest to the latest one. We define the most significant mergers
according to the relative stellar mass ratio at the time of the
merger. The masses of satellite galaxies and stellar merger debris
were calculated by using all star particles associated with a given
object by the AHF halo finder at the present day or at the time
of the infall, respectively. We note that we call mergers signifi-
cant because not all of them can be classified as ‘major’ merg-
ers (and sometimes not even as massive), especially relative to
the present-day mass of galaxies. Nevertheless, using our selec-
tion of HESTIA galaxies allowed us to explore the impact of
5×10⁸−2×10¹⁰ M☉ mergers on the main progenitors’ discs. The
phase-space-chemistry relations of the merger remnants will be
analysed in detail in subsequent works (Paper II and Paper III).

In Fig. 4 one can see that the total number of mergers varies
in the range of ≈10−40 for different galaxies; however, the num-
ber of significant mergers with μ* > 0.2 (i.e. at least 1:5 merg-
ers) is only 1−4. For most of the galaxies the last significant
merger happened >8 Gyr ago, similar to our current expecta-
tions for the MW (see Helmi 2020, for review and references
in the introduction), with the single exception of the M31 ana-
logue in the 09−18 simulation which, however, is in perfect
agreement with the M31 observational data suggesting a major
merger (≈2.5×10¹⁰ stellar mass) event 2.5−4.5 Gyr ago (see, e.g.
<figure><img src="image_6.png" /></figure>

**Fig. 5.** Time evolution of the short-to-long-axis ratio and $V_\phi/\sigma$ relation of the in situ stellar populations. Top: ratio between the mean height and the mean radial extension of the stellar component, calculated as the mean absolute values of the vertical, $\langle|z|\rangle$, and galactocentric radial, $\langle R_{\text{gal}} \rangle$, position of stars of a given age. Bottom: ratio between the mean rotational velocity ($V_\phi$) and the radial velocity dispersion of stars in a given age range. Stars of different ages are shown by different colours, specified by the colour bar on the right, where each line starts from the oldest stars formed in a 2 Gyr age-interval. The grey vertical lines highlight the five most significant mergers. The line length corresponds to the stellar mass ratio of the merger relative to the host at the time of accretion. The black dashed line in the bottom row shows $V_\phi/\sigma_R = 1$. The figure suggests that in the HESTIA simulations even the oldest stars of M31 and MW galaxies were formed as flattened 'discy' components which however formed and remained dispersion-dominated populations. The impact of the mergers is seen in some cases, which, however, do not have a dominant impact on the shape and overall kinematics of the in situ stars.

D'Souza & Bell 2018b; Bhattacharya et al. 2019). On the other hand, there are no significant mergers in the 09–18 MW analogue, where the most massive one has only about $\mu_* = 0.12$. Hence we observe a large diversity in the merger history of the simulated M31 and MW analogues in the HESTIA simulations. We remind the reader that the HESTIA simulations of the LG consist of M31 and MW analogues, aimed at reproducing a wide range of M31 and MW parameters and their satellites system, including the presence of LMC- and SMC-like object as well as an M33-like object (Libeskind et al. 2020).

For the M31 analogues, although the 09–18 simulation reproduces the last major merger very well, in other simulations the last major merger is less massive and happens earlier. This suggests that a significant variance in the satellite/merger mass function and accretion time can be the result of the dwarf merger history and coupling between small- and large-scale seeding in the initial conditions, which is not fully constrained in the HESTIA simulations.

A similar merger diversity is seen across the MW analogues, where the 09–18 and 17–11 simulations show slightly more massive mergers, compared to the estimated masses of different debris (e.g., GSE, Sequoia of $\approx(1-5) \times 10^8 M_\odot$, Kruijssen et al. 2020; Deason et al. 2019). This is likely linked to the longer accretion history of the simulated MW analogues which allows the dwarf galaxy progenitors to evolve longer prior to merging and being completely disrupted. The best candidate, in terms of the mergers parameters, in our sample is the MW analogue from the 37–11 simulation, where most of the mergers happened >8 Gyr ago and they have stellar masses $((1-5) \times 10^8 M_\odot)$, which is in reasonable agreement with the identified merger debris in the MW. Moreover, this galaxy has a few massive objects at $z=0$, similar to the Sgr, LMC and SMC systems. Therefore, the observed diversity of the merger histories in the HESTIA galaxies makes it possible to test both the evolution of the merger remnants in different environments and the response of the in situ stars to the merger events on different timescales.

# 4. Merger impact on the main progenitors

In this section we seek to explore various manifestations of the significant mergers on the in situ stellar populations.

## 4.1. Spatial axis ratio and $V_\phi/\sigma$ evolution of mono-age in situ populations

We begin by describing some general properties of the stellar populations of the HESTIA simulations. Various simulations suggest that stars that form earlier, that is to say out of turbulent gas at higher redshifts (see, e.g. Bournaud et al. 2005; Forbes et al. 2012), tend to be found with a thicker distribution (see, e.g. Minchev et al. 2013; Bird et al. 2013; Grand et al. 2016). Therefore, in addition to the upside-down formation of discs, one could expect some extra effects caused by the mergers. Indeed, Minchev et al. (2013) show that high-redshift stellar samples are both born hot and are additionally heated by mergers, in their hybrid chemo-dynamical model. Moreover, the authors found that merger-induced large-scale radial migration can redistribute hot-born stars (initially concentrated in the inner parts) and populate a more extended thick disc.

In Fig. 5 we show the evolution of the galaxy shape (the ratio between disc thickness and its radial extension) and the $V_\phi/\sigma_R$ relation for stars in different age bins (mono-age populations), where $V_\phi$ is the mean rotational velocity and $\sigma_R$ is the radial velocity dispersion of stars in a given age range. The shape of the galaxies is roughly estimated as the mean ratio between the mean vertical and radial positions of star particles of in situ populations. In Fig. 5 we see that older mono-age populations tend to be thicker and less supported by the regular rotation, however, even the oldest stars show a disc-like behaviour, being flattened towards the galactic plane. However, these older populations do not show substantial rotation compared to the random motions, because the $V_\phi/\sigma_R$ parameter is slightly below unity. With the stellar age decreasing, we see that stellar subsystems become
<figure><img src="image_7.png" /></figure>

**Fig. 6.** Impact of the most significant mergers and close pericentric passages of massive satellites on the star formation history of M 31 and MW HESTIA galaxies. For each galaxy, in the top panels we highlight the time of the close passages (<10 kpc from the host) and the stellar mass ($M_{\text{impact}}$ at the time of the close passage) of the dwarf galaxies which either merged (red) or survived at $z=0$ (blue). The bottom panels show the star formation rate as a function of the lookback time for three galactocentric regions: (0–1) $R_d$ (blue), (1–2) $R_d$ (red), and (2–3) $R_d$ (green), where $R_d$ is the disc scale length at $z=0$. The star formation histories at different galactocentric radii are scaled by the mean star formation rate whose value is mentioned in each panel with the same colour. The star formation histories of the HESTIA galaxies show a number of short timescale peaks in addition to the early rise and slow decrease of the star formation rate. In most cases, the peaks of the star formation correlate with either massive mergers or with the pericentric passages of massive dwarf galaxies. However, not all of the close encounters correlate with the star formation bursts inside the host galaxies.

thinner and more rotationally supported (see, e.g. Martig et al. 2014; Grand et al. 2016; Minchev et al. 2017).

Figure 5 suggests that the mergers moderately affect the axis ratio and kinematic characteristics of the mono-age populations. For most of the merger events we notice only temporal changes in the mono-age populations’ shape, which are seen as spikes in the top panels, corresponding to thickening and radial contraction of the mono-age populations. However, soon after the merger, the shape continues to evolve smoothly. The velocity-to-dispersion ratios also show short timescale variations during the mergers. Nevertheless, in the case of the most significant mergers, we detect a substantial change in the kinematics of the mono-age populations selected in the 2 Gyr range bins. In the M 31 (09–18) galaxy, the most massive merger at ≈3.5–5 Gyr (M5), all intermediate-age mono-age populations (3–7 Gyr old) heated up, which is seen in the decrease in the $V_\phi/\sigma$ ratio. The oldest populations are less sensitive to such an impact. A similar behaviour in the massive mergers’ impact is seen also in other galaxies, except for MW (09–18) and MW (37–11), which did not experience any massive mergers.

To summarise, although the HESTIA galaxies we analyse experienced a diverse assembly and merger history, the general trends in the kinematics and spatial distribution of the in situ mono-age populations are essentially the same. In particu-

lar, older mono-age populations are thicker and, despite some net rotation, are barely supported by the regular rotation. The overall impact of the mergers seems to be rather modest, but more details about the physical properties of the mono-age populations linked to the impact of the mergers are provided in the subsequent sections.

## 4.2. Merger impact on the star formation history, velocity dispersion and orbits of in situ stars

Pericentric passages of satellites may enhance the star formation either via contraction of the ISM or delivering more gas and thus fueling the star formation in the host (Mihos & Hernquist 1994; Hernquist & Mihos 1995). Although some observations of the MW rely on this idea (see, e.g. Ruiz-Lara et al. 2020), some simulations suggest that mergers and galaxy interactions do not always trigger starbursts in the hosts (Di Matteo et al. 2007, 2008; Martig & Bournaud 2008; Fensch et al. 2017; Di Cintio et al. 2021; Annem & Khoperskov 2022). To test this issue in the HESTIA galaxies, in Fig. 6 we show the star formation histories (SFHs) of the M 31 and MW analogues where we plot the age distribution weighted by the initial mass of star particles. Therefore, the mass loss due to stellar evolution does not affect the SFHs. We split the SFHs into three groups
for the stars that are located in different radial bins (0–1) $R_d$ (blue), (1–2) $R_d$ (red), and (2–3) $R_d$ (yellow) at $z = 0$, where $R_d$ is the stellar disc scale length (Libeskind et al. 2020). For each SFH, in the upper subpanels, the vertical lines show the time of close passages for dwarf galaxies (blue) and the mergers (red), where the height of the lines shows the total stellar mass of a perturber at that time ($M_{\text{impact}}$). In this figure, close passages include both pericentric passages and fly-bys of massive systems if they are closer than 10 kpc from the centre of the host galaxy.

One can see a number of peaks in the SFHs where some of them apparently correlate with the mergers or close encounters of dwarf galaxies. In most cases, we see a clear correlation between the mergers or close encounters and the peaks of the star formation in the hosts. However, the amplitude of the bursts does not seem to correlate with the mass of the perturbed, and even more extreme in many cases, passage of a massive system that does not always trigger star formation in the host. This behaviour is somewhat predicted by some other models depending on the interaction parameters and the amount of gas available for star formation (e.g., Di Matteo et al. 2007, 2008; Fensch et al. 2017). Another issue we can address is that most star formation bursts are seen for the stars residing in the inner parts of the galaxies (blue lines); whereas, the outer discs are likely less affected where, however, some bursts can be seen but have lower amplitudes. Overall, we suggest that if the peaks of the star formation are found in the host galaxies, they correlate with close passages of dwarf systems of the mergers. A more detailed analysis of the external impact on the star formation activity in the HESTIA galaxies, including the star formation inside the dwarf galaxies, is a subject of other work (see, e.g. Grand et al. 2020; Dillamore et al. 2022).

Next, we address the impact of the mergers on the orbital composition of in situ stellar populations. In Fig. 7 we show the mean $Z_{\text{max}}$ and mean $Z_{\text{max}}/R_{\text{max}}$ ratio as a function of stellar age. Overall, we see a gradual increase in both $Z_{\text{max}}$ and $Z_{\text{max}}/R_{\text{max}}$ with stellar age. At the same time, we can see that merger times coincide with a sharp increase in the mean $Z_{\text{max}}$ for all of the stars that existed in the host galaxy prior to the merger. Although the latest mergers seem to cause the largest increase in mean $Z_{\text{max}}$, in some galaxies one can still detect a few peaks of the vertical excursion of stars correlating with previous mergers. Interestingly, mergers also affect the $Z_{\text{max}}/R_{\text{max}}$ ratio, which is typically 0.1–0.2 for the youngest populations. However, for the stars that experienced the merger impact, this ratio can reach up to 0.5–0.6. This result suggests that the mergers are relatively more efficient at disc heating in the vertical direction, which inevitably causes disc flaring (Minchev et al. 2015).

Another issue is the impact of the mergers on the stellar velocity dispersion. As mentioned in the Introduction, it is expected that the GSE merger in the MW heated the pre-existing disc, resulting in the Plume/Splash feature (see Di Matteo et al. 2019; Belokurov et al. 2020). In Fig. 8 we show the age-radial velocity dispersion relation, ($\sigma_R$), for stars in different radial bins, as in the previous figures. In addition to the monotonically decreasing velocity dispersion with diminishing age, we can identify two different signals. Firstly, a number of peaks correlate with the times of external perturbations, suggesting that some mergers mostly affect stars that formed during or just before the interaction. Alternatively, these stars could have been formed in the tidal tails caused by the perturbations and thus have a locally higher velocity dispersion. The second type of feature is a step in the velocity dispersion, which results in a flattening for all stars that existed in the disc before the merger took place. This

appears to be the result of the most significant mergers. Such a step in the MW age-velocity relation has long been associated with the effect of the last massive merger (Quillen & Garnett 2001; Freeman & Bland-Hawthorn 2002; Minchev et al. 2018), now known as the GSE event.

# 5. The emergence of the in situ stellar halo

## 5.1. Kinematically defined stellar haloes

In this section, from the overall impact of the mergers on the discs, we move towards the analysis of the in situ stellar halo. Different definitions of the stellar halo have been adopted in the literature and here we use the one purely based on the stellar kinematics (Nissen & Schuster 2010). In particular, we analyse the structure of the Toomre diagram. In Fig. 9 we show the Toomre diagram for all stars (top), in situ only (middle), and all stars from the merger remnants (bottom, accreted only) at $z = 0$. To enhance the separation between the disc and the halo populations, we have omitted the centres of the HESTIA galaxies and consider the stars that are located in (0.5–3)$R_d$, which is accessible for the MW galactic surveys. In all the galaxies in our sample, the distributions on top show the presence of two distinct components, where one is centred near the Local Standard of Rest (LSR) and the others tend to show a slow prograde or even retrograde rotation with larger random velocities. The latter ones are mainly made of stars from disrupted dwarf galaxies, presented in the bottom panel which we explore in detail in Paper II. Therefore, the subsequent focus of our study is to explore the properties of the in situ stars, which, being weighted towards the LSR, show significant extension (see middle panels).

First, we look at the formation path of the in situ stellar halo, defined kinematically in the Toomre diagram. In particular, we investigate when the stars appear outside of the white circle in Fig. 9, which is often used as a threshold between the disc and halo stars in the MW. In Fig. 10, we present the origin of the stellar halo where the X axis corresponds to the time when certain in situ formed stars appear above the white line in Fig. 9, and the Y axis shows the age distribution of these stars. First, we can see that the kinematically defined in situ stellar halo continues to grow over time in all of the HESTA galaxies. This is likely due to constant perturbations from both orbiting dwarf galaxies and ancient mergers. This result should be taken into account in the MW analysis where the GSE-like merger is believed to kick stars out of the disc, and all of the halo stars (both in situ and accreted) that we observe gained halo-like motions at the time of the merger. In HESTIA galaxies, substantial halo stars arrive at the kinematically defined halo region some billion years ago after the mergers. This suggests that the overall mass of the in situ halo cannot be used to constrain the impact of the mergers directly. Nevertheless, in most cases, the massive mergers impact the in situ populations, making them move similarly to halo stars. This is seen as the enhancement of the density (fraction of the overall kinematically defined stellar halo) in Fig. 10 close to the mergers. The bottom panels of Fig. 10 show the cumulative fraction of the in situ stars, which have halo-like kinematics at the time of their formation. This fraction can be rather high (0.2–0.5) at the early phases of the galaxy assembly. However, it decreases rapidly over time, and at the present day, only about 10% of the in situ stars represent a genuine stellar halo. The rest is the result of the heating of the disc stars, thus suggesting that the heating of the pre-existing populations is the major source of the in situ stellar halo formation.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Impact of the most significant mergers and close pericentric passages of massive satellites on the orbital parameters of the in situ stars in M 31 and MW HESTIA galaxies. For each galaxy, in the top panels we highlight the time of the close passages (<10 kpc from the host) and the stellar mass ($M_{\text{impact}}$ at the time of the close passage) of the dwarf galaxies which either merged (red) or survived at $z = 0$ (blue). The two panels below show $Z_{\text{max}}$ and $Z_{\text{max}}/R_{\text{max}}$ as a function of the stellar age for three galactocentric regions: (0–1) $R_d$ (blue), (1–2) $R_d$ (red), and (2–3) $R_d$ (green), where $R_d$ is the disc scale length at $z = 0$.</figcaption></figure>

Another interesting feature we can notice in the Fig. 10 is that the youngest stars are the most sensitive populations to the impact of the mergers. At a given time, more of the contribution to the halo is from the stars that formed recently, which is seen as the overdensity along the one-to-one line. This is not expected for the gradual heating of the disc where stars cannot jump from the cold orbits to the halo-like ones without passing through the phase of moderate heating. An alternative opportunity is the external impact (tidal interaction with satellites and mergers), which affects populations on colder orbits more efficiently. How-

ever, the most massive mergers still affect even warm-orbit stars that formed several Gyr before but already experienced some heating. This is seen best in the case of the M5 merger at ≈3 Gyr in the M 31 (09–18) galaxy, but this is also noticeable in other models.

In order to test whether it is possible to date the mergers by using the in situ stellar halo formation history in Fig. 11 (blue lines), we show the age distribution of the kinematically defined in situ stellar halo (see Fig. 9) while the black lines show the accreted component. Interestingly, our definition of the in situ
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Impact of the most significant mergers and close pericentric passages of massive satellites on the age-radial velocity dispersion ($\sigma_R$) relation of M31 and MW HESTIA galaxies. For each galaxy, in the top panels we highlight the time of the close passages (<10 kpc from the host) and the stellar mass ($M_{\text{impact}}$ at the time of the close passage) of the dwarf galaxies which either merged (red) or survived at $z=0$ (blue). The bottom panels show the stellar velocity dispersion as a function of the stellar age time for three galactocentric regions: (0–1) $R_d$ (blue), (1–2) $R_d$ (red), and (2–3) $R_d$ (green), where $R_d$ is the disc scale length at $z=0$. The growth of the stellar velocity dispersion with the stellar age is remarkable in all the HESTIA galaxies, where in addition to the secular increase likely being caused by the scattering of stars on the spiral arms, we find some age ranges where the velocity dispersion increases faster. These features in stellar kinematics are linked to the impact of the mergers, which heat up the stars that formed before the merger. In the case in which several mergers occurred, the impact of the earliest ones is hard to see because the most recent ones, especially if they are more massive, heat up the disc to the maximal values.</figcaption></figure>

stellar halo results in a wide range of stellar ages where some
narrow peaks are clearly associated with the mergers (vertical
magenta lines). The peaks in the distribution of the stellar ages
are not very prominent because, once we assumed a constant
20% age uncertainty (red dashed lines), they were completely
erased, which is the result of a continuous inflow of stars to the
kinematically-defined region in the Toomre diagram. On aver-
age, the in situ stellar halo of the HESTIA galaxies is rather
young with the mean stellar age of 8–9 Gyr (a single exception
is M31 in the 09–18 simulation with a recent massive merger)
which is ≈1 Gyr younger compared to the accreted component.

Finally, we measured a spatial distribution of the in situ stel-
lar halo in such a way it could be compared to the MW data. In
Fig. 12 we show the relative fraction of the in situ stellar halo,
defined kinematically using the Toomre diagram in Fig. 9. In
the left panel, we present the in situ halo mass fraction without
spatial cuts, while on the right we show only the numbers avoid-
ing the disc region (|z| > 2 kpc). The distributions show that the
in situ (heated) populations represent ≈5–30% in the imminent

vicinity of the disc; this value decreases drastically at distances
larger than 10–15 from the galactic centre and stays almost con-
stant (2–5%) out to 50 kpc. The rest of the stellar halo is made
of merger debris. The right panel can be directly compared to
the results delivered by the H3 survey (Naidu et al. 2020), where
the MW in situ halo fraction decreases from ≈15% at 2–5 kpc
to ≈5% at 15–20 kpc. This suggests that, despite the relatively
more massive accreted systems, the HESTIA galaxies represent
the MW-disc stellar halo and disc behaviour well, and thus the
conclusions we make about the origin of the in situ stellar halo
can be further explored in the MW.

5.2. Searching for Splash/Plume features in HESTIA galaxies

Here, we link our finding to the known impact of the GSE merger
which is believed to be the major heating event in the MW.
First, in Fig. 13 we present the orbital eccentricity distribution
variation as a function of stellar age for all in situ stars (top)
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Toomre diagram for stars located in (1–3) $R_d$ at $z = 0$ ($R_d$ is the disc scale length from Libeskind et al. 2020). The top panels show all of the stars; middle panel corresponds to in situ populations only; and the bottom panel shows all accreted stars. In each panel, the distributions are normalised by the maximum value. White lines correspond to a threshold that is used to distinguish stars with disc-like kinematics from stars with halo-like kinematics (see, e.g. Nissen & Schuster 2010): $\sqrt{V_r^2 + V_z^2 + (V_\phi - V_{LSR})^2} < 180/240 \times V_{LSR}$, where 180 km s⁻¹ is the value typically used in the MW, 240 km s⁻¹ is the LSR in the MW (Reid et al. 2014) and $V_{LSR}$ is the mean rotational velocity in a flat part of the rotation curve in a given HESTIA galaxy at $z = 0$. Even taking some overlap into account, the distribution of accreted and in situ stars is very much so different, where in situ stars show mainly disc-like kinematics while the accreted stars show a diverse behaviour: from strongly counter-rotating (in MW 09–18) to non-rotating (M31 and MW in 17–11) or weakly co-rotating (M31 in 09–18), but also mixed kinematical properties (M31 and MW in 37–11).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Formation of the kinematically defined stellar halo. The top panels show the distribution of the ages of stars that start to contribute to the kinematically defined stellar halo (see Fig. 9) at a given time. Magenta vertical lines highlight five of the most significant mergers (M1–M5), in terms of the stellar mass ratio. Bottom panels show the cumulative fraction of stars that formed with halo-like kinematics. The figure suggests that at a given time, the newly formed stars contribute the most to the kinematically defined stellar halo region of the Toomre diagram. The reason for this type behaviour is that the youngest stars have colder kinematics compared to the pre-existing stars, thus, these recently that formed stars are the most sensitive to the external perturbations of the disc caused by the mergers. Since the fraction of stars formed kinematically hot does not exceed 10–12% of the in situ stars, the heating of the pre-existing populations is the dominant mechanism of the in situ halo formation in the HESTIA galaxies.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Age distribution of the kinematically defined stellar haloes (blue lines; see, Toomre diagram in Fig. 9 for the definition). The red dashed lines correspond to the distributions with the constant 20% age uncertainty. Black lines show the age distribution of the accreted component of the stellar halo. The numbers in the top left corner of each panel show the mean and the standard deviation of the stellar ages for the in situ (blue) and accreted (black) halo. The vertical magenta lines highlight five of the most significant mergers (M1–M5), in terms of the stellar mass ratio. The age distributions show a larger fraction of stars that formed inside the main progenitors close to the times of the mergers. However, the features completely vanished once the age errors were taken into account. Such a distribution transformation is not the result of the error itself but also because a substantial fraction of stars populating the kinematically defined halo between the merger events.</figcaption></figure>

and for in situ stars >1 kpc away from the disc plane (bottom). Similar to the previous figures, the most significant mergers are marked by the vertical lines where the length of the lines corresponds to the merger stellar mass ratio ($\mu_*$) where the panel size is assumed to be unity. The yellow lines show the mean eccentricity trends. Figure 13 clearly shows that the most massive mergers (the longest white vertical lines) result in a sharp increase in the mean eccentricity of the stars. In all of the HESTIA galaxies, the in situ formed halo stars show the breaks of the eccentricity distribution near the significant mergers. More generally, the most massive mergers lead to the substantial heating of the in-plane orbits of stars that can be found on the radial orbits with an eccentricity above 0.8. We note also, that the in situ halo stars have a broad distribution of eccentricities and even some of the oldest stars can be found on nearly circular orbits (see, e.g. Sestito et al. 2019).

Nissen & Schuster (2010) showed that the thick disc stars with halo-like kinematics found in solar vicinity data could be the result of the heating of the disc. More recently, using Gaia data Di Matteo et al. (2019) and Belokurov et al. (2020) quantified the impact of the GSE merger on the pre-existing MW stellar populations that are rotating more slowly, and even counter-rotating relative to disc stars. The excess of the slowly and counter-rotating (chemically defined) in situ stars is called Splash (Belokurov et al. 2020) or Plume (Di Matteo et al. 2019). In order to find similar features in the HESTIA galaxies in

Fig. 14 we show the variations of the azimuthal velocity as a function of stellar age (the dependence on the metallicity is discussed in Paper III) for all in situ stars (top) and in situ stars >1 kpc away from the disc plane (bottom). The picture here is somewhat similar to the eccentricity distribution in Fig. 13 where the merger events heat up the orbits of stars, thus decreasing their net orbital motion, and even pushing some of them on counter-rotating orbits. In Fig. 14, one can see that the HESTIA galaxies reveal a few episodes of disc heating, resulting in the formation of Splash/Plume-like stellar populations. For instance, in the 17–11 M31 galaxy a one-to-one M1 merger ≈10.5 Gyr ago results in the populations with $V_\phi$ ≈ 50 km s⁻¹ net rotation. The second massive merger at ≈7 Gyr ago (M5) creates the stars with $V_\phi$ ≈ 100 km s⁻¹ net rotation. It is important to note, however, that a very massive merger is not required for the formation of Splash/Plume-like populations. In particular, even a 9% stellar mass merger (M3) at ≈8 Gyr ago in the 09–18 MW galaxy is able to create a sharp decrease in the azimuthal velocity which is better seen for the halo stars (bottom row). Among the heated populations, we also notice the presence of stars with negative rotational velocity. Although the counter-rotating stars constitute between 7% and 24% of the stellar halo (see numbers in the bottom row in Fig. 14), these populations represent neither classical counter-rotating discs nor any distinct components of the galaxies while being the tails of the heated Splash/Plume-like stars.

Another interesting feature common in all of the HESTIA galaxies is that even strongly heated stellar populations (low $V_\phi$, or high eccentricity) still demonstrate a substantial net rotation with the direction aligned with most of the disc stars. This means that even the oldest stars inherit disc-like kinematics, suggesting a very fast collapse of the main progenitor at early times of the evolution and thus leaving no space for a non-rotating in situ stellar halo. Therefore, the stellar halo in HESTIA galaxies is predominantly made of disc stars heated by a few significant mergers.

### 5.3. $R_{max}-Z_{max}$ composition of the in situ stellar halo

In this section, we analyse the orbital composition of the stellar component of the HESTIA galaxies. For all star particles, we have reconstructed their orbital parameters by integrating their orbits in the fixed galactic potential at $z=0$ using the AGAMA code (Vasiliev 2019, see Sect. 2.2 for more details). In Fig. 15 we show the density distributions in the $R_{max}-Z_{max}$ plane, where $R_{max}$ is the apocentre and $Z_{max}$ is the maximum height from the mid-plane. Similar to a number of MW studies we present the distribution of stellar density in $R_{max}-Z_{max}$ coordinates (top) (see, e.g., Haywood et al. 2018; Koppelman et al. 2021). In the top panels, we detected a few diagonal overdensities (or ‘wedges’) separated by the prominent gaps. Koppelman et al. (2021) show that in the MW, these structures could be due to resonances and the result of chaotic diffusion while Haywood et al. (2018) suggest that the discrete wedges in the $R_{max}-Z_{max}$ plane – are all reminiscent of some impulsive heating of the early Galactic disc related to some accretion event(s). To clarify the composition of the $R_{max}-Z_{max}$ plane in the HESTIA galaxies, in Fig. 15 (bottom) we show the distribution of angles arctan $Z_{max}/R_{max}$ as a function of stellar age. In most of the HESTIA galaxies, the lowest-angle wedges appear after the latest significant merger. Some other wedges (horizontally aligned overdensities in the bottom panels) are made of stars that formed between some other, previous significant mergers. Therefore, we suggest that in the HESTIA galaxies, we detected a correlation between mergers and wedges in the $R_{max}-Z_{max}$ plane. However, some overlap between these structures and the
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Relative fraction of the in situ stellar halo as a function of the 3D galactocentric distance in different models at $z=0$. The left panel shows the data without any spatial selections, while in the right panel we adopted $|z| > 2$ kpc selection. The in situ halo is defined above the white line in the Toomre diagram in Fig. 9. In all the HESTIA galaxies, similar to the results of the H3 survey of the MW halo (Naidu et al. 2020), the in situ halo contribution rapidly decreases beyond ≈10 kpc and does not exceed 2–5% at larger distances. We note that recent massive mergers in M 31 (09–18) and M 31 (17–11) galaxies (see Fig. 4) result in the overdensity of in situ populations at ≈45 kpc.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Variations of the orbital eccentricity distributions as a function of the stellar age in M 31 and MW HESTIA galaxies. Distributions were normalised by the maximum value at a given age. The top row corresponds to all in situ stars while the bottom corresponds to in situ stars $>1$ kpc away from the galactic plane mostly excluding the disc component. White vertical lines highlight five of the most significant mergers (M1–M5), in terms of the stellar mass ratio. The mean eccentricity values are shown by the yellow (top) and black (bottom) lines. Similar to the rotational velocity behaviour (see Fig. 14), the distributions show several episodes of a sharp increase in the eccentricity which coincide with the mergers.</figcaption></figure>

uncertainties of stellar ages may complicate matching the wedges with accretion events in the MW.

## 6. Summary

We have analysed six M 31 and MW analogues from the HESTIA suite of cosmological hydrodynamics zoom-in simulations tailored to reproduce the properties of the LG galaxies. In this work, we have focussed our analysis on the impact of the mergers on the shape, kinematics and orbital composition of the in situ populations and the origin of the in situ component of the stellar halo. We found that all the M 31 and MW galaxies experi-

enced 10–40 mergers with dwarf galaxies in a stellar mass range of $10^6-2 \times 10^{10} M_{\odot}$; however, only on to four mergers (for different galaxies) are significant, where the stellar mass ratio is 0.2–1 (relative the host at the time of accretion) and all of the mergers (with a single exception) happened ≈7–12 Gyr ago. Our conclusions are the following.

- We found a certain correlation between close pericentric passages of massive satellites and mergers with the evolution of the star formation rate in the host galaxies. In particular, in most of the cases, significant peaks in the SFH correlate with the mergers; however, not all the mergers have the corresponding bursts of star formation in the host galaxy. This
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Structure of the azimuthal velocity as a function of the stellar age in M31 and MW HESTIA galaxies. Distributions were normalized by the maximum value at a given age. The top row shows all in situ stars, while the bottom one corresponds to the in situ stars >1 kpc away from the galactic plane, mostly excluding the disc component. Magenta vertical lines highlight five of the most significant mergers (M1–M5), in terms of the stellar mass ratio, where the length of the line shows the stellar mass ratio of the merger where the vertical panel’s size is assumed to be unity. The mean $V_\phi$ values are shown by the yellow lines. The numbers in white boxes show the mass fraction of counter-rotating stars (negative azimuthal velocity) among the sample presented in a given panel. The azimuthal velocity generally decreases with the stellar age due to secular disc heating; however, at the time of the mergers, it falls faster, resulting in slowly rotating stellar populations similar to the Splash/Plume stars discovered in the MW. We note, however, that in all of the galaxies, except for M31 and MW in the 37–11 simulation, the impact of several mergers is clearly seen in the stair-step-like distributions.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Orbital structure of in situ stellar halo. Top: stellar density distribution in $R_{max}-Z_{max}$ coordinates for all in situ stars (first row) and kinematically defined (see Toomre diagram in Fig. 9) in situ stellar halo (second row). Bottom: distribution of the angles between $R_{max}-Z_{max}$ as a function of the stellar age for all in situ stars (third row) and kinematically defined stellar halo (fourth row). The white vertical lines highlight five of the most significant mergers (M1–M5), in terms of the stellar mass ratio. Similar to a number of recent MW studies (see, e.g., Haywood et al. 2018; Koppelman et al. 2021) the HESTIA galaxies reveal a number of wedges in $R_{max}-Z_{max}$ coordinates. The bottom rows show how these wedges emerge over time. We find that the individual structures in $R_{max}-Z_{max}$ are made of stars predominantly made of stars that formed in between the significant merger events, suggesting a similar origin for similar structures in the MW.</figcaption></figure>
likely depends on the parameters of the interaction (orbit, relative mass) and the amount of gas available inside the host.

- The mergers induce heating of the in situ stars where the individual events are imprinted as the increase in the vertical motions of stars and the increase in the stellar velocity dispersion for the populations that existed in the disc prior to a given merger. We suggest that the latest significant merger makes the heating effect of the previous ones vanish, resulting in a flattening of the age-velocity dispersion relation for older stars.

- All the HESTIA galaxies clearly show the impact of the most massive mergers seen as a sharp increase in the eccentricity (and decrease in the $V_φ$) for all the stars that existed in the main progenitor before the merger thus nicely reproducing the Splash/Plume-like stellar populations recently discovered in the MW using Gaia data. Moreover, each HESTIA galaxy shows a few kinematically distinct (seen in the eccentricity and $V_φ$ distributions) populations caused by different mergers and about 7–14% of these stars have negative angular momentum.

- We have analysed the emergence of the in situ component of the stellar halo, defined kinematically using the Toomre diagram. We have shown that a substantial fraction (30–40%) of the inner stellar halo is formed in situ; while the in situ mass fraction does not exceed 5% of the total stellar halo mass in the outer parts (out to 20–50 kpc) of the HESTIA galaxies. The in situ stellar haloes defined kinematically have a broad range of ages with the mean values of 7–9 Gyr, which is between 1 and 2 Gyr younger than the accreted halo component. Although the age distributions of the in situ stellar halo show a number of prominent peaks correlating with the most significant merger events, its substantial mass emerges in between or well after the mergers.

- We have shown that the diagonal wedges in the $R_{max}-Z_{max}$ distributions can be linked to different merger events. In particular, the lowest angle (arctan $Z_{max}/R_{max}$) wedges correspond to the stars formed after the last significant merger, while the wedges with larger angles are made of stars formed in between the previous significant mergers.

In this work, we showed that mergers play a significant role in the evolution of the HESTIA galaxies, shaping the structure and properties of their in situ stellar populations. The collision and subsequent merging of smaller galaxies lead to the enhanced formation of new stars and the redistribution of pre-existing stellar populations that have been found in the MW thanks to the Gaia data and ongoing large-scale spectroscopic surveys.

**Acknowledgements.** We thank the anonymous referee for their valuable comments. S.K. acknowledges the HESTIA collaboration for providing access to the simulations. F.A.G. acknowledges support from ANID FONDECYT Regular 1211370 and from the ANID BASAL project FB210003. F.A.G. acknowledges funding from the Max Planck Society through a Partner Group grant. A.K. is supported by the Ministerio de Ciencia e Innovación (MICINN) under research grant PID2021-122603NB-C21 and further thanks Field Mice for Emma's house. J.S. acknowledges support from the French Agence Nationale de la Recherche for the LOCALIZATION project under grant agreements ANR-21-CE31-0019. Y.H. has been partially supported by the Israel Science Foundation grant ISF 1358/18. E.T. acknowledges support by ETAg grant PRG1006 and by EU through the ERDF CoE TK133. Software: IPython (Perez & Granger 2007), Astropy (Astropy Collaboration 2013, 2018), NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Pandas (McKinney 2010).

## References

Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 597, 21

Abadi, M. G., Navarro, J. F., & Steinmetz, M. 2006, MNRAS, 365, 747

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

Annem, B., & Khoperskov, S. 2022, MNRAS, submitted [arXiv:2210.17054]

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Atkinson, A. M., Abraham, R. G., & Ferguson, A. M. N. 2013, ApJ, 765, 28

Battaglia, G., Taibi, S., Thomas, G. F., & Fritz, T. K. 2022, A&A, 657, A54

Bell, E. F., Zucker, D. B., Belokurov, V., et al. 2008, ApJ, 680, 295

Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJ, 642, L137

Belokurov, V., Sanders, J. L., Fattahi, A., et al. 2020, MNRAS, 494, 3880

Bhattacharya, S., Arnaboldi, M., Caldwell, N., et al. 2019, A&A, 631, A56

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

Bird, J. C., Kazantzidis, S., Weinberg, D. H., et al. 2013, ApJ, 773, 43

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

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

Bournaud, F., Jog, C. J., & Combes, F. 2005, A&A, 437, 69

Brook, C. B., Kawata, D., Gibson, B. K., & Freeman, K. C. 2004, ApJ, 612, 894

Brook, C. B., Gibson, B. K., Martel, H., & Kawata, D. 2005, ApJ, 630, 298

Chiba, M., & Beers, T. C. 2000, AJ, 119, 2843

Christensen, C. R., Governato, F., Quinn, T., et al. 2014, MNRAS, 440, 2843

Ciucă, I., Kawata, D., Ting, Y. S., et al. 2023, MNRAS, in press https://doi.org/10.1093/mnrasl/slad033

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

Correa Magnus, L., & Vasiliev, E. 2022, MNRAS, 511, 2610

Cunningham, E. C., Sanderson, R. E., Johnston, K. V., et al. 2022, ApJ, 934, 172

Deason, A. J., Belokurov, V., & Evans, N. W. 2011, MNRAS, 416, 2903

Deason, A. J., Belokurov, V., Evans, N. W., & Johnston, K. V. 2013, ApJ, 763, 113

Deason, A. J., Mao, Y.-Y., & Wechsler, R. H. 2016, ApJ, 821, 5

Deason, A. J., Belokurov, V., & Sanders, J. L. 2019, MNRAS, 490, 3426

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451

Di Cintio, A., Mostoghiu, R., Knebe, A., & Navarro, J. F. 2021, MNRAS, 506, 531

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

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

Di Matteo, P., Bournaud, F., Martig, M., et al. 2008, A&A, 492, 31

Di Matteo, P., Lehnert, M. D., Qu, Y., & van Driel, W. 2011, A&A, 525, L3

Di Matteo, P., Haywood, M., Lehnert, M. D., et al. 2019, A&A, 632, A4

Di Matteo, P., Spite, M., Haywood, M., et al. 2020, A&A, 636, A115

Doumler, T., Courtois, H., Gottlöber, S., & Hoffman, Y. 2013a, MNRAS, 430, 902

Doumler, T., Gottlöber, S., Hoffman, Y., & Courtois, H. 2013b, MNRAS, 430, 912

Doumler, T., Hoffman, Y., Courtois, H., & Gottlöber, S. 2013c, MNRAS, 430, 888

D'Souza, R., & Bell, E. F. 2018a, MNRAS, 474, 5300

D'Souza, R., & Bell, E. F. 2018b, Nat. Astron., 2, 737

Duc, P.-A., Cuillandre, J.-C., Karabal, E., et al. 2015, MNRAS, 446, 120

Dupuy, A., Libeskind, N. I., Hoffman, Y., et al. 2022, MNRAS, 516, 4576

Fattahi, A., Deason, A. J., Frenk, C. S., et al. 2020, MNRAS, 497, 4459

Fensch, J., Renaud, F., Bournaud, F., et al. 2017, MNRAS, 465, 1934

Font, A. S., McCarthy, I. G., Crain, R. A., et al. 2011, MNRAS, 416, 2802

Forbes, J., Krumholz, M., & Burkert, A. 2012, ApJ, 754, 48

Freeman, K., & Bland-Hawthorn, J. 2002, ARA&A, 40, 487

Gaia Collaboration (Brown, A. G. A., et al.) 2016, A&A, 595, A2

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gallart, C., Bernard, E. J., Brook, C. B., et al. 2019, Nat. Astron., 3, 932

García de la Cruz, J., Martig, M., Minchev, I., & James, P. 2021, MNRAS, 501, 5105

Grand, R. J. J., Springel, V., Gómez, F. A., et al. 2016, MNRAS, 459, 199

Grand, R. J. J., Gómez, F. A., Marinacci, F., et al. 2017, MNRAS, 467, 179

Grand, R. J. J., van de Voort, F., Zjupa, J., et al. 2019, MNRAS, 490, 4786

Grand, R. J. J., Kawata, D., Belokurov, V., et al. 2020, MNRAS, 497, 1603

Hammer, F., Wang, J., Pawlowski, M. S., et al. 2021, ApJ, 922, 93

Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113

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

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85

Hernquist, L., & Mihos, J. C. 1995, ApJ, 448, 41

Hoffman, Y., & Ribak, E. 1991, ApJ, 380, L5

House, E. L., Brook, C. B., Gibson, B. K., et al. 2011, MNRAS, 415, 2652

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

Ibata, R., Irwin, M., Lewis, G., Ferguson, A. M. N., & Tanvir, N. 2001, Nature, 412, 49

Jean-Baptiste, I., Di Matteo, P., Haywood, M., et al. 2017, A&A, 604, A106

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2
Khoperskov, S., Minchev, I., Libeskind, N., et al. 2023a, A&A, 677, A90
(Paper II)

Khoperskov, S., Minchev, I., Libeskind, N., et al. 2023b, A&A, 677, A91
(Paper III)

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

Koppelman, H. H., Hagen, J. H. J., & Helmi, A. 2021, A&A, 647, A37

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

Libeskind, N. I., Carlesi, E., Grand, R. J. J., et al. 2020, MNRAS, 498, 2968

Lu, Y., Minchev, I., Buck, T., et al. 2022, Nat. Lett., submitted
[arXiv:2212.04515]

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

Mackereth, J. T., Schiavon, R. P., Pfeffer, J., et al. 2019, MNRAS, 482, 3426

Malin, D. F., & Carter, D. 1983, ApJ, 274, 534

Martig, M., & Bournaud, F. 2008, MNRAS, 385, L38

Martig, M., Minchev, I., & Flynn, C. 2014, MNRAS, 443, 2452

Martig, M., Minchev, I., Ness, M., Fouesneau, M., & Rix, H.-W. 2016, ApJ, 831, 139

Martínez-Delgado, D., Gabany, R. J., Crawford, K., et al. 2010, AJ, 140, 962

McCarthy, I. G., Font, A. S., Crain, R. A., et al. 2012, MNRAS, 420, 2245

McConnachie, A. W., Irwin, M. J., Ibata, R. A., et al. 2009, Nature, 461, 66

Mihos, J. C., & Hernquist, L. 1994, ApJ, 425, L13

Minchev, I., Chiappini, C., & Martig, M. 2013, A&A, 558, A9

Minchev, I., Martig, M., Streich, D., et al. 2015, ApJ, 804, L9

Minchev, I., Steinmetz, M., Chiappini, C., et al. 2017, ApJ, 834, 27

Minchev, I., Anders, F., Recio-Blanco, A., et al. 2018, MNRAS, 481, 1645

Minchev, I., Matijevic, G., Hogg, D. W., et al. 2019, MNRAS, 487, 3946

Monachesi, A., Gómez, F. A., Grand, R. J. J., et al. 2016, MNRAS, 459, L46

Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48

Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

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

Pagnini, G., Di Matteo, P., Khoperskov, S., et al. 2023, A&A, 673, A86

Pakmor, R., Springel, V., Bauer, A., et al. 2016, MNRAS, 455, 1134

Perez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21

Pillepich, A., Madau, P., & Mayer, L. 2015, ApJ, 799, 184

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Planck Collaboration XVI. 2014, A&A, 571, A16

Purcell, C. W., Bullock, J. S., & Kazantzidis, S. 2010, MNRAS, 404, 1711

Qu, Y., Di Matteo, P., Lehnert, M. D., & van Driel, W. 2011, A&A, 530, A10

Quillen, A. C., & Garnett, D. R. 2001, ASP Conf. Ser., 230, 87

Quinn, P. J., Hernquist, L., & Fullagar, D. P. 1993, ApJ, 403, 74

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130

Renaud, F., Agertz, O., Read, J. I., et al. 2021, MNRAS, 503, 5846

Ruiz-Lara, T., Gallart, C., Bernard, E. J., & Cassisi, S. 2020, Nat. Astron., 4, 965

Sales, L. V., Helmi, A., Abadi, M. G., et al. 2009, MNRAS, 400, L61

Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016, MNRAS, 457, 1931

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schuster, W. J., Moreno, E., Nissen, P. E., & Pichardo, B. 2012, A&A, 538, A21

Sestito, F., Longeard, N., Martin, N. F., et al. 2019, MNRAS, 484, 2166

Sestito, F., Buck, T., Starkenburg, E., et al. 2021, MNRAS, 500, 3750

Simon, J. D. 2018, ApJ, 863, 89

Sorce, J. G. 2015, MNRAS, 450, 2644

Sorce, J. G., & Tempel, E. 2018, MNRAS, 476, 4362

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289

Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629

Tissera, P. B., Scannapieco, C., Beers, T. C., & Carollo, D. 2013, MNRAS, 432, 3391

Tully, R. B., Courtois, H. M., Dolphin, A. E., et al. 2013, AJ, 146, 86

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22

Vasiliev, E. 2019, MNRAS, 482, 1525

Villalobos, Á., & Helmi, A. 2008, MNRAS, 391, 1806

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261

Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031

Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nat. Rev. Phys., 2, 42

McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, eds. S. van der Walt, & J. Millman, 56

Wetzel, A. R., Hopkins, P. F., Kim, J.-H., et al. 2016, ApJ, 827, L23

White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52

Wilson, M. L., Helmi, A., Morrison, H. L., et al. 2011, MNRAS, 413, 2235

Xiang, M., & Rix, H.-W. 2022, Nature, 603, 599

Zolotov, A., Willman, B., Brooks, A. M., et al. 2009, ApJ, 702, 1058