The MAGPI survey: evidence against the bulge–halo conspiracy Open Access

ABSTRACT

Studies of the internal mass structure of galaxies have observed a ‘conspiracy’ between the dark matter and stellar components, with total (stars|$+$|dark) density profiles showing remarkable regularity and low intrinsic scatter across various samples of galaxies at different redshifts. Such homogeneity suggests the dark and stellar components must somehow compensate for each other in order to produce such regular mass structures. We test the conspiracy using a sample of 22 galaxies from the ‘Middle Ages Galaxy Properties with Integral field spectroscopy’ Survey that targets massive galaxies at |$z \sim 0.3$|⁠. We use resolved, 2D stellar kinematics with the Schwarzschild orbit-based modelling technique to recover intrinsic mass structures, shapes, and dark matter fractions. This work is the first implementation of the Schwarzschild modelling method on a sample of galaxies at a cosmologically significant redshift. We find that the variability of structure for combined mass (baryonic and dark) density profiles is greater than that of the stellar components alone. Furthermore, we find no significant correlation between enclosed dark matter fractions at the half-light radius and the stellar mass density structure. Rather, the total density profile slope, |$\gamma _{\mathrm{tot}}$|⁠, strongly correlates with the dark matter fraction within the half-light radius, as |$\gamma _{\mathrm{tot}} = (1.3 \pm 0.2) f_{\mathrm{DM}} - (2.44 \pm 0.04)$|⁠. Our results refute the bulge–halo conspiracy and suggest that stochastic processes dominate in the assembly of structure for massive galaxies.

1 INTRODUCTION

Since the discovery that the observed kinematics of galaxies do not follow the gravitational potential established by the stellar mass alone, dark matter has been a ubiquitous ingredient in the field of galaxy evolution, and is thought to dominate the Universe’s mass budget (Blumenthal et al. 1984). The theorized assembly history of galaxies has baryonic matter fall into a potential well seeded by dark matter in the early Universe (White & Rees 1978; Wechsler & Tinker 2018). As such, the inclusion of dark matter in the gravitational potential of galaxy dynamical models is now standard. Of particular interest is how dark matter interfaces with the baryonic matter of galaxies, for example the scale at which it becomes the dominant mass component, and what impact it has on the global potential and mass distribution within a galaxy.

Efforts to measure the combined baryonic and dark matter content of galaxies have often assumed the distribution for the total mass density profile, where the combination of the baryonic and dark matter profiles is approximated as a simple power law, of form |$\rho \propto r^{\gamma }$|⁠, with |$\rho$| the mass density, r the galactocentric radius, and |$\gamma$| the slope of the profile (here we adopt the convention of |$\gamma \lt 0$|⁠). Measurements of this slope using dynamical studies (Cappellari et al. 2015; Poci, Cappellari & McDermid 2017; Bellstedt et al. 2018; Li et al. 2019; Veršič et al. 2024), lensing (Etherington et al. 2023), and combinations of the two (Auger et al. 2010; Barnabè et al. 2011; Ruff et al. 2011; Bolton et al. 2012; Sonnenfeld et al. 2013), have found galaxies tend towards mass density profiles with near-isothermal (⁠|$\gamma \sim -2$|⁠) slopes and limited intrinsic scatter, despite a significant range of stellar luminosity profiles. This observation has been dubbed the ‘bulge–halo conspiracy’, as the stellar component and dark matter component must jointly ‘conspire’ to yield isothermal total profiles with low scatter across (mainly early-type) galaxy populations (Dutton & Treu 2014). The conspiracy indicates remarkable homogeneity of mass structures for galaxies from the local Universe to |$z\sim 1$|⁠.

The physical origins of such a conspiracy are unclear. Although both components (baryonic and dark) must exist within the same gravitational potential, the stellar mass structure of galaxies is highly variable across galaxy types. This variability is reflected by a non-universal initial mass functions (IMF) both within (McConnell, Lu & Mann 2016; Conroy, van Dokkum & Villaume 2017; van Dokkum et al. 2017; Martín-Navarro et al. 2019) and between galaxies (Cappellari et al. 2012; Smith 2020; Martín-Navarro et al. 2021; Thater et al. 2023b), and measurements of galaxy stellar structure such as the Sérsic index, which show a large range in values (Graham 2013; Lange et al. 2015; Zahid & Geller 2017). The combination of the stellar component and an independent dark matter component should, naively, increase the observed variance of a population of total density profiles, instead of the small scatter and near isothermal profiles observed.

The conspiracy was explicitly tested by measuring de-coupled stellar and dark matter profiles of |$\mathrm{ATLAS}^{\mathrm{3D}}$| galaxies with Jeans (Cappellari 2008; Cappellari 2020) dynamical modelling (Poci et al. 2017). Poci et al. (2017) found that the intrinsic scatter in the slope of the combined density profiles was only marginally larger than that of the stellar profiles alone, and could therefore make no strong assessments about the validity of the ‘conspiracy’. Their approach included two distinct models; one in which the gravitational potential was treated a single ‘total-mass’ density profile, and another in which the stars and dark matter were modelled explicitly. Jeans modelling method used assumes galaxies have an axisymmetric potential, and Poci et al. (2017) further assumed a cylindrically aligned velocity ellipsoid and a global orbital anisotropy term, which might not reflect the true orbital structures present in galaxies. Crucially, the majority of the |$\mathrm{ATLAS}^{\mathrm{3D}}$| sample has kinematic data only interior to the half-light radius scale, constraining the stellar and dark components on a small physical region. Within the (presumably baryon dominated) half-light radius, it may be expected that the total and stellar mass structures closely trace each other.

In this work, we leverage ‘Middle Ages Galaxy Properties with Integral field spectroscopy’ (MAGPI) Survey¹ data to construct Schwarzschild orbit-based models (Schwarzschild 1979) of massive galaxies at |$z \sim 0.3$|⁠. The benefit of studying this epoch is that the data frequently extends to two or even three half-light radii with a resolution that is comparable to local surveys such as the Sydney Australian Astronomical Observatory Multi-object Integral Field Spectrograph (SAMI) Survey (Croom et al. 2012; Santucci et al. 2022). Furthermore, it is valuable to investigate the dynamical state of galaxies during this period of the Universe’s history. We present triaxial Schwarzschild dynamical models of 22 MAGPI galaxies, constructed using the dynamite code (Jethwa et al. 2020; Thater et al. 2022b), which relaxes some of the assumptions necessary for Jeans dynamical modelling. As an orbit-based technique, Schwarzschild models that accurately reconstruct galaxy observables (flux and higher order stellar kinematics) provide insight into the intrinsic properties of galaxies, such as 3D shape, orbital structure, and dark matter fractions. We use these detailed models to investigate how the internal mass of a galaxy is distributed between its stellar and dark components, and test the predictions of the bulge–halo conspiracy.

In Section 2 we introduce the MAGPI data, and present how the inputs to the Schwarzschild models are derived in Section 3. The construction of the Schwarzschild models is given in Section 4, with our results from these models presented in Section 5. We discuss our results in Section 6 and summarize this work in Section 7.

Throughout this article, we adopt a flat |$\Lambda$|CDM cosmology with |$H_0 = 67.1\ \mathrm{km\, s}^{-1} \mathrm{Mpc}^{-1}$| and |$\Omega _m = 0.3175$|⁠, corresponding to the Planck 2013 cosmology (Planck Collaboration XVI 2014). All scales are converted (angular and physical units) using the angular diameter distance given by the object’s redshift and the above cosmology.

2 DATA SAMPLE

In this work we use the 30 objects for which Jeans dynamical models were constructed by Derkenne et al. (2023), although due to sample cuts described in Section 4.3, we present results for only 22 of these objects. Modelling the same sample with two different techniques provides a valuable comparison point, presented in Section 5.4. Importantly, in Derkenne et al. (2023) we showed that the spatial resolution of this sample is high enough that the recovered mass distribution of the models are not biased due to point-spread function (PSF) effects, which is important to consider when constructing dynamical models of objects at any appreciable redshift. The data are deep enough to provide up to |$\sim 4$| half-light radii of stellar kinematic coverage in some cases.

The sample is drawn from the MAGPI Survey, a VLT/MUSE Large Programme which targets massive galaxies in the Universe’s ‘middle ages’. The survey uses MUSE in wide-field mode, yielding 60 arcsecond on a side fields. Resulting spectra span the range |$4650-9300\,$|Å in steps of |$1.25\,$|Å per pixel, with ground-layer adaptive optics typically resulting in PSF full-width half-maxima (FWHM) of |$\sim 0.6$| arcsec, equivalent to 4.629 kpc at |$z = 0.3$|⁠. Fields are observed for a total on-source integration time of 4.4 h per field. The detailed aims of the MAGPI Survey are presented in Foster et al. (2021). Although the survey officially targets 60 massive primary galaxies, the fields obtained to date also contain many secondary objects that have adequate angular resolution for resolved studies.

The full details of the data reduction process will be presented by Mendel et. al. (in preparation), but we provide an overview here. Data reduction is based on the MUSE reduction pipeline (Weilbacher et al. 2020). The Zurich atmosphere purge (zap) software is used for sky subtraction (Soto et al. 2016). Segmentation maps for all sources are made using profound (Robotham et al. 2018). A ‘minicube’ is constructed for each source identified from the segmentation maps for every field, with the identified source at the centre and minicube extent such that it encompasses the full dilated segment. Images for each field are made from the collapsed spectral cube of each field.

Integrated stellar masses for the sample are calculated using the spectral energy distribution fitting software prospect (Robotham et al. 2020), assuming a Chabrier (2003) initial mass function. The spectral energy distribution is fit using pixel-matched imaging in the ugriZYJHKs bands from the GAMA Survey (Bellstedt et al. 2020), with photometry derived by the MUSE-based profound segmentation maps for each galaxy. The sample spans stellar masses between |$\log _{10}(M_{\star }/\mathrm{M}_\odot) = 10.4$| and |$\log _{10}(M_{\star }/\mathrm{M}_\odot) = 11.6$|⁠, with a sample mean redshift of |$z = 0.304$|⁠. From team visual morphology classifications, our sample is dominated by lenticular and early-type morphologies, with only 4/22 objects showing evidence of spiral arms (Foster et. al., in preparation).

3 INPUTS TO SCHWARZSCHILD MODELS

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 surface brightness. In the following sections we describe our modelling of the galaxy surface brightness and our extraction of the 2D stellar kinematics. Additionally, for our modelling approach, we transform the surface brightness models to stellar mass models via spectral measurements of stellar mass-to-light ratios.

3.1 Galaxy surface brightness models

We parametrize the galaxy light as a series of N concentric 2D Gaussians (Emsellem, Monnet & Bacon 1994) using the Python package mgefit (Cappellari 2002). This multi-Gaussian expansion (MGE) is a projected model of the surface brightness, which is used to describe the intrinsic (deprojected) luminosity density of the stellar tracer within a given gravitational potential. The fitting process is identical to that used by Derkenne et al. (2023), where we presented detailed explorations of how pixel sampling and resolution effects impact MGE models, benchmarked against Hubble Space Telescope results. We use a constant position angle in the MGE fits, as a radially varying position angle was not warranted by the fit quality of the models compared to the data. We choose the r-band images built from the MUSE data cubes for their depth, with the initial image size selected to have sides of 20 |$R_{\mathrm{e}}$| from the profound estimate. The choice of (observed-frame) r-band ensures that the observed frame light used to create the surface brightness models is similar to the spectral window used to measure the galaxy kinematics. Due to the redshift of the targets, the observed light in the r-band corresponds approximately to rest-frame g-band of the sources. Therefore, in our analysis, we calibrate the projected surface luminosity density of each source in the (rest-frame) g-band (equation 4). All neighbouring sources are masked out using the segmentation maps. MGEs are also made of stacked point sources in each field to model the PSF in order to analytically convolve the surface brightness model with the PSF before comparing to observations. In this way, the best-fitting model provides a deconvolved description of the surface brightness, corrected for the effects of PSF.

We require the galaxy surface brightness models to be in units of |$\mathrm{L}_\odot \,\mathrm{pc}^{-2}$| (and eventually the mass models in |$\mathrm{M}_\odot \,\mathrm{pc}^{-2}$|⁠). First, the enclosed counts (C) of each Gaussian component are re-normalized to peak surface brightness (in counts per pixel), |$C_0$|⁠:

with q the axial ratio of the Gaussian and |$\sigma$| the dispersion. The counts are converted to ST magnitudes (⁠|$m_{\mathrm{ST}}$|⁠) via

The g-band peak surface brightness in magnitudes per square arcsecond is given by

where s is the pixel scale, |$A_r$| is the line-of-sight dust correction from Schlafly & Finkbeiner (2011), and k is the k-correction which also converts from r to g band. The k-correction is necessary as the MAGPI sources are sampled at various rest-frame frequencies due to their redshift, which makes comparing across sources without a correction invalid. We follow Hogg et al. (2002) equation (8) to perform the correction, and in practice calculate the difference in magnitudes between the observed frame spectrum and its rest-frame counterpart in the r and g bands, respectively. The final term in equation (3) is the cosmological surface brightness dimming correction, following the derivation given by Whitney et al. (2020) for ST magnitudes. We convert the g-band surface brightness to a projected 2D surface density, |$I^{\prime }$|⁠, in units of |$L_{\odot }\,\mathrm{pc}^{-2}$|⁠, with:

where |$M_g$| is the absolute magnitude of the sun in the g-band from Willmer (2018). Finally, the dispersion of each Gaussian component is transformed from units of pixels to arcseconds by multiplying by the MUSE pixel scale. We show an example of an MGE fit and its residuals in Fig. 1. In the subsequent analysis, we define the half-light radius, |$R_{\mathrm{e}}$|⁠, as the semimajor axis of the best-fitting MGE model ellipse which contains half the model light.

3.2 Stellar kinematics

We extract 2D stellar kinematics using the Python package ppxf (Cappellari & Emsellem 2004; Cappellari 2017), fitting for four kinematic moments; the line-of-sight velocity, velocity dispersion, 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 & Franx 1993). To prepare the minicubes for spectral fitting, we remove all spaxels with a signal-to-noise of less than 1.5 per pixel in the fitted spectral region before Voronoi binning to a signal-to-noise of 10 per bin with the Python package vorbin (Cappellari & Copin 2003). The exact fitted region changes on a galaxy-by-galaxy basis due to the variation of the sample around the |$z=0.3$| target, but is approximately rest-frame g-band.

We use the Indo-US template library (Valdes et al. 2004) due to its high spectral resolution when attempting to measure accurate velocity dispersions with the moderate redshift MAGPI data. The template library was convolved to match the wavelength-dependent line-spread function of the MUSE/MAGPI data, assumed to be spatially invariant in each field. We select an optimal subset of Indo-US templates with an initial aperture spectrum fit for efficiency. We restrict the individual Voronoi bin fits to the top 20 weighted templates from the aperture fit. The aperture spectrum is constructed by co-adding spaxels within the half-light ellipse from the MGE model. For all fits we use fifth degree multiplicative polynomials.

We estimate uncertainties on the fitted kinematic moments by shuffling the ppxf residuals and re-distributing them on the input spectrum before re-fitting across 100 iterations and taking the standard deviation. Since this process adds noise to a data spectrum that already contains noise, we reduce the derived parameter uncertainties by |$\sqrt{2}$| to account for this noise doubling. We note this process is the same as the kinematics presented in Derkenne et al. (2023), and for the same objects, except here we fit for four kinematic moments instead of two, and the template broadening here is wavelength dependent.

As the MGEs and stellar kinematics are derived from the same underlying MUSE data, we rotate the kinematic x and y coordinates by the photometric position angle determined from the MGE fits. This rigorously ensures that the gravitational potential described in Section 4.1 and kinematics are correctly aligned within the Schwarzschild models.

3.3 Stellar mass-to-light ratios

Using the same signal-to-noise cuts and Voronoi bins as in Section 3.2, we extract stellar population measurements using the E-MILES stellar populations models of Vazdekis et al. (2016). These templates allow for the calculation of a spectral mass-to-light ratio whilst also accounting for variations of stellar age and metallicity. In practice, we use the included ppxf E-MILES utilities for the fit, following Cappellari et al. (2013b) equation (2). The mass-to-light ratios are calculated with a Salpeter (Salpeter 1955) IMF, however we note we leave the global scaling of mass in our Schwarzschild models free, so that resulting models are not dependent on any particular IMF. This process assumes the IMF is spatially constant within a given galaxy.

To convert the luminous MGE model into a stellar mass model, an ellipse is constructed for each individual MGE component, using its respective |$\sigma$| and q values. This ellipse represents the region where that MGE component’s contribution is maximal. The mean g-band mass-to-light ratio along this ellipse is associated with that component. All ellipses that lie beyond the coverage of the mass-to-light map are set to have the same mass-to-light ratio as the last ellipse contained within the map. That is, we assume that at large radii the mass-to-light ratio becomes constant, which avoids divergent solutions, and as seen, for example, in the detailed mass-to-light profile of NGC 5120 (Mitzkus, Cappellari & Walcher 2017), and for FCC 47 (Thater et al. 2023b). Each of the N Gaussian surface density components in units of |$L_{\odot }\,\mathrm{pc}^{-2}$| are multiplied by their corresponding mass-to-light ratio, so that the final 2D projected galaxy surface density model is in units of |$\mathrm{M}_\odot \,\mathrm{pc}^{-2}$|⁠. This method for converting surface brightness models to stellar mass follows Poci et al. (2017). We show an example stellar mass-to-light ratio map in Fig. 2.

We do not include uncertainties on the stellar MGE or the stellar population mass-to-light ratios, as perturbing the gravitational potential for error analysis with Schwarzschild models is prohibitively expensive in terms of computational time and energy. While the stellar mass-to-light ratio maps tend to be noisy, they do exhibit structure, and our technique of averaging the mass-to-light ratio along each of the Gaussian components helps reduce the impact of noise. The example we show in Fig. 2 is neither the worst nor best from the sample.

4 SCHWARZSCHILD ORBIT-BASED MODELS

Schwarzschild models are a flexible but computationally expensive method by which to create exceptionally detailed dynamical models of galaxies, recovering intrinsic orbital properties from projected radial velocity and surface density information (Schwarzschild 1979). We use the triaxial dynamite code (Jethwa et al. 2020) to construct the Schwarzschild models presented in this work, following the theoretical outline of van den Bosch et al. (2008) with orbit mirroring corrections included (Quenneville, Liepold & Ma 2022; Thater et al. 2022b). The models are computationally expensive as a large library of orbits that are allowable within 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 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 optimization. In practice, this can require running thousands of models to characterize a parameter space, and large orbit libraries that do not artificially restrict the various orbital structures a galaxy may have. We outline our construction of the gravitational potential, orbit libraries, and parameter space searches in the following sections.

4.1 Constructing the gravitational potential

The 2D axisymmetric MGE model in units of |$\mathrm{M}_\odot \,\mathrm{pc}^{-2}$| can be converted to a triaxial 3D mass density in units of |$\mathrm{M}_\odot \,\mathrm{pc}^{-3}$| using the following expression:

where |$p_i$| = |$B_i/A_i$|⁠, |$q_i$| = |$C_i/A_i$|⁠, and |$A_i$|⁠, |$B_i$|⁠, and |$C_i$| are the major, intermediate, and minor axes of the 3D triaxial Gaussian. The corresponding projected quantities (⁠|$p^{\prime }$|⁠, |$q^{\prime }$|⁠) can be transformed into the intrinsic quantities using the viewing angles of the object (Cappellari 2002; van den Bosch et al. 2008). |$M_i$| represents the mass of the i-th Gaussian component with dispersion |$\sigma _i$|⁠.

Following Zhu et al. (2018), for the dark matter halo we adopt a spherical Navarro–Frenk–White (NFW) halo (Navarro, Frenk & White 1997), where the enclosed mass profile is given as

with r the radius, c the concentration of the halo, |$g(c) = [\ln (1+c)-c/(1+c)]^{-1}$|⁠, and |$M_{200}$| is the virial mass within a radius of |$r_{200}$|⁠. The dark matter halo is then fully described by two free parameters: the concentration and the virial mass. However, we found, as did Zhu et al. (2018) with local Universe data, that for MAGPI data quality the two dark matter halo parameters were unconstrained in our models. Instead we adopted an |$M_{200}-c$| coupled halo. This fixes the concentration according to the relation measured by Dutton & Macciò (2014),

with |$h = 0.671$|⁠. The parameters a and b evolve with redshift as

and

given as equations (10) and (11) in Dutton & Macciò (2014). We fix the |$M_{200}-c$| relation at the mean redshift of the MAGPI sample, |$z = 0.304$|⁠. The critical density at this redshift with the adopted cosmology is |$\rho _c = 1.73 \times 10 ^{-7}\, \mathrm{M}_\odot \mathrm{pc}^{-3}$|⁠. Our models therefore have only one free parameter to describe the dark matter halo, |$M_{200}/M_{\star }$|⁠, where |$M_{\star }$| is the total stellar mass within the virial radius. For completeness, although well below the resolution of MAGPI data, we include a black hole with black hole mass given by the redshift-dependent mass–velocity dispersion relation from Robertson et al. (2006). The total mass (luminous and dark) of the system can then be scaled by the global constant |$\Upsilon _r$|⁠, as the scaled model can be simply evaluated by scaling the orbital velocities without re-integrating the orbit library.

4.2 Building the orbit library

Only certain regular orbital families can exist within any given gravitational potential. Each unique model therefore has its own library of orbits which must be numerically integrated for a sufficient number of periods so as to be representative of the time-averaged orbital properties. The orbits themselves are characterized by the three conserved integrals of motion: binding energy E, z-component angular momentum |$L_z$|⁠, and the third, non-classical integral of motion, |$I_3$| (van der Marel et al. 1998; Cretton et al. 1999).

We construct the orbit library by spanning 21 intervals in E, 15 intervals in |$L_z$|⁠, and 11 intervals in |$I_3$|⁠. We choose this relatively large orbit library size based on the exploration of orbit library size conducted by Quenneville et al. (2022). They found that with |$L_z$| sampled below 15 intervals there were spurious oscillations in the |$\chi ^2$| space of the model, yielding unstable solutions. However, increasing the orbit library to larger sizes does not continuously improve the models, as successively fewer orbits are given significant weighting (Jin et al. 2019). Our orbit library size is a compromise between large enough not to restrict possible orbital structures, and small enough to be computationally efficient. The dynamite code constructs three sets of orbit libraries: box orbits, tube-orbits, and counter-rotating tube-orbits, such that the total orbit library size is |$3 \times 21 \times 15 \times 11$|⁠. We do not regularize the orbits by use of dithering, as dithering vastly increases the computational time required for each model. Dithering seeds |$N^3$| orbits in a mini-grid around each set of initial conditions, where N is the chosen dithering value. We explore the impact of dithering in Appendix  A, and conclude it is not necessary for MAGPI-like data quality.

The minimum and maximum orbital radius is set as |$0.5\sigma _{\mathrm{min}}$| and |$5\sigma _{\mathrm{max}}$|⁠, where |$\sigma _{\mathrm{min}}$| and |$\sigma _{\mathrm{max}}$| are the minimum and maximum Gaussian dispersions from the MGE fit, respectively, following Zhu et al. (2018). These radii set the minimum and maximum orbital energy within a given gravitational potential. We integrate orbits across 200 orbital periods, with 50 000 samples per orbit in the meriodional plane. Again following Zhu et al. (2018), we implement the ‘CRcut’ keyword in dynamite, which avoids the spurious up-weighting of counter-rotating orbits for galaxy regions with high |$v/\sigma$|⁠.

The observable properties of each orbit, stored on an intrinsic coordinate grid, are transformed to projected coordinates. The orbital properties (surface brightness and higher order stellar kinematic moments) are then randomly perturbed by the PSF before being assigned to the observational apertures (the Voronoi bins). In this way, the intrinsic models are convolved with the same PSF as the observational data (van den Bosch et al. 2008).

4.3 Parameter space search and model selection

To summarize the above, our implementation of the Schwarzschild models include five free parameters:

• The axial ratio |$p = B/A$|⁠, where B and A represent the medium and major axes of the 3D triaxial Gaussian, respectively.

• The axial ratio |$q = C/A$|⁠, where C is the minor axis of the 3D triaxial Gaussian.

• The axial ratio |$u = A^{\prime }/A$|⁠, where |$A^{\prime }$| is the projected major axis of the 3D triaxial Gaussian.

• The global mass scaling factor, which scales both stellar and dark components jointly, |$\Upsilon _r$|⁠.

• The dark matter content, parametrized as |$M_{200}/M_{\star }$|⁠, where |$M_{200}$| is the total mass within the virial radius |$r_{200}$|⁠, and |$M_{\star }$| is the total stellar mass.

The values of p, q, and u are all bounded such that a valid deprojection of the 2D stellar mass model exists. These conditions are that |$0 \lt q \le p \le 1$|⁠; |$\mathrm{max}(q/q_{\mathrm{obs}}, p) \lt u$| with |$q_{\mathrm{obs}}$| the minimum axial ratio from the stellar surface brightness MGE model; and |$u \lt \mathrm{min}(p/q_{\mathrm{obs}}, 1)$|⁠. Due to these conditions, the exact bounds used in the models vary for each galaxy.

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 _r$| between 0.1 and 1.5, and at least seven values in dark matter fraction |$M_{200}/M_{\star }$| between -2 and 3 in log space.

Each resulting model is compared to the kinematic observations to determine the best-fitting model via the chi-square value, calculated from the observed kinematics and the Schwarzschild model projected kinematics. The galaxy surface brightness models are used to constrain the linear combination of orbits, but are not used to evaluate the quality of fit between different models, which is done using only the kinematic moments. 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 dark matter fraction, 0.02 for p and q, 0.01 for u (which generally tends very close to unity), and 0.01 in |$\Upsilon _r$|⁠. Practically, the coarse and fine model runs use the ‘FullGrid’ and ‘LegacyGridSearch’ options in the dynamite code, respectively. We allow the iterations to continue for up to 6000 coarse models and 2000 fine models before terminating, although the code can terminate sooner once the |$\chi ^2$| of new models makes no further improvement, quantified as an absolute reduction in the |$\chi ^2$| value of less than 0.1.

As in Zhu et al. (2018), we define the |$1\sigma$| confidence level as |$\Delta \chi ^2 = \sqrt{2N_{\mathrm{kin}}N_{\mathrm{mom}}}$|⁠, where |$N_{\mathrm{kin}}$| is the number of kinematic Voronoi bins and |$N_{\mathrm{mom}}$| is the number of kinematic Gauss–Hermite moments. The coarse and fine modelling runs are combined to establish all models within the |$1\sigma$| confidence region. 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 taken. For example, the upper error is quoted as |$\mathrm{max}({\bf X})-x_0$|⁠, where |${\bf X}$| is the vector of all |$1\sigma$| model parameter values and |$x_0$| is the value from the lowest |$\chi ^2$| model. An example of the parameter space coloured by Schwarzschild model |$\chi ^2$| value is shown in Fig. 3.

Finally, we assess the quality of the resulting models. As noted above, the model iterations were informed by the kinematic |$\chi ^2$| value only, and so there is no requirement for the best-fitting model’s surface brightness distribution to match the observations. In some cases, the kinematics were well reproduced but not the surface brightness, which in turn means the galaxy mass distribution is not well represented by the model. We therefore select only models where the median absolute deviation of the surface brightness residuals (data-model/data) was less than 1 per cent, ensuring only galaxies with near structureless surface brightness residuals are included in the final sample. This strict cut leaves 22 out of a possible 30 objects. We note the cut galaxies are generally spiral systems with bars or kinematic twists, leaving a sample that is dominated by early-type systems. We show two example Schwarzschild models in Figs 4 and 5, which are representative of the selected sample. We show MAGPI object 1 525 170 222 in particular because it has one of the highest reduced |$\chi ^2$| values of the sample at 2.43 (i.e. one of the worst fits) but a model that well represents the observed kinematics, to demonstrate the general fit quality of the models. The remaining 20 Schwarzschild models are all shown in Appendix  B.

5 RESULTS

From our Schwarzschild models we can extract detailed, intrinsic orbital structures of galaxies. In the following sections we present our analysis of the 22 MAGPI galaxies with well-fitting Schwarzschild models. We first examine the bulge–halo conspiracy with our data, showing in particular the large scatter of total density profile slopes and how the dark matter fractions are independent from the stellar mass structure. We then present an analysis on the intrinsic shapes and orbital structures of the sample. Finally we present a comparison of our total mass density profile slopes from the Jeans models presented in Derkenne et al. (2023) and Schwarzschild models derived in this work.

5.1 Testing the bulge–halo conspiracy

We examine two hypotheses of the bulge–halo conspiracy. The first is that populations of galaxies have low intrinsic scatter in their total mass density profile structure, quantified through the total mass density profile slope. From this, we would expect the observed scatter in total mass density slopes to be less than or equivalent to the observed scatter in the stellar mass density profile slopes. Secondly, if the dark and baryonic components conspire through some as of yet unknown mechanism, then the dark matter content of galaxies should correlate with the stellar density structure.

The structure and gradient of the stellar mass profile is derived using the stellar mass models described in Section 3.3 (i.e. the mass-to-light adjusted MGE models), scaled by the best-fitting parameter |$\Upsilon _r$| from the Schwarzschild models. The dark matter profile is determined by the fitted dark halo in the Schwarzschild models, with the total profile being the simple addition of the stellar and dark components. The total, stellar, and dark matter mass density profiles for the whole sample are shown in Fig. 6.

Following the approach of previous works (e.g. Cappellari et al. 2015), we measure the stellar and total mass density slopes – referred to as |$\gamma _{\star }$| and |$\gamma _{\mathrm{tot}}$|⁠, respectively – assuming a single power law of the form:

where |$\rho$| is the mass density and r is the galactocentric radius. As in Derkenne et al. (2023), we fit the slope for the radial region between |$R_{\mathrm{e}}$|/10 and 2|$R_{\mathrm{e}}$|⁠.

From the power-law fits to the MAGPI sample mass density profiles we examine the total profile and stellar profile slope distributions, with the results shown in Fig. 7. The stellar slopes are steeper on average than the total profile slopes, as the baryonic component dominates the central regions but not the outskirts, and the inclusion of a dark matter halo tends to cause the combined profile to become more shallow. We note that the scatter of the total density profiles is larger than for the stellar component alone, comparing |$\sigma _{\mathrm{tot}} = 0.30 \pm 0.03$| to |$\sigma _{\star } = 0.19 \pm 0.02$| (errors here are the standard error on the standard deviation). We find a super-isothermal median total mass density profile slope of |$\gamma _{\mathrm{tot}} = - 2.27 \pm 0.08$| (with the uncertainty characterized as the standard error on the median), and the sample spans |$\gamma _{\mathrm{tot}} \sim -1.6$| to |$\gamma _{\mathrm{tot}} \sim -2.7$| in slope values. For the stellar component, we find a median slope of |$\gamma _{\star } = -2.28 \pm 0.05$| (standard error of the median) and the sample spans between |$\gamma _{\star } \sim -2.2$| to |$\gamma _{\star } -2.85$|⁠. As might naively be expected, the spread of values in the total mass density profile slopes, which combines the stellar light structure with a free dark matter halo, is higher than the spread of the stellar mass profiles alone. By contrast, the bulge–halo conspiracy predicts reduced or equivalent scatter when considering the total mass density profile slopes against the stellar mass profile slopes.

The variability of the total mass density profile slopes presented in this work is in contrast to some previous works. Gravitational lensing studies generally report a scatter of |$\lesssim 0.2$|⁠, such as Barnabè et al. (2011), Li, Shu & Wang (2018), and Etherington et al. (2023). However, Ruff et al. (2011) find comparable scatter to our results with |$\sigma _{\gamma } = 0.25$|⁠. Several dynamics studies have also found comparatively low scatter, with Cappellari et al. (2015) finding |$\sigma _{\gamma } = 0.11$| and Bellstedt et al. (2018) finding |$\sigma _{\gamma } = 0.13$|⁠, but Li et al. (2019) finding a higher value of |$\sigma _{\gamma } = 0.22$|⁠. Poci et al. (2017) explicitly compared the slope distributions between the stellar mass profile slopes and total mass density slopes and found equivalent scatter for the |$\mathrm{ATLAS}^{\mathrm{3D}}$| sample, and could not refute the bulge-conspiracy on that basis. However, the |$\mathrm{ATLAS}^{\mathrm{3D}}$| data typically extends to only a half-light radius, which limits how well the dark matter halo can be constrained. Poci et al. (2017) found |$\sigma _{\gamma } \sim 0.17$| for both the stellar and total mass density profile slope distributions.

We can also test the correlations between the total profile slopes, the dark matter content within the half-light radius, and the stellar profile slopes. We show these results in Fig. 8. We note our results are unchanged if we consider the dark matter profile slopes instead of dark matter fractions. We find a moderate correlation between the stellar mass density profile slope and the total mass density profile slope, in that a galaxy with a steep stellar component (where the stellar mass density drops off quickly with radius) also has a steep total component, with a Pearson correlation coefficient of 0.55, implemented through scipy (Virtanen et al. 2020). This makes sense, as the inner regions of the total slope are generally just the stellar component. For objects with no dark matter, the total mass density slope must equal the stellar mass density slope by construction. The total density profile slopes are strongly correlated with the dark matter fractions within the half-light radius, with a correlation coefficient of 0.83, such that shallower total profiles have a higher dark matter fraction. However, we notice the lack of even a mild correlation between the half-light dark matter fraction and the stellar profile slope, with a correlation coefficient of 0 and p-value of 0.98. This result indicates that the stellar mass structure of a galaxy is independent of its dark matter content, again at odds with the supposition of a ‘conspiracy’ between the stellar and dark matter components.

From our results, we argue that the accumulation of the dark matter halo perturbs the total density slope from steep to shallow values, but that the stellar mass structure does not determine the structure of the dark halo (or the reverse). That the total slope correlates with dark matter fraction is a prediction from the Magneticum (Dolag 2015) simulations. Remus et al. (2017) found that the total slope and dark matter fraction co-evolves for all systems across all redshifts, with a simple relation of the form |$\gamma _{\mathrm{tot}} = 1f_{\mathrm{DM}} -2.52$| with |$f_{\mathrm{DM}}$| the dark matter fraction within the 3D half-mass radius. With the MAGPI data, we find a notably similar result, of |$\gamma _{\mathrm{tot}} = (1.3 \pm 0.2) f_{\mathrm{DM}}-(2.44 \pm 0.04)$|⁠, with |$f_{\mathrm{DM}}$| the dark matter fraction within the half-light radius. Intuitively, the relation between the two parameters makes sense. Galaxies with low dark matter content have total slopes that closely follow the stellar mass distribution. Subsequent (dry) mergers increase the dark matter fractions by increasing the effective radius of the galaxy with stars at large radii, causing the total slope to approach shallower (less negative) values as the extent of the galaxy encompasses more of the (shallower) dark matter profile.

The lack of correlation between the dark matter content and stellar mass structure is in tension with some previous observations. Tortora et al. (2018) investigated the relationship between central dark matter fractions and stellar density using aperture-based dynamical masses from Sloan Digital Sky Survey (SDSS; Abazajian et al. 2009) spectroscopy, and structural parameters from Kilo Degree Survey (KiDS; de Jong et al. 2013) photometry for a sample of |$\sim 3800$| galaxies spanning up to |$z \sim 0.65$|⁠, and found that there exists a strong anticorrelation between the dark matter content of a galaxy and its average central density. Galaxies with the steepest Sérsic indices, of |$n\sim 10$|⁠, were found to have the largest central dark matter fractions. Likewise, an earlier study by Grillo (2010) with a sample of |$1.7\times 10^{5}$| SDSS galaxies found a mild anticorrelation between central dark matter fractions and stellar surface densities, assuming an isothermal power-law density profile.

While it is uncertain what exactly causes the disagreement between those studies and this work, we note that in this work we examine stellar mass profiles instead of stellar light structure, and the Schwarzschild modelling technique implemented applied here to IFS data makes fewer assumptions than aperture based measurements of galaxy parameters. We observe no significant correlation between |$f_{\mathrm{DM}}$| and Sérsric indices with the MAGPI sample, using single-component indices calculated using galfit (Peng et al. 2002). We obtain a p-value of 0.34 and a correlation coefficient of 0.34, in keeping with our results between the lack of conspiracy between the stellar component and central dark matter fractions.

5.2 The dark matter content of massive galaxies

The Schwarzschild models of the MAGPI sample allow us to investigate in detail where each mass component, baryonic or dark, dominates. In Fig. 9 we show the half-light dark matter fractions as a function of stellar mass, and coloured by the radius, |$r_{\mathrm{cross}}$|⁠, where the dark matter density profile exceeds the stellar profile. As with the total density profiles themselves, we find that there is a spread in values for where dark matter becomes the dominant density component. The smallest crossing radius in the sample is MAGPI object 1501196198, where the dark matter density exceeds the stellar density within the half-light radius, and has a dark matter fraction within the half-light radius of |$66_{-8}^{+19}$| per cent. Numerous galaxies experience no transition within the data region and have negligible central dark matter fractions, but from the best-fitting Schwarzschild models have profiles that cross over in excess of 10–20 |$R_{\mathrm{e}}$|⁠, whereas other galaxies are dominated by dark matter in the 1–4 |$R_{\mathrm{e}}$| region. We conclude from this that assuming a common scale radius across a population of galaxies at which the dark matter dominates is problematic, and probably not reflective of realistic galaxy structures. This result is consistent with analysis performed by Harris et al. (2020b), who found remarkable variety in the dark matter content of a sample of 102 early-type galaxies up to 5 |$R_{\mathrm{e}}$| using X-ray gas to probe the density profiles. In their study, some systems are entirely dominated at 5 |$R_{\mathrm{e}}$| by dark matter (⁠|$f_{\mathrm{DM}} \gt 0.95$|⁠), and others have very little even at this large radius (⁠|$f_{\mathrm{DM}} \lt 0.3$|⁠).

For the MAGPI sample, the median dark matter fraction within the half-light radius is 10 per cent with a large sample standard deviation of 19 per cent. This is similar to local Universe measurements from |$\mathrm{ATLAS}^{\mathrm{3D}}$| and ‘Mapping Nearby Galaxies at Apache Point Observatory’ (MaNGA; Bundy et al. 2015) galaxies made using Jeans dynamical modelling, which depending on the dark matter halo adopted, sit between 8 and 17 per cent, with lower dark matter fractions reported for galaxies with higher data quality (Cappellari et al. 2013a; Zhu et al. 2024). Santucci et al. (2022) used Schwarzschild dynamical modelling on a sample of SAMI galaxies to determine their central dark matter fractions, the closest study in methodology to this work, finding a median of 28 per cent with a standard deviation of 20 per cent. However, in our models we have accounted for radially varying stellar population mass-to-light ratios, which might drive some of the observed differences in dark matter fractions. The fitted relations between stellar mass and dark matter fractions as found by Santucci et al. (2022) and Cappellari et al. (2013a) are shown in Fig. 9.

We note that in our analysis we have left the global mass scaling parameter, |$\Upsilon _r$|⁠, free to vary for each galaxy. This can effectively account for stellar IMF variations between galaxies (as well as any other residual mass normalization) as reported in the literature (Cappellari et al. 2012; Smith 2020); however it does not account for IMF variations within galaxies. Since the assumed IMF can adjust the stellar mass distribution, we explore the possible impact of a radially varying IMF in Appendix  C, which yield further uncertainties on our reported dark matter fractions of the order of 10 per cent. In future, and with data that can provide constraints on the IMF, a more complete treatment would include a fully fitted, 2D IMF as part of the modelling process.

5.3 Orbital structure and intrinsic shapes of massive galaxies

The weighting of each orbit in the gravitational potential of the best-fitting Schwarzschild model allows us to examine the orbital structure of galaxies, through the ‘circularity’ space as a function of radius. The circularity, |$\lambda _z$|⁠, is defined as

with |$\overline{L_z} = \overline{xv_y-yv_x}$|⁠, where v represents the instantaneous velocity in a coordinate direction (x, y, z) and |$\overline{L_z}$| is the time-average z-component of the angular momentum. Furthermore, |$r = \overline{\sqrt{x^2 + y^2 + z^2}}$|⁠, and |$\overline{V_c} = \sqrt{\overline{v_x^2 + v_y^2 + v_z^2+ 2v_xv_y + 2v_xv_z + 2v_yv_z}}$|⁠. From this, a ‘cold’ circular orbit in an ordered disc has |$\lambda _z = 1$|⁠. A ‘hot’ box or radial orbit has |$\lambda _z = 0$|⁠. Following Santucci et al. (2022), we define counter-rotating orbits as having |$\lambda _z \le -0.25$|⁠, hot orbits as |$-0.25 \lt \lambda _z \le 0.25$|⁠, warm orbits as |$0.25 \lt \lambda _z \le 0.8$|⁠, and cold orbits as |$\lambda _z \gt 0.8$|⁠. An example of the circularity space can be seen in Fig. 10 for MAGPI 1525170222, which is dominated by warm orbits. From the cuts defined above, we calculate the orbital fractions for the sample within |$R_{\mathrm{e}}$|⁠.

We examine the MAGPI sample for correlations between the orbital fraction (hot, warm, cold, and counter-rotating) within a half-light radius and other galaxy properties. We find a significant correlation between the intrinsic shape of a galaxy and the hot orbit fraction, shown in Fig. 11. We parametrize the intrinsic shape of a galaxy through its triaxiality parameter, defined as |$\mathrm{T}_{R_\mathrm{e}} = (1 - p^2)/(1 - q^2)$|⁠. The shape of a galaxy is not held constant with radius in Schwarzschild models, and so we take the values of p and q at the half-light radius for the triaxiality calculation. We find that the triaxiality of the sample also correlates strongly with the orbital anisotropy at the half-light radius, where we define the anisotropy as

where r is the intrinsic galaxy radius and |$\sigma _t = \sqrt{(\sigma _{\phi }^2 + \sigma ^2_{\theta })}/2$|⁠, with |$\sigma _r$|⁠, |$\sigma _{\phi }$|⁠, and |$\sigma _{\theta }$| the radial, azimuthal, and polar angular velocity dispersions in spherical coordinates, respectively. We choose this definition of anisotropy based on the results of Thater et al. (2022a), who found that the velocity ellipsoids for NGC 6958 are more closely aligned with spherical coordinates than the commonly adopted cylindrical coordinates. In our MAGPI sample, radially anisotropic (⁠|$\beta _r \gt 0$|⁠) objects are likely to be intrinsically triaxial and have mainly hot (box or radial) orbits, with few cold (short-axis tube) orbits.

We show the anisotropy profiles for the sample in Fig. 12, plotted between the PSF FWHM and maximum kinematic extent of the data, where we denote the anisotropy calculated at the half-light radius as |$\beta _{R_{\mathrm{e}}}$|⁠. A common practice for constructing dynamical and lensing models of galaxies is to either hold orbital anisotropy as a global, albeit free, constant (Cappellari et al. 2013a; Poci et al. 2017; Bellstedt et al. 2018; Derkenne et al. 2023), or to assume the systems are isotropic (Auger et al. 2010; Ruff et al. 2011; Sonnenfeld et al. 2013). For the MAGPI sample, we find instead that only the massive galaxies in the sample, with |$\log _{10}M_{\star }/{\rm M}_{\odot } \gtrsim 11.2$| have approximately constant or mildly increasing radial anisotropy profiles as a function of radius. There is remarkable variety in the sample, with some objects showing a dramatic reduction from radially anisotropic (pressure supported) to tangentially anisotropic (rotation supported) systems with radius, and others instead increasing from tangential anisotropy within |$R_{\mathrm{e}}$| towards more isotropic values at larger radii. The central values of anisotropy for the sample are also highly diverse. That the orbital structure is not constant with radius is consistent with our understanding of galaxies having internal photometric structures such as bulges and discs, which clearly must reflect different orbital families.

In Fig. 13 we present the intrinsic shapes of the MAGPI sample on the plane of possible ellipsoid shapes in relation to the triaxiality parameter (de Zeeuw & Franx 1991). We find that the majority of the 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. We construct the 3D, ellipsoidal intrinsic shape for MAGPI object 1 525 170 222 in Fig. 14, which shows how the intrinsic shapes of the best-fitting Schwarzschild model translates to an ellipsoid as seen from different viewing angles. Only 3/22 MAGPI galaxies are formally consistent with oblate spheroids (triaxiality parameter |$\mathrm{T}_{R_\mathrm{e}} \lt 0.1$|⁠). The remainder show a spread of triaxiality, with 10/22 consistent with mild triaxiality (⁠|$0.1 \le \mathrm{T}_{R_\mathrm{e}} \lt 0.3$|⁠), and 9/22 having values of |$\mathrm{T}_{R_\mathrm{e}}$| spanning |$0.3-0.7$|⁠. There are no prolate objects.

Unlike what we find for the MAGPI sample, Santucci et al. (2022) found that 73 per cent of their SAMI sample were oblate spheroids, with only 8 per cent having triaxiality |$\gt 0.3$|⁠, also inferred from Schwarzschild dynamical models. Emsellem et al. (2011) found the vast majority of early-type systems in the |$\mathrm{ATLAS}^{\mathrm{3D}}$| Survey are consistent with being a family of oblate spheroids, with only 12 per cent likely to exhibit mild triaxiality; however, we note this analysis makes inferences from the observed proxy-for-spin parameter and apparent ellipticity rather than dynamical models directly. Additionally, Emsellem et al. (2011) did not explore triaxiality as a free parameter, but considered how consistent the assumption of strict axisymmetry was with the data. There is also a mass dependence to the inferred intrinsic shapes from that work, with lower mass systems being more likely to be regular in their rotation with oblate, disc shapes. Similarly, Jin et al. (2020) found a mass dependence with the intrinsic shapes of galaxies, with more massive galaxies more likely to be prolate.

Some recent results are pointing to galaxies exhibiting triaxial intrinsic shapes rather than axisymmetric ones, in agreement with our MAGPI results. Initial results from Schwarzschild modelling of |$\mathrm{ATLAS}^{\mathrm{3D}}$| galaxies shows that the majority are mildly triaxial or triaxial in shape, with only a few classically oblate systems (Thater et al. 2023a). Intrinsic shapes were inferred from kinematic alignment and ellipticities for 90 galaxies in the MASSIVE Survey (Ma et al. 2014), finding most galaxies in the sample are mildly triaxial (Ene et al. 2018). Complex intrinsic shapes are also observed in the Magneticum simulations, with Valenzuela et al. (2024) finding most galaxies are triaxial. However, that work also found a high precedence of prolate galaxies, of which we have none in the presented MAGPI sample.

5.4 Systematics between Schwarzschild and Jeans modelling

The MAGPI galaxy sample presented here is a subset of the 30 galaxies with total density slopes measured using anisotropic Jeans models and a power law potential by Derkenne et al. (2023), which allows for an examination of the systematics between the two modelling methods on the same objects. Jeans dynamical modelling, under the assumption of a power law total density profile, is a powerful method especially when considering objects at significant redshift. When deriving the total density profile of a galaxy, the absolute calibration of the surface brightness is unimportant; only the relative scaling between MGE components matters to act as the luminous density tracer for the stellar kinematics. As such, assuming a total potential means the models are agnostic to many redshift effects and sources of systematic errors, such as cosmological surface brightness dimming, and how stellar light corresponds to stellar mass. Alternatively, building the gravitational potential from the stellar luminosity density and a separate dark matter halo, as presented in this work, requires more careful calibration of the stellar components (and assumptions around the IMF). Details of the Jeans modelling method can be found in Cappellari (2008), with specific details of how the MAGPI sample was modelled with the Jeans method in Derkenne et al. (2023). In Fig. 15 we show the total density slope comparison of the Jeans and Schwarzschild models, with objects coloured by their triaxiality.

We note a few differences between the modelling approaches. In the Jeans modelling approach used in Derkenne et al. (2023), we assumed axis-symmetry, a global but free orbital anisotropy, a cylindrically aligned velocity ellipsoid, and NFW-like double-power law with a free inner slope. This power-law forms the total gravitational potential, in that we made no attempt to disentangle the baryonic and dark components. In the Schwarzschild method used in this work, the models are fully triaxial and orbital structures can vary freely due to the three-integral orbit-based superposition method of the models. Furthermore, instead of assuming a total density profile shape, we assume the dark matter profile and measured the stellar mass profiles from a combination of photometry and spectral stellar population mass-to-light ratios. Due to this method, it is possible for the fitted total mass density profile to deviate from a power law.

From the direct comparison of these two methods, we can see that in general rank is preserved: steep slopes as measured by Jeans models remain steep when measured using Schwarzschild models, and shallow slopes remain comparatively shallow (with a notable outlier). However, the slopes generally do not agree within the errors. The discrepancy cannot be explained by projection effects, as there is no trend between the triaxiality of the sample and the offsets from the 1–1 line. There is only a mild (and non-significant) correlation between the slopes measured with different methods, with a Pearson correlation coefficient of 0.42 and p-value of 0.05. By contrast, the dynamical masses with an effective half-light radius are well recovered, with almost all the dynamical masses agreeing within the estimated uncertainties, as shown in Fig. 16. Despite a direct galaxy-by-galaxy comparison being inconsistent, the median total density slopes measured by the two methods agree well. From the Schwarzschild models in this work we obtain a median density slope of |$\gamma _{\mathrm{tot}} = -2.27 \pm 0.08$|⁠, and for the overlapping Jeans sample we find |$\gamma _{\mathrm{tot}} = -2.26 \pm 0.04$|⁠.

As argued by Etherington et al. (2023) when comparing a lensing-only analysis of density slopes with a lensing and dynamics approach, a 1–1 correspondence should be expected if the true density profile of the systems are well represented by a power law. That the slopes do not agree, for the most part, within the errors, is evidence that the MAGPI systems deviate from strict power-law potentials. The extra freedom of the Schwarzschild models, in terms of both modelling assumptions and the gravitational potential, yield different slopes to the Jeans derived ones. Interestingly, the variance of the Jeans derived slopes are much lower than that of the Schwarzschild models. The Schwarzschild sample total slopes have a standard deviation of |$0.30 \pm 0.03$|⁠, whereas the Jeans total slopes have a standard deviation of |$0.16 \pm 0.03$|⁠. We suggest this is because the total density profile shape is assumed to be a power law in the Jeans models, but for the Schwarzschild models the total profile is some free combination of the stellar mass profile and a dark matter halo. The additional freedom allowed by the gravitational potential used for the Schwarzschild models permits a larger radial variation of intrinsic density profiles to be fitted. Assuming this diversity is present in real galaxies, this naturally gives rise to larger galaxy-to-galaxy variation in the inferred density slopes.

Recently, a sample of over 10 000 MaNGA galaxies was dynamically modelled by Zhu et al. (2023a) using multiple Jeans method approaches, all of which had a gravitational potential constructed by including the stellar luminosity information with the possible co-addition of different dark halo profile types (in a similar vein to the approach used in this work). Zhu et al. (2023a) found that the total slopes, measured using the different assumed dark matter haloes, were consistent. However, for the Jeans dynamical models constructed by Poci et al. (2017) on the |$\mathrm{ATLAS}^{\mathrm{3D}}$| sample, the mean total slope as measured by assuming a power law total density profile structure was not consistent with the total slope as measured by building the density profile from the de-projected stellar luminosity and an assumed dark halo (comparing the Model I and Model III total density profile slope means in that work).

From the above, we suggest it is not the dynamical modelling method (Jeans or Schwarzschild) per se that causes the systematic difference observed here, but rather the assumption of the structure of the total gravitational potential. It may also be that in assuming a power-law density profile structure when measuring total density slopes, as is done in various dynamical and lensing approaches, artificially reduces the population variance observed. A thorough exploration of the systematics between the modelling techniques and assumptions would be required to fully understand this issue, and is beyond the scope of this work.

6 DISCUSSION

The primary goal of this work is to use the high quality kinematic data from the MAGPI Survey to investigate the stellar and total density profiles of massive galaxies using a highly general dynamical modelling approach at a cosmologically significant look-back time. We have shown no correspondence between central dark matter fractions and stellar density profiles, whilst dark matter is highly correlated to the gradient of the total density profile. These results are in contrast gravitational lensing studies, which consistently report low intrinsic scatter in measured total density slope values, pointing to a homogeneity for (early-type) galaxy mass density structures.

A limitation of our approach has been to use a |$M_{200}-c$| coupled NFW halo, which effectively means that the only freedom in the total density profile is the combination of a globally scaled stellar mass profile and the dark matter fraction. As seen with some galaxies in the sample, it is possible to reproduce the higher order kinematics with only the stellar profile and negligible contributions from the dark halo within the data range probed. Because we have not implemented the effects of IMF, if the true total density profile is even steeper than the measured stellar profiles, our models are bounded at having no dark matter. As a result, it is possible that constraints caused by our modelling approach has impacted the resulting total density slopes returned. In our comparison of the Schwarzschild models against Jeans dynamical models of the same sample in Section 5.4, we show that our Schwarzschild model slopes are more variable than those returned using Jeans models with a completely free power-law slope on the total density profile and no assumption on the form of the stellar or dark matter profiles. This suggests our Schwarzschild models are not too restrictive in their assumptions. Furthermore, if our models lacked significant freedom we would expect to see departures of the modelled kinematics to the observed ones. Instead, we note our models are able to reproduce the higher order kinematic features of many galaxies in the sample, such as MAGPI object 1530197196. This aside, a modelling approach with, for example, a generalized NFW dark matter halo with free inner slope and the inclusion of variable IMF would be invaluable in testing the bulge–halo conspiracy in the future.

A recent meta-analysis of gravitational lensing works by Etherington et al. (2023) examined 48 lensing systems from the Sloan Lens ACS (Bolton et al. 2006; Auger et al. 2010), BOSS Emission Line Lensing (Brownstein et al. 2011), BELLS GaLaxy-Ly|$\alpha$| Emitter Systems (Shu et al. 2016), Strong Lensing Legacy Survey (Gavazzi et al. 2012), and the Lenses, Structure, and Dynamics (Treu & Koopmans 2004) surveys, with total density slopes calculated with two different methods: using a lensing-only approach as presented in Etherington et al. (2022), and a joint lensing and dynamics approach, where the stellar velocity dispersion is used to constrain the mass density profile. Consistent with the bulge–halo conspiracy, Etherington et al. (2023) found a mean, almost isothermal, total density slope for the sample of |$\gamma _{\mathrm{tot}} = -2.075^{+0.023}_{-0.024}$| with an intrinsic scatter of |$\sigma _{\gamma } = 0.172^{+0.022}_{-0.032}$|⁠, assuming a power-law density distribution. This scatter is far smaller than what we find for the MAGPI sample, of |$\sigma _{\gamma } = 0.30 \pm 0.03$|⁠.

It is unclear what drives the differences between the variation in total slope values reported in this work and those found by lensing studies, although it is possibly driven by how the gravitational potential is incorporated in the lensing models compared to the more general potential we have allowed for in this work. Furthermore, Etherington et al. (2023) note that lensing-only analyses are sensitive to the local slope at the Einstein radius, which is more likely to occur at a transitional region between the stellar dominated centre and dark matter dominated outskirts, and so lensing derived slopes are more likely to indicate possible deviations from a power-law density structures for galaxies. Lensing in combination with dynamical methods tends to average the slope between the inner and outer galaxy regions, potentially fueling the conspiracy, as deviations from a power-law structure may be averaged out. Methodologically this might resolve the conspiracy as lensing-only studies are conducted across larger redshift baselines.

The bulge–halo conspiracy has also been explored with 5688 Jeans dynamical models from local Universe MaNGA galaxies (Li et al. 2024) with visual checks performed to ensure model quality. Specifically, Li et al. (2024) explored the correlation between the stellar density slopes, total slopes, and dark matter fraction within the effective radius of their sample. As found in this work, they found a steep stellar density slope corresponds to a steep total density slope, whereas a high dark matter fraction corresponds to a more shallow total density slope (see their figs 6 and 7). Li et al. (2024) found no evidence of the conspiracy for their full sample, in that there is a wide range of slopes observed (between |$-1.2$| and |$-2.8$|⁠). For early-type galaxies with high velocity dispersions, Li et al. (2024) found the median total density slope does trend towards near isothermal values, but that there is still a correspondence between stellar density profiles and total density profiles. That is, dark matter does not compensate to create a conspiracy. Instead, Li et al. (2024) conclude that the outcome of a median isothermal total density slope for early types is simply a coincidence resulting from the mass growth of galaxies.

In an accompanying work to Li et al. (2024), Zhu et al. (2023b) found that galaxies with older stellar populations tend to have steeper total density slopes, and galaxies with young stellar populations have shallow total density slopes (in particular, see their fig. 8). This trend is also coupled to the central velocity dispersions, such that galaxies with high velocity dispersions have comparatively steeper slopes than galaxies with lower velocity dispersions, as initially seen by Poci et al. (2017). Zhu et al. (2023b) further find that the scatter in total density slopes also decreases with increasing central velocity dispersion. The MAGPI galaxies typically have high velocity dispersions, with a range of |$97-328\, \, \mathrm{kms}^{-1}$| and a mean of |$210\, \, \mathrm{kms}^{-1}$|⁠, which is on the high end of the MaNGA sample distribution. Even with a high average velocity dispersion, we still find large scatter in the MAGPI total density slopes.

In exploring the conspiracy with simulation data, Remus et al. (2013) found a power-law approximation to be a reasonable fit for total mass density profiles (see their fig. 2). In the simulations used, merger events caused the total density slope to flatten towards isothermal values, with steep total density slopes likely for systems that have undergone few mergers and therefore have a high fraction of in situ stars. Elliptical galaxies with almost no mergers in their history tended to retain steep density slopes of |$\gamma \sim -3$|⁠, whilst progressive dry mergers cause the slope to converge nearer to |$\gamma \sim -2$| across cosmic time, coupled with increasing dark matter fractions (see Remus et al. 2013, in particular, their fig. 12).

Minor mergers in particular are efficient in building galaxies with high dark matter fractions, as their stars are incorporated into the main progenitor at relatively large radii where dark matter is more dominant, while also reproducing the observed size growth for both elliptical galaxies (Hilz et al. 2012; Hilz, Naab & Ostriker 2013) and disc galaxies (Karademir et al. 2019). In this sense, an isothermal total density slope of |$-2$| might act as a natural attractor due to progressive mergers, as noted by Remus et al. (2017), although values both considerably steeper and more shallow than this are found in both simulations and observations (Poci et al. 2017; Remus et al. 2017; Bellstedt et al. 2018; Derkenne et al. 2023; Zhu et al. 2024). Due to the relation between dark matter fractions, mergers, and in situ stars, Remus et al. (2017) also found an anticorrelation between the dark matter fractions and fraction of in situ stars with the Magneticum simulations. Instead of a conspiracy, the evolution of the galaxy’s total density profile towards shallower values, from steep values in the early Universe, indicates the (stochastic) merger history of a galaxy. Following this line of argument, we would naturally expect the variation of total density profiles we see with MAGPI galaxies, instead of a homogeneous sample of near-isothermal total density slopes. As gas accretion becomes less common at more recent times (de Ravel et al. 2009; Oser et al. 2010; Tacconi et al. 2010), the theorized evolution is for total density slopes to become more shallow, although we note the reverse of this evolution is seen by lensing samples (Etherington et al. 2023).

Finally, there are clearly systematic effects at play between building the gravitational potential under the assumption of a power-law total density profile with Jeans dynamical modelling, and a Schwarzschild modelling approach with a potential constructed from the combination of a stellar mass profile and separate dark matter halo, as shown in Section 5.4. These systematics are despite the fact both methods use almost identical data inputs (MGEs and 2D stellar kinematics). We can only assume that there are further systematic effects to be considered between dynamical modelling in general and lensing, and so advocate as future work the joint analysis of the same systems with the available, complementary, techniques, as done by Poci & Smith (2022), but extended to larger samples of galaxies. Only in this way will we ensure that comparisons across studies, and the conclusions drawn from them, are fairly made.

7 SUMMARY AND CONCLUSIONS

We present Schwarzschild orbit-based triaxial models for a sample of 22 well-resolved galaxies from the MAGPI Survey, exploring in detail the properties of massive galaxies during the Universe’s middle ages. The sample spans stellar masses between |$\log _{10}(M_{\star }/\mathrm{M}_\odot) = 10.4$| and |$\log _{10}(M_{\star }/\mathrm{M}_\odot) = 11.6$|⁠, with a median stellar mass of |$\log _{10}(M_{\star }/\mathrm{M}_\odot) = 11.2$| The MAGPI Survey uses MUSE/VLT in the wide-field mode coupled with ground layer-adaptive optics, resulting in rich and highly resolved fields at a target redshift of |$z \sim 0.3$|⁠. We use 2D, 4-moment stellar kinematics from full spectrum fitting and multi-Gaussian expansions of the stellar light, in combination with spectral mass-to-light ratio maps, to construct the Schwarzschild models. We use the state-of-the-art dynamite code (Jethwa et al. 2020), allowing for fully triaxial Schwarzschild models, and make no assumptions on orbital structure. To our knowledge, this is the first time Schwarzschild modelling has been applied on a sample of galaxies at a cosmologically significant redshift.

Instead of assuming a power law for the gravitational potential, as commonly done in both lensing and dynamical studies of galaxies, we instead measure the stellar surface brightness profile and transform it to a stellar mass profile using spectral stellar mass-to-light ratio maps from MAGPI data. To build the potential we combine the stellar mass profile with a freely scaled NFW dark matter halo, in order to test the so-called bulge–halo conspiracy. Numerous observations have recovered near isothermal total (baryonic and dark) density profiles for galaxy populations, with small intrinsic scatter. This finding requires that the dark and luminous mass components must somehow conspire to yield remarkable homogeneity in observed galaxy mass structures despite varied stellar light structures. From our models we obtain both stellar mass profiles and constrained NFW profiles. In addressing the bulge–halo conspiracy, we find:

• The bulge–halo conspiracy predicts that the variety of structures observed for the stellar component and theorized for the dark component are larger than the variation of the combined, approximately isothermal, profile. Fitting a flexible potential with separate realistic stellar and dark mass profiles, we find a wide range of total density profile slopes, spanning approximately between |$\gamma = -2.7$| and |$\gamma = -1.6$| for |$\rho \propto r^{\gamma }$| and |$\gamma \lt 0$|⁠. Indeed, the dispersion of the total density profile slopes is greater than that for the stellar profile slopes alone, comparing |$\sigma _{\mathrm{tot}} = 0.30 \pm 0.03$| to |$\sigma _{\star } = 0.19 \pm 0.02$| for the total and stellar components, respectively. This finding is inconsistent with the bulge–halo conspiracy.

• While we find a tight correlation (Pearson |$r \gt 0.8$|⁠) between the total density profile slope and the central dark matter fraction, and a moderate correlation between the total density profile slope and the stellar profile slope, we find no significant correlation between the central dark matter fraction and the stellar profile slope. The lack of correlation between the stellar mass structure and the dark matter fraction is also inconsistent with the bulge–halo conspiracy, as the dark matter should systematically adjust with differing stellar profiles if the conspiracy were valid. However, as predicted by the Magneticum simulations, MAGPI galaxies with higher central dark matter fractions have shallower total density profiles. The relation between the dark matter fraction within the half-light radius and density slopes is |$\gamma _{\mathrm{tot}} = (1.3 \pm 0.2) f_{\mathrm{DM}}-(2.44 \pm 0.04)$|⁠, consistent with Magneticum results (Remus et al. 2017).

• From a comparison of total mass density profile slopes from the Schwarzschild and Jeans modelling methods, we suggest that assuming a power law density profile may artificially reduce the observed intrinsic scatter in population density profile slopes.

Our detailed Schwarzschild models also allow us to explore other intrinsic properties of the sample, including orbit types and structures, and the intrinsic shapes of massive galaxies. In exploring the MAGPI sample with these models, we find:

• The median 1 |$R_{\mathrm{e}}$| dark matter fraction for the sample is 10 per cent with a standard deviation of 19 per cent, which is similar to local Universe survey results from MaNGA, but lower than reported from the SAMI Survey.

• The dark matter profile does not tend to dominate at any specific radius across the sample. For some galaxies the dark matter density exceeds the stellar density within the half-light radius, while for others the dark matter only exceeds the stellar mass density at more than 10 |$R_{\mathrm{e}}$|⁠. How baryonic and dark matter is distributed within a galaxy is therefore highly variable.

• In exploring the intrinsic shapes of our MAGPI sample at 1 |$R_{\mathrm{e}}$|⁠, we find that only 3/22 are formally consistent with oblate spheroids (triaxiality parameter |$\mathrm{T}_{R_\mathrm{e}} \lt 0.1$|⁠). The remainder show a spread of triaxiality, with 10/22 consistent with mild triaxiality (⁠|$0.1 \le \mathrm{T}_{R_\mathrm{e}} \lt 0.3$|⁠), and 9/22 having values of |$\mathrm{T}_{R_\mathrm{e}}$| spanning |$0.3-0.7$|⁠.

• The triaxiality of a galaxy correlates strongly with the fraction of hot orbits within an effective radius, as well as the orbital anisotropy of the galaxy. Galaxies that exhibit radial anisotropy are likely to be strongly triaxial in intrinsic shape.

• The anisotropy profiles of the sample are also highly variable. Massive galaxies with radial anisotropy tend to have constant or mildly increasing anisotropy with radius, whereas less massive galaxies that are rotation supported can display either increasing or decreasing tangential anisotropy with radius. The central anisotropy within the half-light radius is also variable for the sample.

From our results, we conclude that the internal mass structure of galaxies tends to be non-homogeneous, varying both between galaxies and within galaxies. Our analysis acts as a caution in making strong assumptions when constructing dynamical (or lensing) models of galaxy internal mass distributions and structure. The more closely we examine galaxies, the more rich in structure they reveal themselves to be – and with the increasing observing power of modern integral field units, we are continuously improving our ability to explore distant galaxies in unprecedented detail.

The detailed orbital structures revealed in this analysis of |$z \sim 0.3$| galaxies is only with a subsample of MAGPI Survey targets, with tens more objects with sufficient data quality for such modelling expected from the completed survey. The advance of both space-based and ground-based instruments, for example JWST and extremely large telescopes such as the E-ELT, mean that it is possible for this kind of analysis to be pushed yet further afield, connecting the local Universe to its middle ages and beyond. The discovery of |$z \sim 2$| lenses presents the opportunity of joint lensing and dynamics studies at lookback times of 10 Gyr (e.g. Mercier et al. 2024; van Dokkum et al. 2024). The instrumentation available is starting to allow for consistently applied methods across unprecedented cosmic time baselines.

Acknowledgement

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 ID 179.A-2004, ID 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 number CE170100013. SMS acknowledges funding from the Australian Research Council (DE220100003). CF is the recipient of an Australian Research Council Future Fellowship (project number FT210100168) funded by the Australian Government. CL, JTM, and CF are the recipients of ARC Discovery Project DP210101945.

Software citations: This work uses the following software packages:

astropy (Astropy Collaboration 2022)

dynamite (Jethwa et al. 2020)

galfit (Peng et al. 2002)

ltsfit (Cappellari et al. 2013a)

matplotlib (Hunter 2007)

mgefit (Cappellari 2002)

numpy (Harris et al. 2020a)

pafit (Krajnović et al. 2006)

pandas (McKinney 2010)

plotbin

ppxf (Cappellari 2008, 2020)

profound (Robotham et al. 2018)

prospect (Robotham et al. 2020)

vorbin (Cappellari & Copin 2003)

scipy (Virtanen et al. 2020)

zap (Soto et al. 2016)

DATA AVAILABILITY

The MAGPI spectral data underlying this work are available from the ESO Science Archive Facility: http://archive.eso.org/cms.html.

Footnotes

References

APPENDIX A: IMPACT OF DITHERING ON MODEL ORBIT WEIGHTS

We chose not to use dithering for our Schwarzschild models. Dithering perturbs the initial starting positions of the orbits in the (E, |$L_z$|⁠, |$I_3$|⁠) space to regularize the final weighted orbits. Dithering increases the orbit library size by a factor of |$N^{3}$|⁠, where N is the chosen dithering of, for example, 3 or 5. For a technical description and motivation for dithering see van den Bosch et al. (2008). Due to the increased orbit library size, dithered models take significantly more computational time to run than undithered models. We explored whether dithering made a substantial impact on the returned model orbit weights by running a test MAGPI galaxy (ID = 1203060081) without dithering, and then with dithering set to 3. We found that there was no significant change in the orbit weights for the test galaxy, however the undithered model took of order 400 h to run, whereas the dithered model took 8000 h of equivalent compute time. We show the orbital weights for the dithered and undithered best-fitting models in Fig. A1. We conclude that for MAGPI-like data quality, dithering is not necessary.

APPENDIX B: SAMPLE SCHWARZSCHILD MODEL KINEMATIC MOMENTS

In Fig. B1 we show the remaining 20 Schwarzschild models for the MAGPI sample used in this work, with models for MAGPI objects 1205093221 and 1525170222 shown in Figs 4 and 5, respectively. For each model we show the surface brightness and line-of-sight kinematic moments up to |$h_4$|⁠. The panels are ordered by MAGPI field number, being the first four digits of the MAGPI ID. Gaps in the fields are for masked regions, e.g. for MAGPI object 1208197197.

APPENDIX C: TESTING THE IMPACT OF VARIABLE IMFS

In our modelling methods we have so far left the thorny question of variable IMFs unanswered. This is partly because the data quality and spectral lines required (e.g. TiO₂ and Na i) to effectively measure the IMF (Spiniello et al. 2012) are unavailable with MAGPI data due to redshift and signal-to-noise. However, we still made an attempt to understand what impact a variable IMF within a galaxy might have on our results. A non-universal IMF between objects is not a concern because of the free re-scaling of all the mass components in the Schwarzschild models, parametrized as |$\Upsilon _r$|⁠, which essentially adjusts the global IMF. Therefore, although all spectroscopic mass-to-light ratios are calibrated to a Salpeter IMF, the mass scaling in the Schwarzschild models is left free and not tied to any particular IMF.

Although a difference of IMF between galaxies is not a concern, an additional gradient on the stellar mass-to-light ratios due IMF structures within galaxies could have an effect on the returned dark matter fractions. Our models attempt to accurately reconstruct the stellar mass profile by use of detailed surface brightness models and stellar population mass-to-light ratios, fitting for age and metallicity variations (which are constrained by the MAGPI data). The mass-to-light ratios adjust the gradient of the stellar surface brightness, which in turn de-distributes the stellar mass. A radially varying IMF would introduce another gradient to the stellar mass profile.

Observations support a non-global galaxy IMF, in particular the 2D resolved view of Fornax3D (Sarzi, Iodice & Fornax3D Collaboration 2022) galaxy FCC 167 showing a dwarf-rich population in the galaxy’s centre, which correspond to higher mass-to-light ratios in the central regions compared to the outskirts (Martín-Navarro et al. 2019). A wealth of other studies support this view of IMF gradients (Martín-Navarro et al. 2014; Conroy et al. 2017; van Dokkum et al. 2017). If we assume that in general massive galaxies tend to be dwarf rich in their centres, then a variable IMF would tend to make the stellar mass profiles steeper. To see what effect this has on Schwarzschild-derived orbits and dark matter fractions, we construct a simple test case. For MAGPI object 1525170222 we artificially steepen the stellar mass profile by increasing the innermost mass-to-light ratio by 15 per cent and decreasing the outermost by the same, and re-run the Schwarzschild models. This adjustment level was chosen because the fraction of mass contained in the dwarf-population FCC 167 changes from 0.66 to 0.5 from the galaxy centre to the outskirts.

For MAGPI object 1525170222 we find the dark matter fraction increases for our test case, to |$0.50^{+0.11}_{-0.16}$| from an original |$0.44^{+0.11}_{-0.06}$|⁠. The global mass scaling |$\Upsilon _r$| parameter correspondingly decreased to |$0.47^{+0.08}_{-0.08}$| from |$0.53^{+0.08}_{-0.04}$|⁠. The reduced |$\chi ^2$| value of the IMF test case and the original model are comparable: 2.45 and 2.43, respectively. Although this is a single test case, we use this example to illustrate the possible uncertainties, of order 10 per cent, on our reported dark matter fractions due to possible IMF gradients.