# The MAGPI survey: orbital distributions, intrinsic shapes, and mass profiles for MAGPI-like eagle galaxies using Schwarzschild dynamical models

Giulia Santucci, Claudia del P Lagos, Katherine E Harborne, Caro Derkenne, Adriano Poci, Sabine Thater, Richard M Mcdermid, J Trevor Mendel, Emily Wisnioski, Scott M Croom, et al.

► To cite this version:

Giulia Santucci, Claudia del P Lagos, Katherine E Harborne, Caro Derkenne, Adriano Poci, et al. The MAGPI survey: orbital distributions, intrinsic shapes, and mass profiles for MAGPI-like eagle galaxies using Schwarzschild dynamical models. Monthly Notices of the Royal Astronomical Society, 2024, 534 (1), pp.502 - 522. 10.1093/mnras/stae2114. insu-04705704

HAL Id: insu-04705704

https://insu.hal.science/insu-04705704v1

Submitted on 23 Sep 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The MAGPI survey: orbital distributions, intrinsic shapes, and mass profiles for MAGPI-like EAGLE galaxies using Schwarzschild dynamical models

Giulia Santucci ⋆, 1,2 Claudia Del P. Lagos ⋆, 1,2 Katherine E. Harborne ⋆, 1,2 Caro Derkenne ⋆, 2,3

Adriano Poci ⋆, 4 Sabine Thater, 5 Richard M. McDermid, 2,3 J. Trevor Mendel, 2,6 Emily Wisnioski ⋆, 2,6

Scott M. Croom ⋆, 2,7 Anna Ferré-Mateu ⋆, 8,9 Eric G. M. Muller ⋆, 2,6 Jesse van de Sande ⋆, 2,10

Gauri Sharma ⋆, 11,12,13,14 Sarah M. Sweet, 2,15 Takafumi Tsukui ⋆, 2,6 Lucas M. Valenzuela ⋆, 16

Glenn van de Ven ⋆ 5 and Tayyaba Zafar ⋆ 3

*Affiliations are listed at the end of the paper*

Accepted 2024 September 6. Received 2024 August 16; in original form 2024 June 3

## ABSTRACT

Schwarzschild dynamical models are now regularly employed in large surveys of galaxies in the local and distant Universe to derive information on galaxies’ intrinsic properties such as their orbital structure and their (dark matter and stellar) mass distribution. Comparing the internal orbital structures and mass distributions of galaxies in the distant Universe with simulations is key to understanding what physical processes are responsible for shaping galaxy properties. However, it is first crucial to understand whether observationally derived properties are directly comparable with intrinsic ones in simulations. To assess this, we build Schwarzschild dynamical models for MUSE-like IFS (integral field spectroscopy) cubes (constructed to be like those obtained by the Middle Ages Galaxy Properties with Integral Field Spectroscopy, MAGPI survey) of 75 galaxies at $z \sim 0.3$ from the EAGLE simulations. We compare the true particle-derived properties with the galaxies’ model-derived properties. In general, we find that the models can recover the true galaxy properties qualitatively well, with the exception of the enclosed dark matter, where we find a median offset of 48 per cent, which is due to the assumed Navarro–Frenk–White profile not being able to reproduce the dark matter distribution in the inner region of the galaxies. We then compare our model-derived properties with Schwarzschild models-derived properties of observed MAGPI galaxies and find good agreement between MAGPI and EAGLE: the majority of our galaxies (57 per cent) have non-oblate shapes within 1 effective radius. More triaxial galaxies show higher fractions of hot orbits in their inner regions and tend to be more radially anisotropic.

**Key words:** galaxies: fundamental parameters – galaxies: kinematics and dynamics – galaxies: stellar content.

## 1 INTRODUCTION

Over the last two decades, the rapidly growing field of integral field spectroscopy (IFS) has significantly deepened our knowledge of stellar structures in galaxies in the local Universe, across a wide range of masses, morphologies, and environments. IFS instruments are the only ones that allow for the stellar and gas phase properties of galaxies to be spatially and simultaneously mapped. IFS surveys, such as ATLAS3D (Cappellari et al. 2011), the Sydney-AAO Multi-Object Integral-Field Spectrograph (SAMI) Galaxy Survey (Croom et al. 2012; Bryant et al. 2015), the Calar Alto Legacy Integral Field Area Survey (Sánchez et al. 2012), MASSIVE (Ma et al. 2014), and the Mapping Nearby Galaxies at Apache Point Observatory survey (Bundy et al. 2015), have contributed to significantly expand our understanding of galaxy kinematics and their connection to intrinsic galaxy properties and their environment (see review by Cappellari

2016). For example, it has been possible to unveil various correlations between the proxy for the stellar spin parameter, $\lambda_{Re}$, which provides a measurement of how rotationally supported a galaxy is, and different galaxy properties. $\lambda_r$ has previously been used to separate slow-rotating galaxies (i.e. galaxies whose kinematics are dominated by random motions) from fast-rotating galaxies (Emsellem et al. 2007, 2011; Cappellari 2016) and has been observed to be strongly correlated with stellar mass, so that the fraction of galaxies with low $\lambda_{Re}$ (slow-rotating systems) increases with increasing stellar mass (Emsellem et al. 2011; Falcón-Barroso et al. 2011; van de Sande et al. 2017; Veale et al. 2017; Brough et al. 2017; Wang et al. 2020), as well as correlated with environment (van de Sande et al. 2021), and age (van de Sande et al. 2018; Croom et al. 2024).

Thanks to these IFS surveys, kinematic maps as well as stellar age and metallicity maps for thousands of nearby galaxies have been made available. By applying orbit-based dynamical models to the IFS data (e.g. Cappellari et al. 2007; van den Bosch et al. 2008; Zhu et al. 2018c; Jin et al. 2020; Santucci et al. 2022, 2023; Thater

<footnote>*E-mail: giulia.santucci@uwa.edu.au</footnote>
et al. 2023a; Derkenne et al. 2024, hereafter D24), we can separate dynamically derived properties of galaxies, such as intrinsic shape, orbital components, velocity anisotropy, and inner mass distribution. These dynamical models can therefore provide insight into the 3D internal kinematic structures underlying the line-of-sight kinematics. Several different implementations of the Schwarzschild method, with varying degrees of symmetry, have been described (e.g. Cretton et al. 1999; Gebhardt et al. 2003; Valluri, Merritt & Emsellem 2004; van den Bosch et al. 2008; Vasiliev & Athanassoula 2015; Vasiliev & Valluri 2020; Neureiter et al. 2021). The Schwarzschild method has been used to study the internal stellar structure of globular clusters (van de Ven et al. 2006), nuclear star clusters (Feldmeier-Krause et al. 2017; Fahrion et al. 2019; Thater et al. 2023b), and galaxies (e.g. Thomas et al. 2007; Cappellari et al. 2007; van de Ven, de Zeeuw & van den Bosch 2008; Feldmeier-Krause et al. 2017; Poci et al. 2019; Jin et al. 2020; Santucci et al. 2022, 2023; Pilawa et al. 2022; Thater et al. 2023a; Mehrgan et al. 2024; Jin et al. 2024), to measure the mass of supermassive black holes (van der Marel et al. 1998; Verolme et al. 2002; Gebhardt et al. 2003; Valluri et al. 2004; Krajnović et al. 2009; Rusli et al. 2013; Seth et al. 2014; Thater et al. 2017, 2019, 2022a; Liepold et al. 2020; Quenneville, Liepold & Ma 2021, 2022; Neureiter et al. 2023; de Nicola et al. 2024), and has also been used to identify accreted galactic components (e.g. Zhu et al. 2020, 2022; Poci et al. 2021).

However, to obtain a complete picture of galaxy evolution and understand how galaxies evolve into slow- and fast-rotating systems we also need observations at higher redshifts. The Middle Ages Galaxy Properties with Integral Field Spectroscopy (MAGPI) survey (Foster et al. 2021) provides observations for a substantial sample of galaxies with both stellar and gas-phase observations at redshifts $z \sim 0.3$, an epoch previously only sparsely studied. The absence of substantial IFS observations in this epoch has limited our understanding of galaxy evolution during the Universe's middle ages when morphology, angular momentum, and star formation activity evolve rapidly, with the environment potentially playing a key role (e.g. Peng et al. 2010; Choi et al. 2018). The higher redshift MAGPI sample, with its wide range in environment, can allow us to better constrain the internal properties of massive galaxies, and shed light on their later accretion history. In particular, thanks to the MAGPI theoretical library (Harborne et al., in preparation), a like-to-like comparison with simulations is now possible. D24 built Schwarzschild dynamical models for 22 galaxies in the MAGPI survey. They found that the inner mass distribution of galaxies tends to be non-homogeneous, varying both between galaxies and within galaxies. They also found that the majority of the galaxies in their sample (over 85 per cent) had non-oblate shapes, and that the shape of a galaxy strongly correlates with the fraction of hot orbits within one effective radius.

MAGPI theoretical library (Harborne et al., in preparation) provides the mock data from galaxy evolution simulations cubes with the same observational conditions and post-process of the MAGPI data. This enables us to assess the robustness of our modelling methods using simulated galaxies with ground truth. Such comparison will provide insights into how stellar orbital properties and mass distribution have evolved as simulations provide a mean to trace the evolutionary paths of galaxies over cosmic time.

In this paper, we build orbit-based triaxial dynamical model fits of 75 MAGPI-like mock-observations of galaxies from the EAGLE simulations. We first compare the model-derived properties with the true properties, and then with the results from D24. In Section 2, we describe the sample of galaxies and the data available for this analysis; we describe how we construct the Schwarzschild dynamical

<table><caption>**Table 1. Specifications of the EAGLE Ref-L100N1504 simulation used in this paper.**</caption><thead><tr><th></th><th>Property</th><th>Units</th><th>Value</th></tr></thead><tbody><tr><td>(1)</td><td>Gas particle mass</td><td>[M<sub>☉</sub>]</td><td>1.81 × 10<sup>6</sup></td></tr><tr><td>(2)</td><td>Dark matter (DM) particle mass</td><td>[M<sub>☉</sub>]</td><td>9.7 × 10<sup>6</sup></td></tr><tr><td>(3)</td><td>Average stellar particle mass at birth</td><td>[M<sub>☉</sub>]</td><td>2.26 × 10<sup>6</sup></td></tr><tr><td>(4)</td><td>Average stellar particle mass at z ~ 0.27</td><td>[M<sub>☉</sub>]</td><td>1.3 × 10<sup>6</sup></td></tr><tr><td>(5)</td><td>Softening length</td><td>[ckpc]</td><td>2.66</td></tr><tr><td>(6)</td><td>Max. gravitational softening</td><td>[pkpc]</td><td>0.7</td></tr></tbody></table>

*Notes.* The rows list: (1) initial particle masses of gas and (2) DM, (3) the average stellar particle mass at birth and (4) at the redshift of interest $z \sim 0.27$, (5) comoving Plummer-equivalent gravitational softening length, (6) maximum physical gravitational softening length. Units are indicated in each row. EAGLE adopts (5) as the softening length at $z \ge 2.8$, and (6) at $z < 2.8$. This simulation has a side length of $L = 100$ cMpc. Here, pkpc and ckpc refer to proper and comoving kpc, respectively, and cMpc to comoving Mpc.

models in Section 3; we present our results in Section 4, and we discuss them in Section 5. Our conclusions are given in Section 6.

Throughout the paper, we adopt a lambda-cold dark matter cosmology with $\Omega_m = 0.3175$ and $H_0 = 67.1$ km s$^{-1}$ Mpc$^{-1}$, corresponding to the Planck 2013 cosmology (Planck Collaboration XVI 2014).

## 2 DATA

In this work, we use the Schwarzschild modelling tool DYNAMITE (Jethwa et al. 2020; Thater et al. 2022b). Our aim in this paper is first to test the robustness of its derived properties. Secondly, we aim to compare the properties derived for simulated galaxies with those obtained for MAGPI galaxies. For a fair comparison, we use the EAGLE – MAGPI theoretical library from the MAGPI Theory Data Release (Harborne et al., in preparation). This library contains kinematic maps of MAGPI-like galaxies that can be analysed with observational methods. We introduce the EAGLE simulations in Section 2.1, the MAGPI theoretical library in Section 2.2, the properties derived from the particle data in Sections 2.3 and 2.4, and the sample selection in Section 2.5.

### 2.1 The EAGLE simulation

The EAGLE simulations suite (Schaye et al. 2015, hereafter S15; Crain et al. 2015) employs a modified version of the N-body Tree-Particle-Mesh smoothed particle hydrodynamics (SPH) code GADGET-3 (Springel 2005) with subgrid recipes for radiative cooling, star formation, stellar evolution, chemical enrichment, stellar, and black hole feedback. The suite consists of a large number of cosmological hydrodynamic simulations with different resolutions, volumes, and subgrid physics models, adopting the Planck Collaboration XVI (2014) cosmology. S15 introduced a reference model, within which the parameters of the subgrid models governing energy feedback from stars and accreting black holes were calibrated to ensure a good match to the $z = 0.1$ galaxy stellar mass function, the size-mass relation of present-day star-forming galaxies and the black hole-stellar mass relation (see Crain et al. 2015, for details). In this work, we take the largest volume with the highest resolution (Ref-L100N1504) presented in S15. Table 1 summarizes the parameters of this simulation.

A major aspect of the EAGLE project is the use of state-of-the-art subgrid models that capture unresolved physics (i.e. processes happening within the resolution limit). The subgrid physics modules adopted by EAGLE are: (i) radiative cooling and photoheating
(Wiersma, Schaye & Smith 2009a), (ii) star formation (Schaye & Dalla Vecchia 2008), (iii) stellar evolution and chemical enrichment (Wiersma et al. 2009b), (iv) stellar feedback (Dalla Vecchia & Schaye 2012), and (v) black hole growth and active galactic nucleus feedback (Rosas-Guevara et al. 2015).

The EAGLE simulations were performed using an extensively modified version of GADGET-3 (Springel 2005). Among those modifications are updates to the SPH technique, which are collectively referred to as ‘Anarchy’ (see Schaller et al. 2015 for an analysis of the impact of these changes on the properties of simulated galaxies compared to standard SPH). We use SUBFIND (Springel et al. 2001; Dolag et al. 2009) to identify self-bound overdensities of particles within haloes (i.e. substructures). These substructures are the galaxies in EAGLE.

Since we are interested in evaluating the robustness of the derived properties from Schwarzschild dynamical models for galaxies in the MAGPI survey, we use snapshot 25 at $z = 0.27$, which is closest to the nominal 0.3 target redshift of the MAGPI survey. Due to the high computational time required for the models to be computed (around 400 CPU hours per galaxy, on average), we select a subsample of EAGLE galaxies in order to match the MAGPI Primary Targets (described in Section 2.2). For each of these galaxies, we follow the steps outlined below:

(i) we calculate the true values for each of the galaxy properties we are interested in from the simulation stellar particle data (Sections 2.3 and 2.4),

(ii) for each mock observation, we emulate the process of fitting a dynamical model as for real data, by:

(a) build Schwarzschild dynamical models using DYNAMITE,

(b) select the model that best matches our data, using a $\chi^2$ comparison,

(c) from the best-fitting model, we derive the galaxy properties we are interested in (such as enclosed matter, intrinsic shape, and orbital distribution),

(iii) we compare the model-derived properties with the properties inferred from the particle data,

(iv) we compare the model-derived properties with the properties derived for MAGPI galaxies by D24.

## 2.2 The MAGPI theoretical library

A suite of cosmological simulations is included in the MAGPI data set. The MAGPI Survey uses existing cosmological hydrodynamical simulations and retrieves data from EAGLE (Schaye et al. 2015; Crain et al. 2015), MAGNETICUM (Teklu et al. 2015; Schulze et al. 2018), HORIZON-AGN (Dubois et al. 2016), ILLUSTRIS-TNG100 (Pillepich et al. 2018; Naiman et al. 2018; Springel et al. 2018; Nelson et al. 2019), and the chemodynamical simulation of Taylor et al. (2015) and Taylor & Kobayashi (2017). For this work, the first in a series, we make use of the MAGPI theoretical data set from EAGLE. In particular, we use the ready-to-use mock kinematic data cubes of the 481 EAGLE galaxies classified as MAGPI Primary Targets (Foster et al. 2021), created using SIMSPIN v2.8.3.<fnref n="1" /> Here, we provide a summary of how these mock observations are created. The full description can be found in Harborne et al. (2020, 2023).

SIMSPIN acts as a virtual telescope wrapper, where the parameters of the virtual telescope can be chosen in a variety of different ways

in order to reproduce observations of any IFS instrument. For the MAGPI theoretical data set, the mock datacubes are created to match the MAGPI observations: a spatial pixel scale of 0.2 arcsec, a velocity pixel scale of 1.25 Å, an average line-spread function of full-width half-maximum = 2.63 Å, and observational ‘noise’. The galaxies are projected to a redshift distance of $z = 0.3$ to match the observation specifications. To best reproduce the variety of inclination angles we see in observations, the mocks are generated at a random inclination angle. The resulting data cubes provide values for the observed flux, mass, $v$, $\sigma$, as well as the higher order kinematics moments $h_3$ and $h_4$ (Harborne et al. 2023).

The MAGPI EAGLE theoretical library also provides different measurements for the effective radius, $R_e$. Here, we take the photometric major axis as measured using the PROFOUNDGETELLIPSES function (Robotham et al. 2018), as our $R_e$.

SIMSPIN also gives the user the 3D shape of the galaxy at $R_e$. The shapes are computed using the method described in the work of Bassett & Foster (2019), which in turn is based on the work from Li et al. (2018) and Allgood et al. (2006). A short summary of the method is presented here, and the full description can be found in Harborne et al. (2023). See Zemp et al. (2011) and Valenzuela et al. (2024) for an overview of different methods.

To start, the initial distribution of stellar particles is assumed to be an ellipsoid with axis ratios $p = q = 1$ (i.e. a sphere, where $p = b/a$ and $q = c/a$, with $a, b$, and $c$ representing the axis lengths in decreasing size such that $a \ge b \ge c$ and $p \ge q$, by necessity). This ellipsoid is grown from the centre position assigned until it contains half the total stellar mass within the subhalo. Stellar particles are then used within the region to measure the reduced inertia tensor. The eigenvalues and eigenvectors of this tensor, $I_{ij}$, give the orientation and distribution of matter within the ellipsoid. The ellipsoid is then deformed to match the distribution of stellar particles, and the procedure is repeated, with a new ellipsoid grown to contain half the mass that reflects the new matter distribution. This process is repeated until the axis ratios $p$ and $q$ stabilize over 10 iterations.

We take these final $p$ and $q$ as our axis ratios at $1 R_e$ and we compute the triaxial parameter $T_{Re}$, calculated at $1 R_e$ and defined in Binney & Tremaine (2008):

$$T_{Re} = (1 - p_{Re}^2)/(1 - q_{Re}^2). \qquad (1)$$

## 2.3 Particle data – mass profiles

The masses of EAGLE galaxies are measured using spherical apertures. This gives similar results to the 2D Petrosian aperture used in observational works, but provides an orientation-independent mass measurement for each galaxy. Following previous studies based on the EAGLE simulations, we define the total stellar mass of a galaxy as the mass of the stellar particles within a 50 pkpc (pkpc refers to proper kpc) aperture centred at the centre of potential (e.g. Schaye et al. 2015). We note that selecting a different aperture (i.e. 30 pkpc) does not change our results. Similarly, dark matter (DM) masses and total masses are computed in spherical apertures (aperture values are publicly available and can be accessed from the EAGLE database<fnref n="2" />). For each of our galaxies, we calculate the enclosed stellar, dark, and total mass within $1 R_e$ by interpolating the values between the two radial apertures closest to $1 R_e$.

We note that due to the resolution of the simulation, dynamical properties of galaxies are not reliable for stellar masses $\lesssim 10^{9.5} M_{\odot}$

<footnote n="1">¹https://kateharborne.github.io/SimSpin/</footnote>

<footnote n="2">²https://virgodb.dur.ac.uk/</footnote>
(see appendix A in Lagos et al. 2017; Ludlow et al. 2023). Moreover, Ludlow et al. (2023) showed that a stellar mass $\gtrsim 10^{10} M_{\odot}$ is needed for stellar kinematics to be resolved at the 3D half-mass radii of EAGLE galaxies; at lower masses, gravitational scattering between stellar and DM particles affects their kinematics and spatial distribution. It is important to note that DM mass resolution and gravitational scattering can also affect the circularity distributions of stellar disc particles by transforming galactic discs from flattened structures into rounder spheroidal systems, causing them to lose rotational support in the process (Wilkinson et al. 2023).

## 2.4 Particle data – orbital distributions

We use the orbital distributions and decomposition of stellar particles in EAGLE galaxies calculated in Santucci et al. (2024). Here, we provide a summary of the method.

To get to the orbital distribution and decomposition of stellar particles in EAGLE galaxies, we first define the z-axis of a galaxy as that parallel to the stellar angular momentum vector of the galaxy computed with all stellar particles of the subhalo with the spatial reference frame being centred on the position of the most bound particle, and the velocity reference being the mean mass-weighted velocity of all stellar particles in the subhalo.

We then calculate the ratio between the z-component of the specific angular momentum of a stellar particle relative to the specific angular momentum of a circular orbit of the same energy, $\lambda_z$ (a.k.a. circularity parameter), following the method introduced by Abadi et al. (2003). With the calculation above, we construct a 2D distribution of $\lambda_z$ versus projected distance (in the xy plane of the galaxy) for each galaxy.

We then separate orbits into four different components as follows:

(i) cold orbits selected as $\lambda_z > 0.8$, are those with near-circular orbits;

(ii) warm orbits selected as $0.25 < \lambda_z < 0.8$, are somewhat disc-like, but with a considerable amount of angular momentum out of the disc (z) plane;

(iii) hot orbits selected as $-0.25 < \lambda_z < 0.25$, are those with near radial orbits;

(iv) counter-rotating orbits are those with $\lambda_z < -0.25$.

These thresholds were adopted to follow the analysis of Zhu et al. (2018a), Jin et al. (2020), Santucci et al. (2022), and D24, and allow for easier comparison with observations.

However, these values of the circularity parameter, $\lambda_z$, are instantaneous, while the model-derived circularities are time-averaged. This can create biases when comparing the circularity distributions from the models and the true ones. The stellar particles in our simulated galaxies do not necessarily conserve $\lambda_z$ when orbiting in the potential, especially those particles on radial/box orbits. These orbits (as mentioned above) are usually characterized by $\lambda_z \sim 0$. This means that in the instantaneous orbital distributions particles on radial/box orbits can appear to have the same orbital circularity as particles on cold/warm orbits (i.e. ‘spreading’ the circularity distribution). A way to obtain the orbits’ circularity in a similar manner to the time-averaged obtained from the dynamical modelling is to freeze the potential, integrate the particle orbits in the potential, and calculate the average values along the orbits. This method, however, is computationally and time expensive. Here for simplicity, since we aim to compare the orbital distributions of simulated galaxies with those calculated in observations and mock observations using Schwarzschild dynamical models, which are time-averaged, we calculate $\lambda_z$ for each particle in three different snapshots (snapshot

$24$ – at $z = 0.36$, snapshot $25$ – at $z = 0.27$, and snapshot $26$ – at $z = 0.18$) and we then average them to get an orbital distribution that is as close as possible to a time-averaged one. A similar approach has also been used by Reiter et al. (in preparation), for TNG50 galaxies, who found this method to be able to better reproduce time-averaged orbits. We note, however, that TNG50 has a high time-resolution between snapshots ($\sim 150$ Myr), while the time range between two EAGLE snapshots is close to 1 Gyr.

From the orbital distributions we calculate the fraction of orbits in each component, within $1 R_e$, by summing the mass of each stellar particle within $1 R_e$ that is in the corresponding orbital family, and dividing by the total stellar mass within $1 R_e$.

## 2.5 EAGLE sample selection

We select our EAGLE galaxies to match the MAGPI primary targets. This initial selection gives us 481 galaxies with $10.0 < \log M_{\star,50kpc}/M_{\odot} < 11.54$. We exclude four of these galaxies where $R_e < 0.6$ arcsec, due to their spatial size being smaller than the mock instrumental spatial resolution. Since our focus is on the dynamical properties within $1 R_e$, we select galaxies with robust kinematic data (with at least 50 particles in each spaxel; Jiménez et al. 2023, appendix and Harborne et al. in preparation) available up to at least $1 R_e$. We carefully visually inspect the 350 selected galaxies and we exclude all galaxies whose kinematics are influenced by mergers, or that have a bright secondary object within $1 R_e$ in their stellar velocity field. This leaves us with 259 galaxies with $10.12 < \log M_{\star,50kpc}/M_{\odot} < 11.54$. After further data quality requirements imposed in the below sections, we have a final sample of 75 galaxies.

# 3 SCHWARZSCHILD MODELS – DYNAMITE

The MAGPI theoretical data set was created with the intention to be used as a like-to-like comparison for MAGPI observational data. As such, we strive to apply methods as closely as possible to those used for the MAGPI data in order to reduce any possible bias due to different methods.

The two basic requirements when constructing a Schwarzschild model are 2D stellar kinematics with as many higher order moments as feasible (dependent on the quality of the spectral data), and a description of the galaxy’s surface brightness. In the following sections, we describe our modelling of the galaxy surface brightness and our extraction of the 2D stellar kinematics. Additionally, since models for the MAGPI sample were obtained by transforming the surface brightness models to stellar mass models via spectral stellar population measurements to infer stellar mass-to-light ratios, we also choose to use stellar mass maps instead of luminosity maps to construct our mass distribution.

## 3.1 Data preparation

### 3.1.1 Multi-Gaussian expansion profiles

We parametrize the galaxy mass using multi-Gaussian expansion (MGE) method, which consists of a series of concentric 2D Gaussians (Emsellem, Monnet & Bacon 1994), using the PYTHON package MGEFIT (Cappellari 2002). This method enables us to take the point spread function into account; given a value of the inclination and assuming an intrinsic shape, the MGE model can be deprojected analytically, and used to describe the intrinsic mass density of the stellar tracer within a given gravitational potential.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The stellar mass model, for example, galaxy EAGLE 9498015. The top panel shows the isophotes of the galaxy in black and of the best MGE model overlaid in red. The bottom panel shows the residuals of the model as (data model) divided by data image.</figcaption></figure>

We fit the MGEs by following the process described in Derkenne et al. (2023) and D24 as closely as possible. Similarly to them, we require all components to have the same position angle, as a change in position angle was not warranted by the fit quality of the models compared to the data. The fits are applied to the mass images created by SIMSPIN. We convert the dispersion of each Gaussian component from units of pixels to arcseconds by multiplying by the mock-MUSE pixel scale (0.2 arcsec). We show an example of an MGE fit in Fig. 1.

### 3.1.2 Stellar kinematics

As mentioned in Section 2.2, stellar kinematic maps are one of the outputs from SIMSPIN provided for the MAGPI theoretical data set. These include line of sight velocity ($v$), velocity dispersion ($\sigma$), and $h_3$ and $h_4$, the Gauss-Hermite moments which describe the skewness and kurtosis of the profile, respectively (Gerhard 1993; van der Marel 1993).

We apply Voronoi binning to the kinematic maps – using the number of particles ($N_p$) as a proxy for signal-to-noise ratio. We select a target $N_p = 200$ and apply the Voronoi binning with the PYTHON package VORBIN (Cappellari & Copin 2003). We find that this threshold is the best compromise to achieve robust values for the

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Number of Voronoi bins within $1 R_e$ for galaxies that meet our quality criteria as a function of stellar mass within 50 kpc for our EAGLE sample (Section 2.5). We exclude from our sample galaxies that have fewer than 80 Voronoi bins within $1 R_e$ (horizontal grey dashed line). Galaxies that do not meet the required number of Voronoi bins within $1 R_e$ are shown as open circles. The histograms on the top and on the right panels show the distribution in mass and Number of Voronoi bins, respectively, for the galaxies in our sample before (hatched histograms) and after (coloured histograms) the cut in Voronoi bins. The average number of Voronoi bins for each galaxy in the final sample is 99.</figcaption></figure>

kinematic maps (ensuring the outer regions, where the uncertainties are higher and the kinematic less reliable, have less weight) while keeping a high enough number of bins for our galaxies to be properly spatially sampled (as demonstrated by Harborne et al., in preparation). For each Voronoi bin, $v$, $\sigma$, $h_3$, and $h_4$ are taken as the weighted average (by number of particles) of the values in each spaxel belonging to the Voronoi bin. Spaxels with more than 200 particles are left unbinned. We show the number of bins within $R_e$ as a function of stellar mass in Fig. 2. To ensure we have enough resolution elements in our galaxies to construct robust models, we further exclude all the galaxies with less than 80 Voronoi bins within $1 R_e$.<fnref n="3" /> This cut leaves us with 75 galaxies. It is important to note that this final selection excludes most of the smallest galaxies from our sample, as well as low-mass edge-on (with angles as defined in Section 2.2) galaxies, since they have fewer bins within $1 R_e$ than more face-on galaxies of similar size and mass.

## 3.2 Construction of Schwarzschild models

Schwarzschild (1979) models are a flexible but computationally expensive method that allows us to model triaxial stellar systems in three steps: first, we construct a model for the underlying gravitational potential of the galaxy; secondly we calculate a representative library of orbits in these gravitational potentials; and then we find a combination of orbits that can reproduce the observed kinematic maps and luminosity distribution. These steps are fully described in van den Bosch et al. (2008) and Zhu et al. (2018b). Here, we

<footnote n="3">³We find that the properties derived from models built using maps with fewer than 80 bins are not well constrained, even when thousands of models are calculated.</footnote>
use the triaxial DYNAMITE code<fnref n="4" /> (Jethwa et al. 2020) to construct the Schwarzschild models, with orbit mirroring corrections included (Quenneville, Liepold & Ma 2022; Thater et al. 2022b). The models are computationally expensive as a large library of orbits in a given gravitational potential are integrated for a set number of orbital periods (to ensure reliable convergence), and the resulting surface brightness and kinematic moments are compared to observations. Any change to the given potential therefore results in re-integrating the orbit library to determine the orbital weights that best reproduce the observations via non-negative least-squares (NNLS) optimization. This can require running thousands of models to characterize a parameter space, and large orbit libraries in order to encapsulate the true diversity of orbital structures that may be present in galaxies. For the construction of our models, we set-up the gravitational potential, orbit libraries, and parameter space searches in a similar way to D24. We present an outline of our configuration process below.

The underlying gravitational potential of the galaxy is constructed with a triaxial stellar component, a spherical DM halo, and a central supermassive black hole. The triaxial stellar component mass is calculated from the best-fitting 2D MGE. Each of these MGEs mass densities is deprojected assuming the orientation in space of the galaxy, described by three viewing angles (θ, φ, ψ), to obtain a 3D mass density. The space orientation (θ, φ, ψ) can be converted directly to the intrinsic shape, described by the parameters $p = B/A$, $q = C/A$, and $u = A'/A$, where A, B, and C represent the major, medium, and minor axes of the 3D triaxial Gaussian, respectively, and A' is the projected major axis of the 3D triaxial Gaussian. We leave p, q, and u as free parameters to allow intrinsic triaxial shapes to be fitted.

A spherical Navarro–Frenk–White (NFW; Navarro, Frenk & White 1996) halo is adopted. The mass, $M_{200}$ (mass enclosed within a radius, $R_{200}$, where the average density is 200 times the critical density), in an NFW DM halo is determined by two parameters in our modelling: the concentration parameter, c, and the ratio of DM within $R_{200}$, $M_{200}/M_*$ (where $M_*$ is the total stellar mass). Due to the MAGPI and MAGPI mocks data quality, we adopt a $M_{200} - c$ coupled halo. This fixes the concentration according to the relation measured by Dutton & Macciò (2014), leaving the DM-to-stellar mass ratio within $R_{200}$, $M_{200}/M_*$, free.

For completeness, we also include a central supermassive black hole component. However, the spatial resolution of our data cannot resolve the influence radius of the black hole, so its mass leaves no imprint on the stellar kinematic maps of our galaxies and therefore does not affect our results. For this reason, similarly to D24, we include a black hole with fixed mass, derived by the redshift dependent mass–velocity dispersion relation from Robertson et al. (2006). The total (luminous plus dark) mass of the system can then be scaled by the global constant ($\Upsilon_*$), as the scaled model can be evaluated by scaling the orbital velocities without re-integrating the orbit library.

Our Schwarzschild implementation creates initial conditions for the orbits by sampling from the three integrals of motion (energy $E$, second integral $I_2$, and third integral $I_3$). Each set of orbits is sampled across the three integrals with the following number of points: $n_E \times n_{I_2} \times n_{I_3} = 21 \times 15 \times 9$. We test different samplings of the second integral $I_2$, for example, increasing $n_{I_2}$ from 15 to 18 and to 40, and refitting the models. The best-fitting values retrieved by this search are consistent, within the $1\sigma$ confidence level, with the best-fitting values retrieved by our regular runs. We use three sets

of 21 × 15 × 9 orbits: a set of box orbits, tube orbits, and counter-rotating tube orbits. Following D24, we also do not dither the orbits, since dithering is not necessary for MAGPI-like data (see appendix A of D24), nor do we regularize<fnref n="5" /> the orbital weights.

Once the orbit libraries are created, the orbit weights are determined by fitting the set of orbit-superposition models to the projected and deprojected luminosity density (from MGE fits) and the 2D line-of-sight stellar kinematics. The modelled and the observed kinematics are then divided by the observational error so that a $\chi^2$ comparison is achieved. The weights are determined by the van den Bosch et al. (2008) implementation, using the Lawson & Hanson (1974) NNLS implementation. The best-fitting model is defined as the model with minimum kinematic $\chi^2$. To ensure that we get the best deprojection fit, we first run what we call a ‘coarse’ grid of models across the allowed parameter space for the five free parameters, sampling a minimum of three evenly spaced values within the bounds for p, q, and u, at least six values in $\Upsilon_*$, and at least seven values in DM-to-stellar mass ratio $M_{200}/M_*$ between -2 and 3 in log space.

The best-fitting model from the coarse run is then used to initiate a dynamic ‘fine’ parameter search with increasingly smaller step sizes around the best-fitting solution, up to a minimum of 0.05 in log space for DM fraction, 0.02 for p and q, 0.01 for u (which generally tends very close to unity), and 0.1 in $\Upsilon_*$. The coarse and fine model runs use the ‘FullGrid’ and ‘LegacyGridSearch’ options, respectively, in the DYNAMITE code.

As in Zhu et al. (2018b) and D24, we define the $1\sigma$ confidence level as $\Delta\chi^2 = \sqrt{2N_{\text{kin}}N_{\text{mom}}}$, where $N_{\text{kin}}$ is the number of kinematic Voronoi bins and $N_{\text{mom}}$ is the number of kinematic Gauss–Hermite moments. We choose this definition to ease the comparison with previous work, and particularly with D24. The quoted parameter value (e.g. p) is calculated from the model with the lowest $\chi^2$ value. To estimate the uncertainties, parameter values from all $1\sigma$ models are calculated, and the upper and lower bounds of that distribution are taken. For example, the upper error is quoted as $\max(X) - x_0$, where X is the vector of all $1\sigma$ model parameter values and $x_0$ is the value from the lowest $\chi^2$ model. We show the best-fitting kinematic maps and parameter space coloured by Schwarzschild model $\chi^2$ value for example galaxy EAGLE 9498015 in Figs 3 and 4.

From these models we derive the internal orbital structure, inner mass distribution, intrinsic shape, and velocity anisotropy for each galaxy. Intrinsic shapes can also be reconverted into viewing angles. The viewing angles are typically well constrained, with a median uncertainty of $7^\circ$. In addition, following Zhu et al. (2018a), Jin et al. (2020), and Santucci et al. (2022, 2024), we separate orbits into four different components (cold, warm, hot, and counter-rotating) according to their orbit circularity $\lambda_z$ as described in Section 2.4, with $\lambda_z \equiv J_z/J_{\text{max}}(E)$ around the short z-axis, normalized by the maximum of a circular orbit with the same binding energy $E$. We calculate the fraction of orbits in each component within $1 R_*$.

We derive dynamically based intrinsic shapes using the triaxial parameter at $1 R_*$, $T_{\text{Re}} = (1 - p^2_{\text{Re}})/(1 - q^2_{\text{Re}})$, where $p_{\text{Re}}$ and $q_{\text{Re}}$ are the intermediate-to-long and short-to-long axis ratios at $1 R_*$. We separate galaxies into three groups according to this parameter: oblate ($T_{\text{Re}} < 0.3$), triaxial ($0.3 \le T_{\text{Re}} < 0.8$), and prolate ($T_{\text{Re}} \ge 0.8$).

<footnote n="5">⁵Lipka & Thomas (2021) presented a novel technique to optimize the amount of regularization used for orbit-weight solving, which can have a significant impact on constraining physical parameters such as galaxy mass, anisotropy, and viewing angle. We did not implement this technique due to the data being of poorer resolution than the one employed in Lipka & Thomas (2021), and hence adding regularization would not yield more precise results in our case.</footnote>

<footnote n="4">⁴Version: 3.0.0</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Example galaxy EAGLE 9498015. This galaxy has log $M_{*,50kpc}/M_{\odot} = 11.19$ and $R_e = 6.08$ kpc. Columns show 2D maps for, from left to right, flux, velocity, velocity dispersion, $h_3$ and $h_4$. The first row shows the observed maps, the second row shows the best-fitting maps derived from the Schwarzschild modelling, and the third row shows the residuals, calculated as the difference between the observation and the model, divided by the observational uncertainties (or divided by the model, in the case of the surface brightness). The best-fitting model maps ($\chi^2_{red} = 0.72$) accurately reconstruct the structures seen in the observations.</figcaption></figure>

A detailed analysis of these properties, and how they compare to the particle-derived properties is presented in the following section.

# 4 RESULTS

In this section, we present the results we obtain modelling our sample of 75 MAGPI-like EAGLE galaxies with the Schwarzschild modelling technique. First, we compare the properties derived from the models with the true values from the simulations in Section 4.1. We then present analysis of the intrinsic shapes and orbital structures of the sample, comparing our results with results from the MAGPI Survey (D24) in Section 4.2.

## 4.1 Comparison with particle data

### 4.1.1 Mass profiles

The total mass ($M_{\text{tot}}$) radial distribution is one of the fundamental parameters of the Schwarzschild model, which includes a stellar component and a DM component. As mentioned in Section 3.2, a black hole mass component is included as well, but not discussed here as its contribution to the total mass distribution is negligible. Fig. 5 shows the one-to-one comparison of the model recovered and the true value of the enclosed (within $1 R_e$) total mass (top panel), and DM mass (bottom panel).

In general, the total mass is recovered well, with a median offset of 9 per cent. However, we find that for almost half of the galaxies in our sample the retrieved value of DM within $1 R_e$ has a median offset of 48 per cent, in the majority of cases underestimating the true value. While there are degeneracies in the models that could lead to a not well-constrained DM profile, we find that this discrepancy is mostly due to the assumed NFW profile not being able to reproduce the DM distribution in the galaxy – likely because the DM distribution is not spherical, for example being more concentrated than an NFW. Moreover, for lower-mass galaxies, the models are constrained up to about $1 R_e$, which is not enough to well constrain the DM content.

Schaller et al. (2015) investigated the mass distributions of EAGLE central galaxies, fitting their halo densities with NFW profiles,

focusing on the outer regions of the haloes. We compare the enclosed DM mass from their NFW fits with our results, for the galaxies in both samples. In Fig. 6, we show the enclosed DM mass profiles from the particle data, the NFW fit on that data (from Schaller et al. 2015), and the enclosed DM mass retrieved by our modelling. The top panel shows the DM mass profiles for a galaxy with a DM distribution close to the NFW one, while the bottom panel shows the enclosed DM mass for a galaxy where the NFW profile is not a good fit in the inner region (although it is a good fit in the outer region of the halo, which it is not shown here). As seen, even though our fit does not reproduce the real DM distribution, it is consistent with the NFW fit on that galaxy. The retrieved value of DM mass within $1 R_e$ is therefore affected by the fact that the assumed NFW profile is not able to reproduce the DM distribution in the inner regions of the galaxy. Schaller et al. (2015) also showed that while the EAGLE haloes are well described by the NFW density profile at large radii, closer to the centre the presence of stars can produce steeper profiles. We believe these differences in the inner density profiles, as well as the possible non-spherical shape of the DM haloes, can lead to the discrepancies we see in the enclosed DM mass profiles.

### 4.1.2 Intrinsic shapes

We show the one-to-one comparison between the model-derived and the true axis ratios within $1 R_e$, $p_{Re}$, and $q_{Re}$, in Fig. 7. We find that the model-derived $p_{Re}$ is consistent, within the errors, with the true value, with a median overestimation of 0.04. The model-derived $q_{Re}$ is qualitatively consistent with the true value, but shows a large scatter around the one-to-one line, with a median overestimation of 0.18.

Fig. 8 shows the one-to-one comparison of the true and model-recovered triaxiality at $1 R_e$ of our galaxies. We find that the true value of triaxiality is recovered with a median offset of 0.07 (underestimating the true values). However, there is large scatter around the one-to-one line, mostly at the low and at the high end of the triaxiality parameter. Similar results were also found by Jin et al. (2019), looking at nine galaxies in the highest resolution simulation, Illustris-1, of the Illustris simulations. It is interesting to note that
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example galaxy EAGLE 9498015: model parameter grid. There are five free parameters: the intrinsic shape of the flattest Gaussian component ($p_{\text{min}}$, $q_{\text{min}}$, and $u_{\text{min}}$), the DM-to-stellar mass ratio, log $M_{200}/M_*$, and a global mass scaling factor, which scales both stellar and dark components jointly, $\Upsilon_*$. The diamonds represent the parameters explored, with the best-fitting model highlighted with a black cross. Models within the best-fitting region are colour-coded according to their $\chi^2$ values shown in the colour bar. The best-fitting values are well constrained. All 8336 models are shown, which is a typical number of models run for each EAGLE object.</figcaption></figure>

even though we find less agreement between the model-derived and the true $q_{\text{Re}}$, it is $p_{\text{Re}}$ that has the largest impact on $T_{\text{Re}}$, so that even a small difference (~ 0.02, e.g.) between the model-derived and the true $p_{\text{Re}}$ can lead to very different values of $T_{\text{Re}}$ (by construction, see equation 1). Moreover, the triaxiality parameter and the axis ratios are closely connected with the inclination angle of the observed galaxy – which can sometimes be difficult to constrain, in particular for galaxies where there is no clear disc structure nor triaxial features.

In order to understand the impact of the discrepancy between modelled and true shape, we separate galaxies into three groups according to their triaxiality: oblate ($T_{\text{Re}} < 0.3$), prolate ($T_{\text{Re}} \ge 0.8$), and triaxial ($0.3 \le T_{\text{Re}} < 0.8$). We find that 73 per cent of the galaxies in our sample that are classified as oblate, triaxial, or prolate using the true values, maintain the same classification when using the model-derived values. This shows that even though the exact values of individual galaxies could be off, in studying a bulk population we are able to get a statistical view of their shapes. Eight of the true oblate

galaxies (11 per cent of the total and 26 per cent of the true oblate galaxies) are classified as triaxial, 10 (14 per cent) of the true triaxial galaxies are classified as oblate, and two of the true prolate galaxies (3 per cent of the total) are classified as triaxial. In the majority of cases, the reason behind the wrong classification of the 18 oblate and triaxial misclassified galaxies is the low ($< 30^{\circ}$) inclination angle, which can be difficult to constrain in the modelling.

### 4.1.3 Fraction of orbits

As mentioned in Section 2.4, stellar orbits can be characterized by two main properties: the radius $r$ (in this case time-averaged), representing the size of each orbit, and their circularity $\lambda_z$. In the Schwarzschild models, these orbital distributions are time-averaged over 200 orbital periods. Taking the radius, $r$, and the circularity, $\lambda_z$, of each orbit, and considering their weights given by the solution from the best-fitting model, we can use the orbit circularity distribution in
<figure><img src="image_7.png" /><figcaption>**Figure 5.** One-to-one comparison of the true and model-recovered total (top panel) and dark (bottom panel) mass of EAGLE galaxies within 1 $R_e$. Fractional residuals [defined as (model − true)/true] of the total and DM mass, respectively, are shown in the bottom region of each panel. The one-to-one line is shown in grey. The model-reconstructed enclosed total mass well matches the true one, within the errors. The derived enclosed DM mass is in general in agreement with the true one, however there are several galaxies where the DM mass is underestimated.</figcaption></figure>

the phase space to obtain the probability density of orbits within 1 $R_e$, for each galaxy. An example of the circularity space can be seen in Fig. 9, for example, galaxy EAGLE 9498015, which is dominated by hot orbits. The panel on the left shows the particles’ instantaneous orbital distribution, while the middle panel shows the snapshot-averaged one and the left panel shows the orbital distribution from the models (time-averaged over 200 orbital periods). The fractions of mass in each orbital component within 1 $R_e$ are annotated in each panel. These figures show how the difficulty is not just in recovering

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Enclosed DM mass profiles from the particle data (in red), the NFW fit on that data (in dark orange, from Schaller et al. 2015), and the enclosed DM mass retrieved by our modelling (shaded area delimits the $1\sigma$ region). The top panel shows the DM mass profiles for a galaxy with a DM distribution close to the NFW one, while the bottom panel shows the enclosed DM mass for a galaxy where the NFW profile is not a good fit. As seen, even though our fit does not reproduce the real DM distribution, it is consistent with the NFW fit on that galaxy.</figcaption></figure>

the orbital distributions, but also in finding the best way to compare them between simulations and models, since the techniques used to derive them are different. Instantaneous orbital distributions show stellar particles covering almost the entire range in $\lambda_z$. As mentioned in Section 2.4, the stellar particles in our simulated galaxies do not necessarily conserve $\lambda_z$ when orbiting in the potential. This is particularly true for particles on radial/box orbits (which have, averaging over time, $\lambda_z \sim 0$). Instead, some of these particles will appear to have the same orbital circularity as particles on cold/warm
<figure><img src="image_9.png" /><figcaption>**Figure 7.** One-to-one comparison of the true and model-recovered axis ratios within $1R_e$, $p_{Re}$ (top panel), and $q_{Re}$ (bottom panel). We find that the model-derived $p_{Re}$ is consistent, within the errors, with the true value. $q_{Re}$ shows a large scatter around the one-to-one line, with a median overestimation of 0.18.</figcaption></figure>

orbits (i.e. 'spreading' the circularity distribution). This difference between the instantaneous (from the particle data) and the time-averaged (from the models) orbital distributions means that we cannot directly compare them (see Appendix A for more details). For this reason, we use the snapshot-averaged orbital distributions.

From the cuts defined in Section 3.2, we calculate the orbital fractions for the sample within $R_e$ and then compare them to their true values averaged over three snapshots. Fig. 10 (as described in Section 2.4) shows the one-to-one comparison of the true and model-recovered fraction of the cold, warm, and hot orbits.

We find that, for $\lambda_z$, the comparison between the model-derived and true fractions of cold and counter-rotating orbits are in general

<figure><img src="image_10.png" /><figcaption>**Figure 8.** One-to-one comparison of the true and model-recovered triaxial parameter $T_{Re} = (1 - p_{Re}^2)/(1 - q_{Re}^2)$. Galaxies with $T_{Re} < 0.3$ are classified as oblate, galaxies with $T_{Re} \ge 0.8$ as prolate and those in-between as triaxial. Grey dashed lines represent $T_{Re} = 0.1, 0.3$, and $0.8$. Galaxies that are classified as oblate using both their true and model-derived triaxiality are shown in orchid, galaxies classified as triaxial are shown in purple and galaxies classified as prolate are shown in blue. Galaxies that are misclassified using the model-derived triaxiality are shown in grey. In general, we find that, even when the values of triaxiality are not close to the one-to-one line, 73 per cent of the galaxies in our sample are correctly classified using their model-derived triaxiality.</figcaption></figure>

agreement, within the errors. The model-derived fractions of hot orbits are higher (especially at low fractions) than those derived from the particle data, while the warm fractions are generally lower (especially at high fractions). We find that the model-derived fractions of cold orbits have a median overestimation of 0.07, while the warm fractions are underestimated by 0.14, and the counter-rotating fractions of 0.02. The hot orbit fractions have a median overestimation of 0.06. We stress, however, that the discrepancies that we see are most likely not caused by the Schwarzschild models not being able to reproduce the orbital distributions well, but by the two different methods used to derive the orbits' circularity.

We also show the one-to-one comparison of the instantaneous and model-recovered fractions of the cold, warm, hot, and counter-rotating orbits in Appendix B (Fig. B1), for completeness. We find that the instantaneous orbits and the model-derived orbits have large scatter and offsets, compared to the snapshot-averaged values. In particular, as expected, the model-derived fractions of hot orbits are generally higher than the instantaneous fractions, while the model-derived fractions of warm orbits tend to be lower than the instantaneous fractions. As mentioned in Section 2.4, this is not surprising since the model-derived orbital distributions are time-averaged.

## 4.2 Dark matter content, intrinsic shapes, and orbital distributions of massive galaxies in EAGLE and MAGPI

We now explore the distribution of the model-derived galaxy properties (DM fractions, intrinsic shapes, fraction of orbits, and velocity anisotropy) with stellar mass and as a function of other galaxy
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Example galaxy EAGLE 9498015: orbit circularity $\lambda_z$ versus radius (in kpc). The colour indicates the density, w, of the orbits on the phase space, and the three horizontal black dashed lines indicate $\lambda_z = 0.8, 0.25$, and $-0.25$, dividing the orbits into four regions (cold, warm, hot, and counter-rotating orbits). Fractions of mass within $1R_e$ in each component are annotated in the figures. The vertical grey dashed line is located at $1R_e$. Left panel shows the true instantaneous orbital distribution of the stellar particles; middle panel show the average orbital distribution calculated over three snapshots; and right panel shows the time-averaged orbital distribution derived from the dynamical modelling. Averaging over the three snapshots helps to reconstruct the orbital distribution and makes it possible to compare the particles’ circularity with what is derived by DYNAMITE.</figcaption></figure>

properties, and we compare our results with those derived for 22 MAGPI galaxies by D24.

We show in Fig. 11 the distributions in the size–mass plane for our sample of EAGLE galaxies, in light blue, and the MAGPI galaxies in D24, in grey. EAGLE galaxies are larger than the galaxies in the MAGPI sample. This is likely a result of our Voronoi binning selection criterion, which excluded the most compact galaxies from our sample. The two samples also differ in stellar mass: our EAGLE sample has relatively more galaxies in the stellar mass range $10.8 < \log M_{*,50kpc}/M_{\odot} < 11.2$ than the MAGPI sample.

We define the fraction of DM within $1R_e$, $f_{DM}$, as the ratio between the enclosed mass of DM within $1R_e$ and the enclosed total mass within $1R_e$. We show in Fig. 12, the model-derived DM fractions for our 75 EAGLE galaxies (in light blue) and for the 22 MAGPI galaxies (in grey) modelled by D24. To better compare the results, we also show the mean values of $f_{DM}$ for different bins in mass. We find that the median DM fraction within the $1R_e$ is 26 per cent, with a standard deviation of 15 per cent for the EAGLE sample. This fraction is larger, although still consistent within the errors, than what is found for MAGPI galaxies ($f_{DM} = 10$ per cent with a standard deviation of 19 per cent). This discrepancy is mainly driven by the large fractions of DM found for EAGLE galaxies with $\log(M_{*,50kpc}/M_{\odot}) < 11.2$. For galaxies with stellar masses above $\log(M_{*,50kpc}/M_{\odot}) \sim 11.2$, we find comparable DM fractions for EAGLE and MAGPI galaxies of similar stellar masses ($f_{DM} = 26$ per cent with a standard deviation of 9 per cent and $f_{DM} = 27$ per cent with a standard deviation of 21 per cent, respectively). We note, however, that the uncertainties in these fractions are large, and they are not taken in consideration when calculating the median values of DM in the different mass bins in Fig. 12.

In Fig. 13, we present the intrinsic shapes as the distribution of axis ratios $q_{Re}$ versus $p_{Re}$ for all galaxies in the EAGLE (light blue circles) and MAGPI (grey squares) samples. The grey dashed curves separate the galaxies into oblate ($T_{Re} < 0.3$), triaxial ($0.3 \le T_{Re} < 0.8$), and prolate ($T_{Re} \ge 0.8$) galaxies. We find that the majority of galaxies are consistent with being ellipsoids of varying degrees of triaxiality, with a few tending towards circular (i.e. highly flattened, oblate) discs in their intrinsic shapes. Only 13/75 EAGLE galaxies are classically defined oblate spheroids, with a triaxiality parameter of less than 0.1. The rest of the galaxies in the sample are oblate (19/75) and triaxial (35/75), with seven prolate objects. MAGPI galaxies show different distributions in both $p_{Re}$ and $q_{Re}$, with values generally higher, in

particular for $q_{Re}$, than EAGLE galaxies. Matching the samples in both mass range and size does not decrease the differences. However, we find that, despite these differences, MAGPI galaxies show a similar range in triaxiality as the EAGLE sample. The two samples have similar percentages of oblate spheroids: 18 per cent $\pm 4$ per cent for EAGLE and 14 per cent $\pm 7$ per cent for MAGPI. The MAGPI sample has a higher fraction (77 per cent $\pm 9$ per cent) of non-oblate ($T_{Re} \ge 0.3$) galaxies than EAGLE (57 per cent $\pm 6$ per cent).

We also check whether we find a similar close connection for EAGLE galaxies between triaxiality and the fraction of hot orbits as found for MAGPI galaxies (D24). We show the triaxiality as a function of the fraction of orbits within $1R_e$ for both EAGLE and MAGPI galaxies in Fig. 14. We find that galaxies with higher triaxiality tend to have higher fractions of hot orbits and lower fractions of warm and cold orbits, while we do not see any trend with the counter-rotating orbits. To test whether these trends are statistically significant, we use the Kendall’s correlation coefficient $\tau$, using the PYTHON package `scipy.stats.kendalltau` (Virtanen et al. 2020). This correlation coefficient is robust to small sample sizes. A $\tau$ value close to 1 indicates strong correlation, whereas a value close to -1 indicates strong anticorrelation. For the fractions of hot orbits, we find a value of $\tau = 0.36$, with a p-value of $4.29 \times 10^{-6}$. Negative correlations are found for the fraction of cold and warm orbits, with $\tau = -0.28$ and p-value of $4.57 \times 10^{-4}$ and $\tau = -0.41$ and p-value of $2.16 \times 10^{-7}$, respectively.

D24 found in their MAGPI sample that radially anisotropic ($\beta_r > 0$) objects are likely to be intrinsically triaxial and have mainly hot (box or radial) orbits, with few cold (short-axis tube) orbits. To check whether this is also true in our simulation sample, we define the velocity anisotropy parameter $\beta_v$, in the standard spherical coordinates $(r, \theta, \phi)$, following Binney & Tremaine (2008):

$$ \beta_r = 1 - \frac{\Pi_{tt}}{\Pi_{rr}}, \qquad (2) $$

with

$$ \Pi_{tt} = \frac{\Pi_{\theta\theta} + \Pi_{\phi\phi}}{2}, \qquad (3) $$

$$ \Pi_{kk} = \int \rho \sigma_k^2 d^3 x = \sum_{n=1}^{N} M_n \sigma_{k,n}^2 \qquad (4) $$

with $\sigma_k$ the velocity dispersion along the direction k at a given location inside the galaxy. The summation defines how we computed
<figure><img src="image_12.png" /><figcaption>**Figure 10.** One-to-one comparison of the true and model-recovered fraction of the cold (in blue), warm (in orange), hot (in red), and counter-rotating (in green) orbits. The true values are calculated from the average of the orbital distribution in three different snapshots. We find that, for $\lambda_z$, the comparison between the model-derived and true fractions of orbits are generally in agreement, within the errors. However, fractions of hot orbits tend to be overestimated and fractions of warm orbits tend to be underestimated by the models.</figcaption></figure>

this quantity from our Schwarzschild models. $M_n$ is the mass contained in each of the $N$ polar grid cells in the meridional plane of the model, and $\sigma_{k,n}$ is the corresponding mean velocity dispersion along the direction $k$. We calculate the value of $\beta_r$ at $1 R_e$, and we denote it $\beta_{Re}$. $\beta_{Re} > 0$ indicates radial anisotropy, $\beta_{Re} < 0$ indicates tangential anisotropy and $\beta_{Re} = 0$ indicates isotropy.

We show the velocity anisotropy $\beta_{Re}$ as a function of the fraction of hot orbits, for EAGLE (light blue circles) and MAGPI (grey squares) galaxies in Fig. 15. We find a correlation between the fraction of hot orbits and $\beta_{Re}$ with $\tau = 0.38$, and $p$-value $= 2.04 \times 10^{-6}$, so that galaxies with higher fractions of hot orbits are more likely to be

close to isotropy ($\beta_{Re} = 0$) or radially anisotropic ($\beta_{Re} > 0$), similar to what was found for MAGPI galaxies by D24 ($\tau = 0.47$, and $p$-value $= 1.71 \times 10^{-3}$).

# 5 DISCUSSION

In general, we find that Schwarzschild models of mock observations from EAGLE can, in most cases, reconstruct the true values of each parameter quite accurately, within the errors. However, there are parameters where the retrieved values are offset. One of these is the enclosed DM mass. We find that, while there are degeneracies in the
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Distribution of the EAGLE galaxies in our sample (light blue circles) and of the MAGPI galaxies in D24 (grey squares) in the size-stellar mass plane. We show the marginalized distributions (calculated using kernel density estimates) of mass and size of the galaxies in each sample in the top and right panels of the figure, respectively. Light blue lines are for the EAGLE galaxies in our sample, while grey lines are for the galaxies in the MAGPI sample. EAGLE galaxies are larger than the galaxies in the MAGPI sample. The two samples also differ in stellar mass: our EAGLE sample has relatively more galaxies in the stellar mass range $10.8 < \log M_{\star,50\,\text{kpc}}/M_{\odot} < 11.2$ than the MAGPI sample.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Fraction of DM to total mass ($f_{\text{DM}} = M_{\text{dark}}/M_{\text{tot}}$) within $1 R_e$ from dynamical models as a function of stellar mass for EAGLE galaxies (light blue circles) and MAGPI galaxies (grey squares). Bold points show the mean values for each mass bin, and the error bars show their $1\sigma$ scatter. We find a median DM fraction of 26 per cent, with a standard deviation of 15 per cent. EAGLE galaxies with stellar masses below $\log(M_{\star,50\,\text{kpc}}/M_{\odot}) \sim 11.2$ show higher DM fractions than MAGPI galaxies, while EAGLE galaxies at higher stellar masses have values of DM fractions comparable to those of MAGPI galaxies.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** The distribution of axis ratios $q_{\text{Re}}$ versus $p_{\text{Re}}$ at $1 R_e$ from the best-fitting models for all galaxies in our EAGLE sample (light blue circles) and in the MAGPI sample (grey squares; D24). The grey dashed curves from right to left indicate $T_{\text{Re}} = 0.3$ and 07; the orange line indicates $T_{\text{Re}} = 1.0$, and the red line is for $T_{\text{Re}} = 0.0$. Perfectly oblate spheroids have zero triaxiality. The blue and grey daggers in the top left corner show the median (asymmetric) errors on the intrinsic axial ratios for EAGLE and MAGPI, respectively. We show the marginalized distributions (calculated using kernel density estimates) of $q_{\text{Re}}$ and $p_{\text{Re}}$ of the galaxies in each sample in the top and right panels of the figure, respectively. Light blue lines are for the EAGLE galaxies in our sample, while grey lines are for the galaxies in the MAGPI sample. The majority of the galaxies in both the EAGLE sample and the MAGPI sample have non-oblate shapes.</figcaption></figure>

models for which the enclosed DM mass can be harder to constrain, the majority of the discrepancy can be explained by the fact that we have assumed a spherical NFW distribution in the Schwarzschild modelling which, in some cases, it is not representative of the DM distribution in the simulated galaxies. We find a median underestimation of the enclosed DM mass of 48 per cent. It is interesting to note that previous results for the Illustris simulation (Jin et al. 2019) have found that the model-derived DM within $1 R_e$ was generally overestimated (as opposed to the underestimation we find), while they found the stellar mass to be underestimated. This difference is most likely due to the fact that we are using stellar mass maps to construct our MGEs, instead of having light-based MGEs (as in Jin et al. 2019), helping to constrain the stellar mass. Moreover, Valenzuela et al. (2024) found that the shape of the DM haloes in the MAGNETICUM hydrodynamical cosmological simulation and found that, in the inner region of galaxies, the DM follows the stellar component in shape and orientation. Since DYNAMITE allows for different types of DM distributions to be assumed, it would be interesting to test the systematics each assumption causes and their effect on the enclosed mass distributions in future work.

We also find that the derived intrinsic shapes within $1 R_e$ have large uncertainties, especially for galaxies with face-on viewing angles. 27 per cent of the galaxies in our sample are misclassified, most of them being oblate galaxies classified as triaxial, or triaxial galaxies classified as prolate. One possible reason for this is that a small variation in the axis ratios $p$ and $q$ can lead to very different values of $T_{\text{Re}}$ for nearly spherical systems (due to how $T_{\text{Re}}$ is defined). In
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Triaxiality, $T_{Re}$ as a function of the fraction of orbits within 1 $R_e$ for both EAGLE (coloured circles) and MAGPI (grey squares) galaxies. We find that the triaxiality parameter increases with increasing fraction of hot orbits and decreasing fraction of cold and warm orbits. No clear correlation is found between the fraction of counter-rotating orbits and a galaxy's triaxiality.</figcaption></figure>

In particular, if the reconstructed values of *p* are different from the true ones, even though still consistent within the errors, they can lead to a large change in $T_{Re}$. Moreover, the range in the shape parameters that we cover is determined by the minimum value of the axis ratio from the projected surface brightness/mass in the MGEs. For face-on views, it is very possible for the modelling to overestimate *q*, as seen in Fig. 7. In addition to that, as mentioned in the previous paragraphs, we take the DM component to be spherical, while we assume a triaxial stellar component. This means that the degeneracy between stellar and DM mass could influence the intrinsic shapes. Even though a high

percentage (27 per cent) of galaxies are misclassified, we note that $T_{Re}$ never changes as much as to have prolate galaxies misclassified as oblate galaxies or oblate galaxies as prolate.

For the fraction of orbits within 1 $R_e$, we find that the cold, hot, and counter-rotating mass fractions are reconstructed qualitatively well, with median offset of 0.07, 0.06, and 0.02, respectively. We find that the model-derived fractions of the warm orbits are underestimated by 0.14. These results are in agreement with Jin et al. (2019). The differences in the orbital distributions, however, can be at least partially attributed to a difference in measurement techniques. Even
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Velocity dispersion anisotropy in spherical coordinates within 1 $R_e$, $\beta_{Re}$, as a function of the fraction of hot orbits for EAGLE galaxies (light blue circles) and MAGPI galaxies (grey squares). $\beta_{Re} > 0$ indicates radial anisotropy, $\beta_{Re} < 0$ indicates tangential anisotropy and $\beta_{Re} = 0$ indicates isotropy. We find that $\beta_{Re}$ increases with increasing fraction of hot orbits, so that galaxies with higher fraction of hot orbits are more radially anisotropic. We also find a good agreement between the MAGPI and EAGLE samples.</figcaption></figure>

though we have averaged the instantaneous orbital distribution from the particle data using three different snapshots, these distributions are not the same as the time-averaged distributions we get from the models (which are calculated, by default, over 200 orbital periods, meaning that, especially at larger radii, the orbital distributions can be averaged over a *Hubble* time). Moreover, galaxy interactions or interactions with the environment can have an effect on the orbital distributions, and can lead to snapshot-averaged values being very different to what we would have if we calculated the circularity by freezing the potential, integrating the particle orbits in the potential, and then calculating the average values along the orbits (which is, however, computationally and time expensive).

After checking the comparison between the model-derived properties and the true ones, we compare our results with the orbital properties of MAGPI galaxies derived using Schwarzschild dynamical modelling by D24. We find that our selected sample from the EAGLE galaxies reproduces the galaxy properties derived from observed MAGPI galaxies qualitatively well, showing similar trends in orbital fractions, triaxiality (although not individual axis ratios), and velocity anisotropy, although with a larger scatter in all of the relations explored. EAGLE galaxies, however, seem to have higher fractions of DM, compared to their MAGPI counterparts, and are less round in shape. The galaxies in the EAGLE and MAGPI samples are mostly non-oblate objects, which are likely to be radially anisotropic ($\beta_{Re} > 0$), and to have mainly hot (box or radial) orbits, with few cold (short-axis tube) orbits.

To understand the connections between the different galaxy properties we analyse here and to understand the correlations we find for the galaxies in the EAGLE and MAGPI samples, we will now explore the connection between the merger history of our EAGLE galaxies and their intrinsic properties at $z \sim 0.3$.

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Medians (symbols connected by lines) of the fractions of each orbital component, as labelled, within 1 $R_e$, as a function of the total number of minor and/or major mergers experienced by the galaxies in our sample from $z=2$ up to 0.27. Error bars show the 1$\sigma$ scatter around the median. Note that we perturb the x-axis position of the symbols slightly to avoid overlapping error bars. The fractions of hot orbits increase with increasing number of mergers, while the fractions of warm and cold orbits decrease with increasing number of mergers.</figcaption></figure>

## 5.1 Intrinsic properties and merger history

We define mergers and merger histories in the same way as in Santucci et al. (2024), but based on the merger history of galaxies at $z=0.27$: galaxies with ‘no mergers’ are those that did not have a minor or major merger since $z=2$, with minor and major mergers being those with a stellar mass ratio between the second most massive and most massive progenitors between [0.1–0.3] and [0.3,1], respectively.

The majority of the galaxies in our sample ($\sim 69$ per cent) have experienced at least one merger (either minor or major) in their life, and 31 per cent of the galaxies have experienced two or more mergers, with the majority of these being minor mergers. In order to constrain the correlation between galaxy mergers and a galaxy’s orbital properties, we show in Fig. 16 the median values of the orbital fractions within 1 $R_e$ as a function of the total number of minor and/or major mergers each galaxy has undergone from $z=2$ up to 0.27. We do not find any correlation between the fractions of counter-rotating orbits and the number of mergers. Even though the overall trends are consistent with constant, we see suggestions of possible correlations between the other fraction of orbits and the number of mergers: the fractions of hot orbits increase with increasing number of mergers, while the fractions of cold and warm orbits decrease, with galaxies that have experienced three or more mergers having $\sim 51$ per cent more mass in hot orbits, $\sim 57$ per cent less mass in warm orbits, and 66 per cent less mass in cold orbits than galaxies with $\le 1$ merger.

Several theories have been proposed where different types of mergers can have different effects on a galaxy’s structure and composition, leading to the formation of slow- or fast-rotating galaxies. For example, Cappellari (2016) described how early-type galaxies can form through two main formation channels: fast-rotating galaxies start as star-forming discs and evolve through a set of
processes dominated by gas accretion, bulge growth, and quenching. By comparison, slowly rotating galaxies assemble near the centre of massive haloes via intense star formation at high redshift (Naab, Johansson & Ostriker 2009; Clauwens et al. 2018), and evolve from a set of processes dominated by gas-poor mergers, resulting in more triaxial shapes. Gas-poor and gas-rich mergers can also have very different effects on the spin of galaxies (as shown by Naab et al. 2014, e.g.). More recently, Schulze et al. (2018) studied galaxies in the MAGNETICUM simulation and found that at least 30 per cent of the population of slow rotators at $z \sim 0$ underwent the transition from fast to slow rotator after a significant merging event (with a mass ratio between 0.2 and 1). Lagos et al. (2022) used mock observations of galaxies with $M_* > 10^{10} M_\odot$ from the EAGLE simulation and found that slow-rotating galaxies with triaxial or prolate shapes are usually formed as a consequence of major and minor mergers, with prolate shapes being driven by gas-poor mergers. Slow-rotating galaxies that only had very minor mergers are associated with triaxial shapes while those that formed in the absence of mergers are oblate.

To explore how orbital fractions are influenced by different types of mergers, we classify the mergers into gas-rich (wet) and gas-poor (dry) using the ratio of gas to stellar mass of the merger ($f_{gas,merger}$, derived by Lagos et al. 2018). The latter is the sum of the gas masses of the two merging galaxies divided by the sum of the stellar masses of the two galaxies. In this paper, we classify mergers with a gas fraction below 0.2 as dry or gas-poor, and mergers with gas fractions above this value as wet or gas-rich (as in Santucci et al. 2024). We find that the majority of the mergers experienced by the galaxies in our sample were wet mergers, with ~35 per cent of the galaxies experiencing a dry merger – of these, 1/3 had undergone one major dry merger. We show the median values of the orbital fractions within 1 $R_e$ as a function of the number of wet (top panel) and dry (bottom panel) minor and/or major mergers each galaxy has undergone from $z=2$ up to 0.27 in Fig. 17. Note that each panel does not include galaxies that have experienced both wet and dry mergers. We find that the change in the fractions of hot and warm orbits that we see with number of mergers is driven by dry mergers, while no change in the fractions of orbits is found with increasing number of wet mergers.

This scenario and the changes in orbital fractions we find with number of mergers are consistent with the correlation we found between fraction of orbits and number of wet/dry mergers in Santucci et al. (2024), where the effects of wet mergers on the orbital fractions are weak, but appear to accumulate with the number of wet mergers (so that no clear difference is found after one wet merger), while the effects of dry mergers are stronger, so that even one dry merger can increase the fraction of hot orbits and decrease the fractions of warm and cold orbits.

The high fraction of wet mergers can also explain the high number of triaxial galaxies and low number of prolate galaxies that we see in our sample. Previous studies found prolate systems with minor-axis rotation produced through collisionless equal-mass merger of disc galaxies (Moody et al. 2014; Rodriguez-Gomez et al. 2015). Li et al. (2018) found that almost all prolate galaxies in the Illustris simulation have experienced at least one dry merger (mainly major), with a radial merger orbit, and they find nearly no prolate galaxy with stellar mass $M_* < 3 \times 10^{11} M_\odot$. On the other hand, Valenzuela et al. (2024) found that this relation is not as clear for galaxies in the MAGNETICUM simulation. Lagos et al. (2022) found, for galaxies in the EAGLE simulation, that prolate shapes are driven by gas-poor mergers. Li et al. (2018) also found that the majority of the galaxies they studied in the Illustris simulations had triaxial shapes between $z = 2.32$ and 0.46, due to the high number of wet mergers, which can often make

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Medians (symbols connected by lines) of the fractions of each orbital component, as labelled, within 1 $R_e$, as a function of the number of minor and/or major wet (top panel) and dry (bottom panel) mergers experienced by the galaxies in our sample from $z=2$ up to 0.27. Error bars show the $1\sigma$ scatter around the median. Note that we perturb the x-axis position of the symbols slightly to avoid overlapping error bars. The fractions of hot orbits increase with increasing number of dry mergers, while the fractions of warm orbits decrease with increasing number of dry mergers. No change is found in the fraction of orbits as a function of wet mergers.</figcaption></figure>

the galaxy shapes irregular. Pulsoni et al. (2020) studied the shape of the stellar haloes of galaxies from the IllustrisTNG simulations and found a non-negligible fraction of galaxies with prolate shapes, in particular at $\log M_*/M_\odot > 11$, with the triaxiality increasing with increasing radius. They did not find any clear difference between the shapes of the stellar haloes (when looking at the outskirts) of fast- and slow-rotating galaxies.
For the future, it will be interesting to compare the galaxy populations at $z \sim 0.3$ and $\sim 0$, to understand how galaxy shapes evolve with time, by comparing our results with observations, such as SAMI (Santucci et al. 2022) and ATLAS3D (Thater et al. 2023a), with shapes derived using triaxial Schwarzschild models.

Santucci et al. (2024) found good agreement between the orbital distributions of SAMI and EAGLE galaxies at $z \sim 0$. The fact that we also find good agreement between the orbital distribution of MAGPI and EAGLE galaxies is interesting on its own because it shows that, at least up to $z \sim 0.3$, the EAGLE simulations seem to reproduce the dynamical properties of the galaxies in the observed Universe qualitatively well.

# 6 SUMMARY

We build Schwarzschild dynamical models for 75 MAGPI-like mock observations of EAGLE galaxies. We use the state-of-the-art DYNAMITE code (Jethwa et al. 2020), allowing for fully triaxial Schwarzschild models, and make no assumptions on the orbital structure. We derive their intrinsic shapes and orbital distributions and compare them with the true values from the particle data. We then explore the distribution of the model-derived galaxy properties (DM fractions, intrinsic shapes, fraction of orbits, and velocity anisotropy) with stellar mass and as a function of other galaxy properties, and we compare our results with those derived for 22 MAGPI galaxies by D24.

For our EAGLE sample, we find that:

(i) the enclosed DM mass within $R_e$ is well reproduced when the DM has a distribution that is well represented by an NFW profile. However, enclosed DM can be underestimated in the other cases, with a median offset of 48 per cent;

(ii) 73 per cent of the galaxies in our sample have model-derived shapes consistent with their intrinsic ones, with a median offset in triaxiality of 0.07. Even though the exact values of individual galaxies can be off, in studying a bulk population we are able to get a statistical view of their shapes;

(iii) the cold, hot, and counter-rotating mass fractions are reconstructed qualitatively well, with median offsets of 0.07, 0.06, and 0.02, respectively, while the model-derived fractions of warm orbits have a median offset of 0.14;

(iv) our selected sample from the EAGLE galaxies tends to reproduce the trends found from observed MAGPI galaxies (by D24), showing similar trends in orbital fractions, triaxiality and velocity anisotropy, although with a larger scatter in all of the relations explored. EAGLE galaxies, however, seem to have higher fractions of DM, compared to their MAGPI counterparts, in particular at $\log M_*/M_{\odot} < 11.2$. The galaxies in the EAGLE and MAGPI samples are mostly non-oblate ($T_{Re} \ge 0.3$) objects, which are likely to be radially anisotropic ($\beta_{Re} > 0$), and to have mainly hot (box or radial) orbits, with few cold (short-axis tube) orbits;

(v) over half of the galaxies in our sample ($\sim 69$ per cent) have experienced at least one merger (either minor or major) in their life. This seems to lead to an increase in the fractions of hot orbits with increasing number of mergers, and a slight decrease in the fractions of cold and warm orbits with increasing number of mergers. The effect, however, is cumulative and seems to have some impact only on the orbital fractions of galaxies with two or more mergers;

(vi) the majority of the mergers experienced by the galaxies in our sample were wet mergers – with 35 per cent of the galaxies in our sample experiencing a dry merger – explaining the weak, but cumulative effects of mergers on the fraction of orbits.

Given the array of simulations present in the MAGPI theoretical library, it would be interesting to see how our results compare to properties derived from mock observations of different simulations. This comparison could also help to understand whether the differences we see between simulations and observations come from the different ‘recipes’ used to create them, or from the underlying physics.

# ACKNOWLEDGEMENTS

We warmly thank the reviewer for their comments which helped improve this work. We wish to thank the ESO staff, and in particular the staff at Paranal Observatory, for carrying out the MAGPI observations. MAGPI targets were selected from GAMA. GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. GAMA was funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. GAMA photometry is based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 179.A-2004 and 177.A-3016. The MAGPI team acknowledge support by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project no. CE170100013. This work was performed on the OzSTAR national facility at Swinburne University of Technology. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government, and from the Victorian Higher Education State Investment Fund (VHESIF) provided by the Victorian Government. GS thanks Luca Cortese and Bodo Ziegler for their suggestions, which improved the manuscript. GS and KH acknowledge funding from the Australian Research Council (ARC) Discovery project DP210101945. GS acknowledges funding from the UWA Research Collaboration Awards 2022. CL has received funding from the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project no. CE170100013. AP acknowledges support from the Hintze Family Charitable Foundation. AFM acknowledges support from RYC2021-031099-I and PID2021-123313NA-I00 of MICIN/AEI/10.13039/501100011033/FEDER, UE, NextGenerationEU/PRT. CL, JTM and CF are the recipients of ARC Discovery Project DP210101945. LMV acknowledges support by the German Academic Scholarship Foundation (Studienstiftung des deutschen Volkes) and the Marianne-Plehn-Program of the Elite Network of Bavaria. SMS acknowledges funding from the Australian Research Council (DE220100003).

# DATA AVAILABILITY

The data underlying this article are available through http://icc.dur.ac.uk/Eagle. Data products created for this article are available upon request to the corresponding author.

Measurements from Schwarzschild models of MAGPI galaxies are available by contacting Richard McDermid.

# REFERENCES

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

Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781

Bassett R., Foster C., 2019, MNRAS, 487, 2354

Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton University Press, Princeton, NJ

Brough S. et al., 2017, ApJ, 844, 59
Bryant J. J. et al., 2015, MNRAS, 447, 2857

Bundy K. et al., 2015, ApJ, 798, 7

Cappellari M., 2002, MNRAS, 333, 400

Cappellari M., 2016, ARA&A, 54, 597

Cappellari M., Copin Y., 2003, MNRAS, 342, 345

Cappellari M. et al., 2007, MNRAS, 379, 418

Cappellari M. et al., 2011, MNRAS, 413, 813

Choi H., Yi S. K., Dubois Y., Kimm T., Devriendt J. E. G., Pichon C., 2018, ApJ, 856, 114

Clauwens B., Schaye J., Franx M., Bower R. G., 2018, MNRAS, 478, 3994

Crain R. A. et al., 2015, MNRAS, 450, 1937

Cretton N., de Zeeuw P. T., van der Marel R. P., Rix H.-W., 1999, ApJS, 124, 383

Croom S. M. et al., 2012, MNRAS, 421, 872

Croom S. M. et al., 2024, MNRAS, 529, 3446

Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140

de Nicola S., Thomas J., Saglia R. P., Snigula J., Kluge M., Bender R., 2024, MNRAS, 530, 1035

Derkenne C. et al., 2023, MNRAS, 522, 3602

Derkenne C. et al., 2024, MNRAS, 533, 1300

Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497

Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948

Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359

Emsellem E., Monnet G., Bacon R., 1994, A&A, 285, 723

Emsellem E. et al., 2007, MNRAS, 379, 401

Emsellem E. et al., 2011, MNRAS, 414, 888

Fahrion K. et al., 2019, A&A, 628, A92

Falcón-Barroso J., Sánchez-Blázquez P., Vazdekis A., Ricciardelli E., Cardiel N., Cenarro A. J., Gorgas J., Peletier R. F., 2011, A&A, 532, A95

Feldmeier-Krause A., Zhu L., Neumayer N., van de Ven G., de Zeeuw P. T., Schödel R., 2017, MNRAS, 466, 4040

Foster C. et al., 2021, Publ. Astron. Soc. Aust., 38, e031

Gebhardt K. et al., 2003, ApJ, 583, 92

Gerhard O. E., 1993, MNRAS, 265, 213

Harborne K. E., van de Sande J., Cortese L., Power C., Robotham A. S. G., Lagos C. D. P., Croom S., 2020, MNRAS, 497, 2018

Harborne K. E. et al., 2023, Publ. Astron. Soc. Aust., 40, e048

Jethwa P., Thater S., Maindl T., Van de Ven G., 2020, Astrophysics Source Code Library, record ascl:2011.007

Jiménez E., Lagos C. d. P., Ludlow A. D., Wisnioski E., 2023, MNRAS, 524, 4346

Jin Y., Zhu L., Long R. J., Mao S., Xu D., Li H., van de Ven G., 2019, MNRAS, 486, 4753

Jin Y., Zhu L., Long R. J., Mao S., Wang L., van de Ven G., 2020, MNRAS, 491, 1690

Jin Y., Zhu L., Zibetti S., Costantin L., van de Ven G., Mao S., 2024, A&A, 681, A95

Krajnović D., McDermid R. M., Cappellari M., Davies R. L., 2009, MNRAS, 399, 1839

Lagos C. d. P., Theuns T., Stevens A. R. H., Cortese L., Padilla N. D., Davis T. A., Contreras S., Croton D., 2017, MNRAS, 464, 3850

Lagos C. d. P. et al., 2018, MNRAS, 473, 4956

Lagos C. d. P., Emsellem E., van de Sande J., Harborne K. E., Cortese L., Davison T., Foster C., Wright R. J., 2022, MNRAS, 509, 4372

Lawson C. L., Hanson R. J., 1974, Solving Least Squares Problems. Englewood Cliffs, Prentice-Hall

Li H., Mao S., Emsellem E., Xu D., Springel V., Krajnović D., 2018, MNRAS, 473, 1489

Liepold C. M., Quenneville M. E., Ma C.-P., Walsh J. L., McConnell N. J., Greene J. E., Blakeslee J. P., 2020, ApJ, 891, 4

Lipka M., Thomas J., 2021, MNRAS, 504, 4599

Ludlow A. D., Fall S. M., Wilkinson M. J., Schaye J., Obreschkow D., 2023, MNRAS, 525, 5614

Ma C.-P., Greene J. E., McConnell N., Janish R., Blakeslee J. P., Thomas J., Murphy J. D., 2014, ApJ, 795, 158

Mehrgan K., Thomas J., Saglia R., Parikh T., Neureiter B., Erwin P., Bender R., 2024, ApJ, 961, 127

Moody C. E., Romanowsky A. J., Cox T. J., Novak G. S., Primack J. R., 2014, MNRAS, 444, 1475

Naab T., Johansson P. H., Ostriker J. P., 2009, ApJ, 699, L178

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

Naiman J. P. et al., 2018, MNRAS, 477, 1206

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

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

Neureiter B. et al., 2021, MNRAS, 500, 1437

Neureiter B., Thomas J., Rantala A., Naab T., Mehrgan K., Saglia R., de Nicola S., Bender R., 2023, ApJ, 950, 15

Peng Y.-j. et al., 2010, ApJ, 721, 193

Pilawa J. D., Liepold C. M., Delgado Andrade S. C., Walsh J. L., Ma C.-P., Quenneville M. E., Greene J. E., Blakeslee J. P., 2022, ApJ, 928, 178

Pillepich A. et al., 2018, MNRAS, 475, 648

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

Poci A., McDermid R. M., Zhu L., van de Ven G., 2019, MNRAS, 487, 3776

Poci A. et al., 2021, A&A, 647, A145

Pulsoni C., Gerhard O., Arnaboldi M., Pillepich A., Nelson D., Hernquist L., Springel V., 2020, A&A, 641, A60

Quenneville M. E., Liepold C. M., Ma C.-P., 2021, ApJS, 254, 25

Quenneville M. E., Liepold C. M., Ma C.-P., 2022, ApJ, 926, 30

Robertson B., Hernquist L., Cox T. J., Di Matteo T., Hopkins P. F., Martini P., Springel V., 2006, ApJ, 641, 90

Robotham A. S. G., Davies L. J. M., Driver S. P., Koushan S., Taranu D. S., Casura S., Liske J., 2018, MNRAS, 476, 3137

Rodriguez-Gomez V. et al., 2015, MNRAS, 449, 49

Rosas-Guevara Y. M. et al., 2015, MNRAS, 454, 1038

Rusli S. P. et al., 2013, AJ, 146, 45

Sánchez S. F. et al., 2012, A&A, 538, A8

Santucci G. et al., 2022, ApJ, 930, 153

Santucci G. et al., 2023, MNRAS, 521, 2671

Santucci G. et al., 2024, MNRAS, 528, 2326

Schaller M. et al., 2015, MNRAS, 451, 1247

Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210

Schaye J. et al., 2015, MNRAS, 446, 521

Schulze F., Remus R.-S., Dolag K., Burkert A., Emsellem E., van de Ven G., 2018, MNRAS, 480, 4636

Schwarzschild M., 1979, ApJ, 232, 236

Seth A. C. et al., 2014, Nature, 513, 398

Springel V., 2005, MNRAS, 364, 1105

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726

Springel V. et al., 2018, MNRAS, 475, 676

Taylor P., Kobayashi C., 2017, MNRAS, 471, 3856

Taylor E. N. et al., 2015, MNRAS, 446, 2144

Teklu A. F., Remus R.-S., Dolag K., Beck A. M., Burkert A., Schmidt A. S., Schulze F., Steinborn L. K., 2015, ApJ, 812, 29

Thater S. et al., 2017, A&A, 597, A18

Thater S., Krajnović D., Cappellari M., Davis T. A., de Zeeuw P. T., McDermid R. M., Sarzi M., 2019, A&A, 625, A62

Thater S. et al., 2022a, MNRAS, 509, 5416

Thater S. et al., 2022b, A&A, 667, A51

Thater S., Jethwa P., Lilley E. J., Zocchi A., Santucci G., van de Ven G., 2023a, preprint (arXiv:2305.09344)

Thater S., Lyubenova M., Fahriion K., Martín-Navarro I., Jethwa P., Nguyen D. D., van de Ven G., 2023b, A&A, 675, A18

Thomas J., Saglia R. P., Bender R., Thomas D., Gebhardt K., Magorrian J., Corsini E. M., Wegner G., 2007, MNRAS, 382, 657

Valenzuela L. M., Remus R.-S., Dolag K., Seidel B. A., 2024, preprint (arXiv:2404.01368)

Valluri M., Merritt D., Emsellem E., 2004, ApJ, 602, 66

van de Sande J. et al., 2017, ApJ, 835, 104

van de Sande J. et al., 2018, Nat. Astron., 2, 483

van de Sande J. et al., 2021, MNRAS, 508, 2307
## APPENDIX A: INSTANTANEOUS AND TIME-AVERAGED ORBITAL DISTRIBUTIONS

In order to justify the time-averaging of the orbital distributions within the simulation for comparison with our model-derived results from DYNAMITE, we create an *N*-body model of the EAGLE galaxy 9498015 using the Schwarzschild-derived best-fitting potential. Following the approach of Zhu et al. (2022), we keep the potential fixed (similar to what is being done using DYNAMITE) and allow the stellar particles to evolve within it. We then measure the instantaneous orbital distributions of the stellar particles at different times, as well as the time-averaged distribution. The time-average distribution is calculated from time-averaged positions and velocities of the stellar orbits.

In Fig. A1, we show the comparison between the instantaneous orbital distribution from the *N*-body model after 10 Gyr (left panel), the time-averaged distribution of the particles in the *N*-body model (middle panel), and the orbital distribution from the best-fitting Schwarzschild dynamical model (right panel). We also indicate the values of the orbital fractions within $1 R_e$ for each component.

The instantaneous orbital distribution for the *N*-body model is, as expected, very close to the true orbital distribution (shown in Fig. 9), and the time-averaged distribution aligns closely with what we obtain from the Schwarzschild models, with very similar orbital fractions within $1 R_e$. This suggests that time averaging is necessary to compare orbital fractions derived from simulations with the Schwarzschild model-derived values in a consistent manner.

van de Ven G., van den Bosch R. C. E., Verolme E. K., de Zeeuw P. T., 2006, A&A, 445, 513

van de Ven G., de Zeeuw P. T., van den Bosch R. C. E., 2008, MNRAS, 385, 614

van den Bosch R. C. E., van de Ven G., Verolme E. K., Cappellari M., de Zeeuw P. T., 2008, MNRAS, 385, 647

van der Marel R. P., 1993, in Danziger John I., Zeilinger W. W., Kjaer K., eds, The Velocity Dispersion Anisotropy of Elliptical Galaxies European Southern Observatory Conf. Workshop Proceedings.European Southern Observatory (ESO), Garching, p. 79

van der Marel R. P., Cretton N., de Zeeuw P. T., Rix H.-W., 1998, ApJ, 493, 613

Vasiliev E., Athanassoula E., 2015, MNRAS, 450, 2842

Vasiliev E., Valluri M., 2020, ApJ, 889, 39

Veale M., Ma C.-P., Greene J. E., Thomas J., Blakeslee J. P., McConnell N., Walsh J. L., Ito J., 2017, MNRAS, 471, 1428

Verolme E. K. et al., 2002, MNRAS, 335, 517

Virtanen P. et al., 2020, Nat. Method, 17, 261

Wang B., Cappellari M., Peng Y., Graham M., 2020, MNRAS, 495, 1958

Wiersma R. P. C., Schaye J., Smith B. D., 2009a, MNRAS, 393, 99

Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009b, MNRAS, 399, 574

Wilkinson M. J., Ludlow A. D., Lagos C. d. P., Fall S. M., Schaye J., Obreschkow D., 2023, MNRAS, 519, 5942

Zemp M., Gnedin O. Y., Gnedin N. Y., Kravtsov A. V., 2011, ApJS, 197, 30,

Zhu L. et al., 2018a, Nat. Astron., 2, 233

Zhu L. et al., 2018b, MNRAS, 473, 3000

Zhu L., van de Ven G., Méndez-Abreu J., Obreja A., 2018c, MNRAS, 479, 945

Zhu L. et al., 2020, MNRAS, 496, 1579

Zhu L. et al., 2022, A&A, 660, A20

<figure><img src="image_20.png" /><figcaption>**Figure A1.** Example galaxy EAGLE 9498015: orbit circularity $\lambda_z$ versus radius (in kpc). The colour indicates the density, *w*, of the orbits on the phase space, and the three horizontal black dashed lines indicate $\lambda_z$ = 0.8, 0.25, and -0.25, dividing the orbits into four regions (cold, warm, hot, and counter-rotating orbits). Fractions of mass within $1 R_e$ in each component are annotated in the figures. The vertical grey dashed line is located at 1 $R_e$. Left panel shows the instantaneous orbital distribution of the stellar particles in the *N*-body model after 10 Gyr; middle panel show the time-averaged orbital distribution of the stellar particles in the *N*-body model; and right panel shows the time-averaged orbital distribution derived from the dynamical modelling. Time-averaging helps to reconstruct the orbital distribution in a similar manner to what is derived by DYNAMITE.</figcaption></figure>
## APPENDIX B: FRACTION OF ORBITS FROM INSTANTANEOUS DISTRIBUTIONS

We show in Fig. B1 the one-to-one comparison of the instantaneous and model-recovered fractions of the cold, warm, and hot orbits. We find that, for $\lambda_z$, the comparison between the model-derived and true fractions of cold and counter-rotating orbits are generally in

agreement, within the errors. However, model-derived fractions of hot orbits are higher than the instantaneous fractions (with a median offset of 0.28), while the model-derived fractions of warm orbits tend to be lower than the instantaneous ones (with a median offset of 0.23). As mentioned in Section 2.4, this is not surprising since the model-derived orbital distributions are time-averaged.

<figure><img src="image_21.png" /><figcaption>**Figure B1.** One-to-one comparison of the true and model-recovered fraction of the cold (in blue), warm (in orange), hot (in red), and counter-rotating (in green) orbits. The true values are calculated from the instantaneous orbital distribution. We find that, for $\lambda_z$, the comparison between the model-derived and true fractions of cold and counter-rotating orbits are in general agreement, within the errors. However, model-derived fractions of hot orbits are generally higher than the instantaneous fractions, while the model-derived fractions of warm orbits tend to be lower than the instantaneous ones.</figcaption></figure>
¹International Centre for Radio Astronomy Research (ICRAR), M468, Uni-
versity of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia

²ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia

³Research Centre for Astronomy, Astrophysics, Astrophotonics, School of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109, Australia

⁴Astrophysics Subdepartment, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK

⁵Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17, A-1180 Wien, Austria

⁶Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia

⁷Sydney Institute for Astronomy, School of Physics, University of Sydney, NSW 2006, Australia

⁸Instituto de Astrofísica de Canarias, c/ Vía Láctea s/n, E-38205-La Laguna, Tenerife, Spain

⁹Departamento de Astrofísica, Universidad de La Laguna, E-38205-La Laguna, Tenerife, Spain

¹⁰School of Physics, University of New South Wales, Sydney, NSW 2052, Australia

¹¹Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, F-67000 Strasbourg, France

¹²Institute for Advanced Study, University of Strasbourg, 5 allée du Général Rouvillois, F-67083 Strasbourg, France

¹³Department of Physics and Astronomy, University of the Western Cape, Cape Town 7535, South Africa

¹⁴SISSA International School for Advanced Studies, Via Bonomea 265, I-34136 Trieste, Italy

¹⁵School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

¹⁶Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 München, Germany

This paper has been typeset from a T\_EX/L\_T\_EX file prepared by the author.