Dark matter measurements combining stellar and HI kinematics: 30% $1-\sigma$ outliers with low dark matter content at $5R_\mathrm{e}$

We construct the Schwarzschild dynamical models for 11 early-type galaxies with the SAURON and Mitchell stellar IFUs out to $2-4 R_\mathrm{e}$, and construct dynamical models with combined stellar and HI kinematics for a subsample of 4 galaxies with HI velocity fields out to $10 R_\mathrm{e}$ obtained from the Westerbork Synthesis Radio Telescope, thus robustly obtaining the dark matter content out to large radii for these galaxies. Adopting a generalised-NFW dark matter profile, we measure an NFW-like density cusp in the dark matter inner slopes for all sample galaxies, with a mean value of $1.00\pm0.04$ (rms scatter $0.15$). The mean dark matter fraction for the sample is $0.2$ within $1 R_\mathrm{e}$, and increases to $0.4$ at $2 R_\mathrm{e}$, and $0.6$ at $5 R_\mathrm{e}$. The dark matter fractions within $1 R_\mathrm{e}$ of these galaxies are systematically lower than the predictions of both the TNG-100 and EAGLE simulations. For the dark matter fractions within $2 R_\mathrm{e}$ and $5 R_\mathrm{e}$, 40% and 70% galaxies are $1-\sigma$ consistent with either the TNG-100 or the EAGLE predictions, while the remaining 60% and 30% galaxies lie below the $1-\sigma$ region. Combined with 36 galaxies with dark matter fractions measured out to $5 R_\mathrm{e}$ in the literature, about 10% of these 47 galaxies lie below the $3-\sigma$ region of the TNG-100 or EAGLE predictions.


INTRODUCTION
In the context of cold dark matter cosmology (Blumenthal et al. 1984;Davis et al. 1985), galaxies form in the centres of dark matter halos.The distribution of dark matter within galaxies is affected not only by the characteristics of dark matter particles (e.g.Kaplinghat et al. 2016;Robertson et al. 2018) but also by the presence of baryonic matter during galaxy formation.Star formation is regulated by a series of feedback mechanisms, such as active galactic nucleus feedback (e.g.Di Matteo et al. 2005;Hopkins et al. 2006;Debuhr et al. 2010) and supernova feedback (e.g.Strickland & Heckman 2009;Ostriker & Shetty 2011).These processes transfer energy to dark matter particles through dynamical friction, thus altering the distribution of dark matter within galaxies (e.g.Navarro et al. 1996a;Read & Gilmore 2005;Cole et al. 2011).In shallow gravitational potential wells, such as those of dwarf galaxies, baryonic feedback has a noticeable effect on the dark matter distribution, which can explain the observed "cored" dark matter distribution (Read et al. 2019).However, for massive galaxies, baryonic feedback might not be enough to significantly modify the distribution of dark matter (Tollet et al. 2016;Hopkins et al. 2018;Bullock & Boylan-Kolchin 2017).
Most studies that assess the dark matter mass fraction of large samples with dynamical modelling rely on stellar kinematic data from spatially resolved integral field unit (IFU) spectrographs, such as Atlas 3D (Cappellari et al. 2011), CALIFA (Sánchez et al. 2012), SAMI (Croom et al. 2012), and SDSS-IV MaNGA (Bundy et al. 2015).However, these data only cover the brightest central 1-2 ef-fective radii ( e ) of galaxies, where the contribution of dark matter is minimal.Therefore, using only these stellar kinematic data to measure the properties of dark matter halos (including the dark matter fraction and density profile) can lead to considerable uncertainties (Cappellari et al. 2013;Zhu et al. 2018;Santucci et al. 2022;Zhu et al. 2023).
In order to obtain more precise measurements of dark matter, we require kinematic data extending to the outskirts of galaxies.Although stellar kinematics using extended IFU spectrographs and long-slit spectrographs (Weĳmans et al. 2009;Forestell & Gebhardt 2010;Cappellari et al. 2015;Boardman et al. 2016;Comerón et al. 2023) can provide such measurements, the amount of data is limited by the required long exposure time.Therefore, other extended tracers are often used to measure dark matter distributions.Planetary nebulae (PNe) and globular clusters (GCs) are commonly employed for dark matter measurements independently (Romanowsky et al. 2003;Côté et al. 2001Côté et al. , 2003;;Deason et al. 2012;Alabi et al. 2016Alabi et al. , 2017;;Dumont et al. 2023), or in combination with stellar kinematics (de Lorenzi et al. 2008(de Lorenzi et al. , 2009;;Napolitano et al. 2009Napolitano et al. , 2011Napolitano et al. , 2014;;Das et al. 2011;Morganti et al. 2013;Zhu et al. 2014Zhu et al. , 2016;;Li et al. 2020;Yuan et al. 2022), mostly for the most massive early-type galaxies, where there are enough of these tracers to make measurements.On the other hand, the H i gas rotation curve is a powerful tracer of dark matter mass in late-and early-type galaxies (de Blok et al. 2008;Oh et al. 2011;Lelli et al. 2016).Combining H i gaseous kinematics and stellar IFUs can provide strong constraints on the dark matter fraction of galaxies from the inner to outer regions (Yang et al. 2020).
Current measurements of dark matter fractions in galaxies, in comparison to cosmological numerical simulations (Lovell et al. 2018), demonstrate the following: (1) Different cosmological numerical simulations do indeed produce varying dark matter fractions within 5 e ; (2) Observations that reach out to 5 e are inadequate and display a considerable scatter, with few measurements in galaxies with stellar masses below 10 10  ⊙ ; (3) Dark matter fractions measured using stellar kinematics only covering 1 − 2 e may be affected by systematic biases due to the degeneracy between dark matter and baryonic matter.Consequently, the existing observational data are not enough to impose stringent constraints on cosmological numerical simulations.To effectively constrain cosmological numerical simulations, precise measurements of dark matter fractions are needed from large samples extending to large radii, while breaking the degeneracy between dark matter and baryonic matter.
This paper aims to measure the dark matter profiles of a selection of intermediate mass early-type galaxies with extended kinematics.We structure this paper as follows: Section 2 introduces the sample and data used in our dynamical models.Section 3 outlines our dynamical modelling method.Section 4 presents our results, including the dark matter fractions and the total density slopes.In Section 5, we discuss our results and compare them with simulations and observations in the literature.Finally, Section 6 summarises our work.

SAMPLE AND DATA
We choose the sample of 12 ETGs presented in (Boardman et al. 2017), with the exception of NGC 3626 due to its strong bar (Salo et al. 2015).Details of the 11 galaxies in our sample are shown in Table 1.This sample has a stellar mass range of 10 10 − 10 11  ⊙ and a range in effective radius of 1 − 4 kpc, and includes 3 slowrotators dominated by random motions and 8 fast-rotators dominated by circular motions, as classified in Cappellari et al. (2011).We show the mass-size relation of the sample in Figure 1, and compare  Cappellari et al. (2013).The grey dots are all the galaxies in the Atlas 3D Survey (Cappellari et al. 2011).The gold and violet solid curves are the mass-size relations of the early-type galaxies in the TNG-100 and EAGLE simulations, and the corresponding dark and light shaded areas cover the 16th to 84th fractiles and 0.15th to 99.85th fractiles, respectively.
it with the full Atlas 3D sample.We also plot the mass-size relations of the early-type galaxies in the TNG-100 simulation (Naiman et al. 2018;Marinacci et al. 2018;Pillepich et al. 2018;Nelson et al. 2018;Springel et al. 2018) and the EAGLE simulation (Schaye et al. 2015;Crain et al. 2015), in which the  e of a simulated galaxy is obtained by projecting this galaxy along a random direction to obtain its surface brightness (Genel et al. 2018), and the early-type galaxies are defined as central galaxies with a specific star formation rate (sSFR) lower than 10 −11 /yr in both simulations.Cosmological simulations have not been able to accurately replicate the mass-size relation seen in observations, especially at the high-mass end (Lovell et al. 2018), which could lead to discrepancies when using  e to measure enclosed dark matter fractions.
We used the -band images obtained with the MegaCam instrument on the Canada-France-Hawaii Telescope (CFHT) as part of the MATLAS survey (Duc et al. 2015).This instrument has a wide fieldof-view of 1 degree 2 , and these images have a limiting magnitude of ∼ 29 mag arcsec −2 in the -band, which allows us to trace the surface brightness out to 5 e , beyond the range of our kinematic data sets.As the MATLAS -band images of most sample galaxies (all except NGC 2764, NGC 6798 and UGC 03960) are saturated in the galaxy centre, we also use their -band images taken from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) for the correction of saturation, which are stacked from short-exposure images to reach a limiting magnitude of 23.2 in the -band.We refer the reader to Chambers et al. (2016) and the references therein for further information on Pan-STARRS.
This sample of galaxies has two sets of stellar kinematic data.The central stellar kinematics were observed with the SAURON integral field spectrograph on the William Herschel Telescope (Bacon et al. 2001) in the Atlas 3D Survey (Cappellari et al. 2011).The SAURON integral field spectrograph has a field-of-view of 33 × 41 arcsec 2 with lenslets of 0.94 arcsec 2 and an instrumental velocity disper-Table 1. Basic information of sample galaxies.All values are taken from Cappellari et al. (2011), except that the K-band stellar mass  * , is estimated from M  according to Equation 2 in Cappellari (2013), and  e in kpc is calculated from  e in arcsec and distance.FR / SR stands for fast-rotators / slow-rotators as classified in Cappellari et al. (2011) sion of 105 km/s, and covers the central 1 e of the sample galaxies.
The spectra were Voronoi binned (Cappellari & Copin 2003) to a signal-to-noise ratio of 40, and fitted with the penalised pixel fitting method (pPXF; Cappellari & Emsellem 2004) to obtain the stellar kinematics (velocity, velocity dispersion and Gaussian-Hermite moments ℎ 3 , ℎ 4 ) 1 .The extended stellar kinematics was observed with the Mitchell integral field spectrograph (Hill et al. 2008) on the Harlan J. Smith telescope, as presented in Boardman et al. (2017).
The Mitchell integral field spectrograph has a large field-of-view of 1.68 × 1.68 arcmin 2 , a large fibre size of 2.08 arcsec and a maximum instrumental velocity resolution of 42 km/s, and covers 2 − 4 e for the sample galaxies.The stellar kinematics were obtained with the pPXF method which includes the upgrades described in (Cappellari 2017), on spectra that were Voronoi binned to a target signal-to-noise ratio of 20.All of the sample galaxies were observed with the Westerbork Synthesis Radio Telescope (WSRT) and showed extended H i discs or rings in the intensity maps (Serra et al. 2012).However, only a subsample of four galaxies shows regular rotating features in their H i velocity fields in deeper observations with the WSRT: NGC 2685, NGC 4203, NGC 5582 and NGC 6798 (see Appendix A for these velocity maps).The observing time for these follow-up observations was 9 × 12 hours for each galaxy.The raw data were processed and calibrated using the MIRIAD software package (Sault et al. 1995), resulting in a datacube with a spectral resolution of 16 km/s.We note that the observations of NGC 2685 and NGC 4203 were presented in Józsa et al. (2009) and Yıldız et al. (2015), respectively, and refer the readers to Yıldız et al. (2015) for more details about the data reduction.

METHODS
For the entire sample, we construct the Schwarzschild model with two-aperture stellar kinematics, the SAURON data and the Mitchell data.For the subsample with available WSRT data, we also build dynamical models combining two-aperture stellar kinematics and H i velocity fields: the stellar kinematics is modelled with 1 The stellar kinematics of NGC 1023 and NGC 2685 were observed and first presented in Emsellem et al. (2004), and these observations were re-reduced as part of the Atlas 3D Survey a Schwarzschild model, while the H i velocity fields are modelled as a thin ring misaligned with the stellar disc.In this section, we use NGC 4203 from the subsample to demonstrate our methods.We begin by constructing the gravitational potential of the galaxy.We then illustrate the stellar orbit library used in the Schwarzschild models and the misaligned gas model for the H i discs of the subsample.Subsequently, we select the best-fitting model and determine the 1 −  confidence levels.Finally, we obtain the mass budget, particularly the dark matter fraction, from the models.

Gravitational potential
The gravitational potential of the galaxy is contributed by three components: stars, a dark matter halo, and a central black hole.One key point to break the stellar-dark matter degeneracy is to describe the stellar mass distribution precisely.
The stellar mass distribution of a galaxy is equal to its surface brightness times its stellar mass-to-light ratio.For the sample galaxies with no saturation in the central region of the MATLAS -band image, we adopt the MATLAS -band image as the surface brightness.For the rest of the sample galaxies, we first match and calibrate the Pan-STARRS images referring to the MATLAS -band image; we then replace the saturated central region in the MATLAS -band image with the Pan-STARRS one to generate a combined image; we finally adopt the combined image as the surface brightness, and mask foreground stars and their saturation.We fit the surface brightness with the sum of multiple 2-dimensional Gaussian profiles (MGE; Emsellem et al. 1994;Cappellari 2002), and show the MGE fitting of NGC 4203 in the top panel of Figure 2.
The sample galaxies display strong gradients in their -band massto-light ( * /  ) ratio obtained from the SAURON spectra Poci et al. (2017) and the -band mass-to-light ratio ( * /  ) obtained from the Mitchell spectra (Boardman et al. 2017), which are obtained with the stellar population synthesis models.Therefore, we take the massto-light ratio gradient into account instead of assuming a constant mass-to-light ratio (Yang et al. 2020).We first rescale the Mitchell  * /  linearly to match the SAURON  * /  .We then generate a set of mass MGEs, ensuring that the corresponding mass-to-light ratio obtained from the mass and surface brightness MGEs is consistent with the observational one.We show the fitting of the mass-to-light gradient of NGC 4203 in the right panel of Figure 2.
The stellar mass and the corresponding mass-to-light ratio ob- tained through stellar population models are affected by the selection of the initial mass function (IMF) (Santini et al. 2012).To account for this, we introduce a scaling factor  to indicate the total mass scale due to the IMF:  = 1.0 is set to the stellar mass produced by the Salpeter IMF (Salpeter 1955), while the Kroupa IMF (Kroupa 2001) produces a stellar mass roughly with  = 0.6 (Speagle et al. 2014).
We then deproject the 2-dimensional MGE model to a triaxial MGE model with a set of viewing angles (, , ).The viewing angles are linked to the three ratios of the intrinsic shape of the galaxy ( , , ) (van den Bosch et al. 2008), which are the free parameters of the gravitational potential.We set the minor-to-major axis ratio  and the intermediate-to-major axis ratio  as free parameters to allow triaxiality, and fix  = 0.9999 to reduce model freedom.
We choose a spherical generalised-NFW profile (Navarro et al. 1996b;Zhao 1996) for the dark matter halo, .
(1) which contains four parameters:   is the scale density,   is the scale radius,  is the inner slope (when  = 1, it reduces to the NFW profile), while  controlling the turning point is fixed to 2. The central black hole is regarded as a point source with a softening length.Because it is not resolved by the kinematic data and has little effect on the result, we fix its mass at 10 6  ⊙ and its softening length at 10 −3 pc.
In total, we have six free parameters in the dynamical model: the total mass scale , the minor-to-major axis ratio  and the intermediateto-major axis ratio , the scale density   , the scale radius   and the inner slope  of the dark matter halo.

Model of stellar kinematics
We construct a Schwarzschild model for each set of parameters using the DYNAMITE package (Jethwa et al. 2020;Thater et al. 2022), which is a modified version of a triaxial dynamical modelling code (van den Bosch et al. 2008) written in Python.We generate a general orbit library and an additional box library based on the orbit sampling method described in (Yang et al. 2020;van den Bosch et al. 2008), and adopt a sampling number of each orbit library (  ×   ×   for the general orbit library and   ×   ×   for the box library) of 45 × 15 × 13 with a dithering of 3 in every dimension to fit the extended stellar kinematics.
The stellar kinematics have two apertures: the SAURON aperture has a higher resolution within a smaller field-of-view that is completely covered by the Mitchell one.To prevent the fitting from being biased toward the central region, we adjust the weights of each aperture in two steps: (1) we reduce the total weights of the SAURON aperture by a factor of  bins,S / bins,M , where  bins,S and  bins,M are the number of bins for the SAURON and Mitchell apertures, respectively; (2) we increase the weights of the outer regions of the Mitchell aperture to make sure that their weights are comparable to the central region.

Model of misaligned gas kinematics
For each galaxy in the subsample, we measure the position angle ( 1 ) of the Mitchell stellar velocity field, and the position angle ( 2 ) and inclination ( 2 ) of the H i velocity field using the package KINEMETRY (Krajnović et al. 2006), as shown in Table 2.There is a significant misalignment between the stellar disc (or the equatorial plane of the stellar component) and the H i disc for NGC 4203, and a clear misalignment for NGC 2685 and NGC 6798 (about 15°ignoring counter-rotation).The inclination of the stellar disc  1 is actually the viewing angle  introduced in Section 3.1, which is bound to the stellar minor-to-major axis .Considering that  is a free parameter in the model, it is possible that there is a misalignment even for NGC 2685 when taking the misalignment in the inclinations into account.Therefore, we build a misaligned gas model to describe the motions of these H i discs for all galaxies in the subsample.
As the gas discs still show nearly regular rotation in their velocity fields, we assume that a particle on the gas disc circulates around its own short axis with a circular velocity   , and it reaches nearly dynamical equilibrium in the radial direction of the stellar disc.We project   along the tangential direction of the stellar disc, where  0 is its distance to the centre of the gas disc,  0 is the distance to the short axis of the stellar disc, and  is the angle between the axes of the stellar and gas disc. is connected to the position angle and inclination of the stellar disc ( 1 ,  1 ) and those of the gas disc The ± sign takes into account the rotation direction of the inclination.
In the most simplified case where two discs have the same position angle and inclination (, ), the + sign represents the two discs fully aligned such that the model reduces to the H i model described in Yang et al. (2020), and the − sign represents the rare case where two discs rotate an angle of  in opposite directions from the line of sight, so that  = 2 × .
The equation for the dynamical equilibrium is where AccR  0 is the gravitational acceleration along the radial direction of the stellar disc at the distance  0 calculated from the gravitational potential.Given the mass distribution of a galaxy and the position on the gas disc, we can deduce the circular velocity   of the gas disc.
The modelled light-of-sight velocity of the gas disc is given by where  ′ is the azimuthal angle from the observed major axis of the gas disc.The modelled velocity is then convolved with the beam of the observation for a direct comparison with the observational velocity.
We will discuss this model of the misaligned gas kinematics further in Yang et al. in preparation.

Confidence level
We have a total of six free parameters.There are two parameters for the stellar component sampled on linear grids: (i)  is the minor-to-major axis ratio,  ∈ [0.06,  max ] with a step size of 0.02, where  max is the maximum intrinsic flattening allowed by the observed axis ratio of the surface brightness for each galaxy; (ii)  is the intermediate-to-major axis ratio,  ∈ [0.90, 0.99] with a step size of 0.01 for a nearly oblate intrinsic shape with some triaxiality still allowed. is further released to lower values in the case that the best-fitting value hits the boundary.
There are three parameters for the dark matter profile: (i)  s is the central density sampled on a logarithmic grid, log[ s /(M ⊙ • pc −3 )] ∈ [−5.00, −1.00] with a step of 0.25; (ii)  s is the scale radius sampled on a logarithmic grid, log( s /pc) ∈ [3.75, 5.50] with a step of 0.25; (iii)  is the inner slope,  ∈ [0.0, 2.0] on a linear grid with a step of 0.1; The last is the total mass scale factor sampled on a linear grid: (i)  represents the total mass produced by different IMFs,  ∈ [0.4,1.6] with a step of 0.1.
We explore the parameter grid by first creating the initial model from the given initial parameters and then computing the models adjacent to the initial model in the parameter grid.We select some models with the lowest  2 and calculate new models next to them, repeating this step until no more new models are added.We finally choose the model with the least  2 in the kinematic residual maps as the best-fitting model, with The fluctuation in the observational data results in a fluctuation in  2 .Considering that we include weights in the stellar kinematics, we adopt the bootstrapping method to measure this fluctuation for stellar kinematics.We add Gaussian perturbations whose standard deviations equal the observational 1 −  uncertainties to the observational stellar kinematics and measure its corresponding  2 bts .We repeat this process and take Δ 2 star =  2 bts + (  2 bts ) −  2 star,min as the 1 −  confidence level for stellar kinematics, where  2 bts and (  2 bts ) are the mean and standard deviation of  2 bts and  2 star,min is the minimum of  2 star .We adopt Δ 2 gas ∼ √︃ 4 2 gas,min − 2 gas as the 1 −  confidence level for gas kinematics (see Appendix B), where  2 gas,min is the minimum of  2 gas and  gas is the number of gas data points.In the general case where the model is a simplification of the data, Δ 2 gas ∼ 2 √︃  2 gas,min , consistent with the bootstrapping result in Yang et al. (2020).Therefore, the 1 −  confidence level is We show the best-fitting model with the combined stellar and H i kinematics of NGC 4203 in Figure 3 as an example of our fit, in which the major structures of the stellar kinematics and H i velocity fields are recovered.

Mass budget
We generate the enclosed mass profiles of the stellar, dark matter, and total mass from the model parameters for all sample galaxies to exhibit their mass distribution.For each galaxy, we first select all models within its 1 −  confidence level.We then obtain the enclosed mass profiles of stellar, dark matter and total mass for every model.We finally adopt the mean profile of all models within the 1 −  confidence level and take the standard deviation as the 1 −  uncertainties for the enclosed mass profiles of stellar, dark matter and total mass, respectively.
We obtain the dark matter fraction profiles, defined as the ratio of the enclosed dark matter and total mass within a certain radius, for all sample galaxies in a similar way.We first calculate the dark matter fraction profiles for every model within the 1 −  confidence level of each galaxy.We then adopt the mean dark matter fraction profile, and take the standard deviation as the 1 −  uncertainties.
Figure 4 shows the enclosed mass profile and dark matter fraction profile of NGC 4203 from two sets of models: one with combined stellar and H i kinematics and the other with stellar kinematics only.The figure extends to 5 e , and the coverage of stellar kinematics is marked with a red dashed line.The total enclosed mass profiles from the two models are consistent within the 1 −  uncertainties within the coverage of stellar kinematics: for NGC 4203, the model with stellar kinematics reports a total mass of 2.2 +0.6 −0.5 × 10 10  ⊙ while the model with combined stellar and H i kinematics reports a total mass of 2.8 +0.6 −0.6 × 10 10  ⊙ , and we also obtain similar consistency for the other galaxies in the subsample.However, the mass beyond is only accurately determined by the model with combined stellar and H i kinematics.For NGC 4203, the stellar mass dominates the enclosed mass within 2 e (approximately 5kpc), while dark matter dominates outside of 5 e (approximately 10 kpc).The models constrained by stellar kinematics have considerable uncertainties in the dark matter mass, especially when extrapolating to the outer regions beyond the data coverage, even though both SAURON and Mitchell kinematics were used.

RESULTS
We list the best-fitting parameters using the two-aperture kinematics for all sample galaxies in Table 3.We also include the results of combining the stellar kinematics and the H i misaligned discs for the subsample in additional columns for comparison and provide the orientation of the H i discs calculated using Equation 3 in Table 4.

Stellar mass
The total stellar masses of the sample galaxies are affected by the choices of the IMFs used in the stellar population models, and we have introduced a scaling factor  to indicate this mass scale. ranges from ∼ 0.5 to ∼ 1.0, suggesting a wide IMF distribution from Kroupa-like IMFs to Salpeter-like IMFs, and is not affected much by including the H i kinematics.Previous research has found that there is an IMF variation in early-type galaxies, which correlates with factors such as the stellar mass-to-light ratio (Cappellari et al. 2012) and velocity dispersion (Li et al. 2017;Zhou et al. 2019).However, we find no clear correlation between the preference of IMFs and other stellar or kinematic features within our sample, which is probably the result of a limited sample size.
Taking the scaling factor  into consideration, our full sample has a mean -band stellar mass-to-light ratio of 3.5 ⊙ / ⊙ at the galaxy centre with a root mean square (rms) scatter of 1.2 ⊙ / ⊙ , and 2.4 ⊙ / ⊙ at the galaxy outskirts with an rms scatter of 1.0 ⊙ / ⊙ .

Dark matter fraction
We obtain the dark matter fraction within 1 e , 2 e , and 5 e (  DM ( e ),  DM (2 e ), and  DM (5 e )) from the dark matter fraction profiles derived in Section 3.5, for all sample galaxies modelled with two-aperture stellar kinematics and for galaxies in the subsample modelled with combined stellar and H i kinematics, as listed in Table 5.The full sample modelled with stellar kinematics has a mean dark matter fraction of 0.19 ± 0.03 with an rms scatter of 0.10 within 1 e , 0.35 ± 0.04 with an rms scatter of 0.15 within 2 e , and 0.62 ± 0.05 with an rms scatter of 0.18 within 5 e .The subsample modelled with combined stellar and H i kinematics has a mean dark matter fraction of 0.10 ± 0.02 with an rms scatter of 0.03 within 1 e , 0.21 ± 0.03 with an rms scatter of 0.06 within 2 e , and 0.47 ± 0.03 with an rms scatter of 0.13 within 5 e , which is lower than the full sample on average but consistent with it within a 1 −  scatter.
We plot the relation between the dark matter fractions and stellar mass for the sample galaxies in the top panel of Figure 5.The results with two-aperture stellar kinematics are shown in empty red circles, whereas the results with stellar and H i kinematics are shown in full blue squares, connecting to their counterpart results in black dashed lines.Among the four galaxies in the subsample, the dark matter fraction of NGC 2685 remains unchanged throughout all radii when include H i kinematics, while the dark matter fraction of NGC 4203 increases, and the dark matter fractions of NGC 5582 and NGC 6798 decrease, especially at  DM (5 e ).This is due to the coverage of stellar kinematics.As Figure 4 shows, the total enclosed mass modelled with stellar kinematics is consistent with that modelled including H i kinematics within the coverage of stellar kinematics, but the extrapolation to large radii is biased.The disparity in the fractions of stellar and dark matter further enlarges the bias in the dark matter fraction.Therefore, the dark matter fractions modelled with the stellar kinematics of a single galaxy can be randomly biased at large radii beyond the scope of the data, leading to unchanged, higher, or lower dark matter fractions.Nevertheless, the analysis of a large sample of galaxies can still provide some observational constraints on the dark matter fractions.
For comparison, we obtain analogous relations for the early-type galaxies in the TNG-100 and EAGLE simulations.We calculate its dark matter fractions, defined as the ratio of the mass of dark matter particles and that of all particles within a certain enclosed radius scaled by  e , using the same method as that published by the TNG-100 team (Lovell et al. 2018).We divide all simulated early-type galaxies within a mass range of 10 8 − 10 12  ⊙ into mass bins with a logarithmic step of 0.5 and take the median dark matter fraction of each mass bin as the value, the 16% and 84% fractiles as the 1 −  scatter, and the 0.15% and 99.85% fractiles as the 3 −  scatter.Galaxies in the EAGLE simulation have dark matter fractions systematically higher than the TNG-100 simulation, but with a much larger scatter towards low dark matter fractions.There is a discrepancy between the observed dark matter fractions and the simulation predictions.Compared to the TNG-100 simulation, approximately 30%, 40%, and 70% of the galaxies in the sample show  DM ( e ),  DM (2 e ), and  DM (5 e ) that is consistent within a 1 −  scatter, while the remaining galaxies lie below the 1 −  region, and 20% galaxies have dark matter fractions lying below the 3 −  region even at 5 e .Compared to the EAGLE simulation, still 30%, 40% of the galaxies in the sample show  DM ( e ) and  DM (2 e ) that is consistent within a 1 −  scatter.However, only 40% galaxies have  DM (5 e ) which 1 −  is consistent with the EAGLE simulation, and almost no galaxy has  DM (5 e ) lying below the 3 −  region.
Since the mass-size relation of this sample is not fully representative of the mass-size relations in the cosmological simulations, we make a similar plot using the dark matter fraction within 1kpc, 5kpc and 10kpc (  DM (1kpc),  DM (5kpc) and  DM (10kpc)) to avoid the bias caused by  e measurements, as shown in the bottom panel of Figure 5.There is a turnover around 10 10.5  ⊙ in the relations from the simulations when the radius is scaled by  e in the top panel, but this turnover is not present when the absolute radii are used in the bottom panel, except for some galaxies with extremely low dark matter fractions around 10 10.5  ⊙ in the EAGLE simulation.
The dark matter fractions observed using absolute enclosed radii also show similar discrepancies to the relations predicted by simulations, whose strength decreases with radius.The dark matter fractions within 1kpc, 5kpc, and 10kpc show similar discrepancies between observations and simulations as those within 1 e , 2 e , and 5 e .This similarity suggests that the difference of  e between the observed and simulated galaxies does not account for all the differences of enclosed dark-matter fractions between our sample galaxies and simulations.
Dark matter fractions have been estimated using dynamical models with multiple tracers in previous studies.Stellar IFUs tracing stellar kinematics are typically used to measure dark matter fractions of galaxies within 1 e , such as those of the Atlas 3D sample (Cappellari et al. 2013), the SAMI sample (Santucci et al. 2022), the early-type galaxies presented in Jin et al. ( 2020) from the SDSS-IV MaNGA, early-type galaxies of the SDSS-IV MaNGA Dynpop sample (Zhu et al. 2024) ,and the Fornax3D project (Ding et al. 2023).However, measurements can be extended to 5 e when combining stellar IFU with globular clusters (GCs), such as in the case of M 87 (Li et al. 2020), or with GCs and planetary nebulae (PNe), such as in NGC 5846 (Zhu et al. 2016) and M 31 (Yuan et al. 2022), or with H i kinematics, such as in the case of NGC 2974 (Yang et al. 2020).GCs alone can also be used to measure dark matter fractions beyond 5 e , as seen in the SLUGGS survey (Alabi et al. 2017).We plot the correlation of dark matter fraction and stellar mass inferred from these measurements together with our results of this study, and compared it with analogous trends in the TNG-100 simulation and the EAGLE simulation, as shown in Figure 6.
The observational dark matter fractions show a large scatter, and the  DM ( e ) of Ding et al. (2023) and this work are higher than the median of the results obtained from the central stellar IFU data (Atlas 3D , SAMI and SDSS-IV MaNGA).All the observational data together still follow the same pattern as the simulations, with the dark matter fraction decreasing and then increasing as the stellar mass increases and reaching its minimum at a stellar mass of approximately 10 10.5  ⊙ .However, there is still an offset between the median dark matter fractions of the observational data and the relation of the simulations.This offset decreases as the enclosed radius increases: Various studies have measured  DM ( e ), most of which have mean values below the 3 −  scatter of the relations predicted by the simulations.The sample of galaxies with  DM (2 − 5  e ) in the literature is relatively small, and the galaxy fractions with dark matter fractions that are 1 −  consistent with simulation predictions are comparable to our sample.There are about 10% of 47 galaxies with  DM (5 e ) lying below the 3 −  region compared to the TNG-100 or the EAGLE simulations.

Dark matter halo properties
The inner slopes of the dark matter profiles provide information to constrain the properties of dark matter.Our full sample modelled with stellar kinematics has uniform NFW-like dark matter cusps, which have a mean value of 1.00 ± 0.04 with an rms scatter of 0.15.Including H i kinematics does little change to the result: The subsample modelled with combined stellar and H i kinematics has a mean value of 1.12 ± 0.03 with an rms scatter of 0.15. Figure 7 displays the probability density distribution of the inner slopes of the dark matter profiles for our sample, which represents the normalised histogram of the inner slopes of all models within the 1− confidence level for each galaxy.Although the NFW-like dark matter cusps are favoured by our full sample modelled with stellar kinematics, there is still a possibility of dark matter cores for certain galaxies.
The subsample modelled with combined stellar and H i kinematics provides stronger constraints on the inner slopes of the dark matter profiles and largely eliminates other profiles except for the NFW-like cusps.Both models rule out steeper cusps with inner slopes greater than 1.5.Previous measurements indicate that density cusps ( ∼ 1) have been widely discovered in massive early-type galaxies (e.g.Sonnenfeld et al. 2015;Wasserman et al. 2018), and also in dwarf spheroidals (Read et al. 2018;Hayashi et al. 2020).Cuspy haloes ( ∼ 0.5) have also been reported in brightest cluster galaxies of relaxed galaxy clusters (e.g.Sand et al. 2004;Newman et al. 2013) and nearby dwarf galaxies (e.g.Adams et al. 2014;Cooke et al. 2022), but rarely reported in intermediate-mass early-type galaxies (Yang et al. 2020).Density cores usually exist in dwarf galaxies (e.g.Simon et al. 2003;Battaglia et al. 2008;Walker & Peñarrubia 2011) and low surface brightness galaxies (e.g. de Blok et al. 2001;Kuzio de Naray et al. 2008), but also in a few massive early-type galaxies (Forestell & Gebhardt 2010;Oldham & Auger 2018).
The virial mass  200 (defined as the enclosed mass within radius  200 , where the average density within  200 is 200 times the critical density) and the concentration  (defined as the ratio of  200 and the scale radius  s ) are commonly used to describe dark matter haloes, whose relation shows the evolution of dark matter haloes.We plot the observational relation between the virial mass  200 and the concentration  of dark matter haloes with that shown in the simulations, as shown in Figure 8.The datasets are the same as those used in Figure 6, except for the Fornax3D project which adopted the  200 - relation obtained with a simulation in Planck cosmology as presented in Dutton & Macciò (2014).Therefore, we compare this relation with the  200 - relations in the TNG-100 and EAGLE simulations, as well as with the relation presented in Dutton & Macciò (2014).The .Relation of dark matter fraction and stellar mass from observations and simulations.From left to right: the dark matter fractions of the sample galaxies within 1 e , 2 e , and 5 e .Observations: this work (subsample using stellar IFU + H i , remaining galaxies using stellar IFU), the galaxies in the Fornax3D project (Ding et al. 2023), the Atlas 3D sample (Cappellari et al. 2013), the SAMI sample (Santucci et al. 2022), the early-type galaxies presented in et al.
Figure 7. Normalised probability density distribution of dark matter inner slope  for all sample galaxies (subsample using stellar IFU + H i , remaining galaxies using stellar IFU).
observational  200 - relation overlaps those in the simulations but shows a much larger scatter.

Total mass density
Observations show that the total mass density of early-type galaxies follows an isothermal profile with a density slope of  tot ∝  −2 and a small intrinsic scatter (e.g.Cappellari et al. 2015;Serra et al. 2016;Bellstedt et al. 2018), known as the 'bulge-halo' conspiracy, which is also reproduced in cosmological simulations (Wang et al. 2020).This isothermal profile is not a natural consequence in cold dark matter cosmology, and can be used to gain insight into the joint action of dark matter and stars in the galaxy formation process.Therefore, we obtain the total density profiles for our sample galaxies.
For each model within its 1 −  confidence level of a galaxy, we generate its total mass density profile from the derivative of its enclosed mass profile.We also fit the overall slope, the inner slope (between 0.1 e and  e ) and the outer slope (between  e and 4 e ) for each profile.We then adopt the mean profile of all models within the 1 −  confidence level and take their standard deviation as the 1 −  uncertainties for the total mass density profile.For the full sample modelled with stellar kinematics, we calculate the mean values and scatters of the overall, inner, and outer slopes using these slopes of all models within the 1 −  confidence level, respectively.The total density profiles of all sample galaxies follow a nearly isothermal profile that  tot ∝  −2 , and have a mean logarithmic slope  tot = 2.19 ± 0.04 with an rms scatter of 0.27.There is a small difference between the inner and outer slopes: the mean inner slope is 2.27±0.04 with an rms scatter of 0.31, and the mean outer slope is 1.84 ± 0.09 with an rms scatter of 0.24.Considering that only the galaxies in the subsample reach 4 e , we calculate the average outer slope for the subsample modelled with combined stellar and H i kinematics, 2.15 ± 0.05 with a scatter of 0.41, finding no obvious systematic bias due to data coverage.We plot the distribution of the total mass density for sample galaxies in Figure 9.
Our results are consistent with previous dynamical results using SAURON IFU and GC kinematics extended to 4 e , which reports  tot = 2.19 ± 0.03 with an rms scatter 0.11 in Cappellari et al. (2015) and  tot = 2.24 ± 0.05 in Bellstedt et al. (2018), as well as dynamical results using H i kinematics extended to 6 e in Serra et al. (2016), which reports  tot = 2.18 ± 0.03 with an rms scatter 0.11.Boardman et al. (2016) modelled NGC 3998 with the same stellar kinematics data sets but a constant mass-to-light ratio and an NFW dark matter profile and reported  DM ( e ) = 0.07 0.08 −0.07 .In this work, we measure  DM ( e ) = 0.05±0.06for NGC 3998 taking into account the mass-to-light gradient, consistent with the previous result.Since the stellar mass dominates the enclosed mass within 1 e in NGC 3998, the measurements of  DM ( e ) are expected to be close to 0, regardless of the assumptions made.

DISCUSSION
The accuracy of dark matter fraction measurements depends strongly on the kinematic tracers available.For the full sample modelled with stellar kinematics, the dark matter fractions have relatively large uncertainties due to the low signal-to-noise ratio at the outskirts of the Mitchell stellar kinematic data, while combining H i kinematic data can significantly reduce statistical errors for the subsample.There may be some systematic error caused by our misaligned gas model used to fit the H i kinematics.This misaligned gas model is a first-order disc model in which the radial motions on the disc are not taken into account, which will be discussed further in the work of Yang et al. in preparation.The assumption that the gas disc is in nearly dynamical equilibrium is also a choice in the case of lacking boundary conditions in the H i observations.Nevertheless, since this misaligned model reduces to the full aligned disc model when the angle between the stellar and H i discs is reduced to zero, we believe that this assumption would have a minimal effect on the measurements of dark matter fractions.In Section 4.2, we find that the dark matter fractions in observations are on average lower than the relations in the TNG-100 and the EAGLE simulations.Before concluding on a true discrepancy between observations and simulations, we would like to consider some impacts that may result in bias in the measurements of the dark matter fraction and then partially explain the discrepancy between the observational and simulated results.
(i) Sample bias.The galaxies in the sample shown in Figure 6 are mainly early-type galaxies with intermediate masses, which may not reflect the overall characteristics of nearby early-type galaxies in terms of sizes, morphologies, star formation rate, etc.This could lead to a sample bias such that the results of the study may not be applicable to all nearby early-type galaxies.
(ii)  e bias.Dark matter fractions are defined as the enclosed dark matter fraction within a certain radius, which is usually scaled with the half-light radii of the galaxy  e .However, galaxies in simulations are usually larger than those observed, particularly for early-type galaxies at the high-mass end (Genel et al. 2018;Zanisi et al. 2021).This bias in  e may lead to a bias in the measurements of the dark matter fraction.In this work, we measured dark matter fractions on an absolute scale, as shown in the bottom panel of Figure 5 and found a similar discrepancy between the observed and simulated galaxies, suggesting that the bias in  e is not a major factor in the measurements of the dark matter fraction.
(iii) Systematic bias of dynamical models.The errorbars provided in these dynamical measurements contain only statistical errors, while systematic errors are not included.Various assumptions used in dynamical models result in systematic bias: (a) The stellar mass-to-light ratios of galaxies tend to vary with the radius.Assuming a constant stellar mass-to-light ratio gradient instead of taking into account the actual gradient may lead to a significant underestimate of the dark matter fraction (Li et al. 2020).
(b) The choice of the formality of dark matter haloes may cause a systematic bias on the measurements of the dark matter fraction.(Cappellari et al. 2013;Zhu et al. 2016).One major contribution regards the inner slope of the dark matter profile: when the total enclosed mass of a galaxy is constrained at a certain radius, a density core with shallower inner slope will result in a lower dark matter fraction than a density cusp within the enclosed radius.Tests on mock core-like galaxies generated from simulation show that assuming an NFW dark matter halo, the dark matter fraction is overestimated by 38%, while this overestimation reduces to 18% if a generalised-NFW profile is used (Jin et al. 2019).
(c) The exact orbital anisotropies of tracers in galaxies are not known in advance.Assuming different orbital anisotropies can also lead to a systematic effect on the measurements of the dark matter fraction (Alabi et al. 2017).
We want to emphasise that our measurements have overcome most of these problems related to the systematic bias by considering stellar mass-to-light ratio gradient carefully, adopting a generalised-NFW profile for dark matter haloes, and creating a triaxial orbit-based Schwarzschild model.But dark matter fractions in the literature are measured with different assumptions and different data quality, which might lead to some systematic bias.
A fair comparison between simulations and observations could be achieved in the future by applying our method to a large sample of galaxies with uniformly high-quality data and ensuring that the galaxy samples from observations and simulations are representative of the nearby universe.

SUMMARY
In this paper, we map the mass distributions of 11 early-type galaxies with two-aperture stellar kinematics that extended to 2 − 4 e using the Schwarzschild model, and we map the mass distributions of a subsample of 4 galaxies with H i discs using dynamical models that combine two-aperture stellar and H i gas kinematics extended to 10 e .Our main findings are summarised below: • We adopt the stellar mass-to-light ratio obtained by stellar population synthesis assuming an Salpeter IMF, but allowing for a free scale parameter.The distribution of this scale parameter shows no preference for any particular IMF in our sample, suggesting a wide IMF distribution from Kroupa-like IMFs to Salpeter-like IMFs.
• We robustly obtain the enclosed dark matter fraction within ∼ 5 e for the 11 galaxies.The mean dark matter fraction within 1 e is 0.2 ± 0.1 (rms), and increases to 0.4 ± 0.2 (rms) within 2 e , and 0.6 ± 0.2 (rms) within 5 e for the full sample modelled with stellar kinematics, while for the subsample modelled with combined stellar and H i kinematics, the dark matter fractions are 0.5 ± 0.1 (rms) within 5 e .The dark matter fractions of the sample galaxies are generally lower than the predictions of the TNG-100 and EAGLE simulations within 1 e , and there are 40% and 70% sample galaxies with dark matter fractions that are 1 −  consistent with either the TNG-100 or the EAGLE simulations at 2 e and 5 e , respectively, while the remaining 60% and 30% galaxies lie below the 1− region.Combined with dark matter fraction measurements out to 5 e in the literature, there are about 10% of 47 galaxies lying below the 3 −  region of the TNG-100 or the EAGLE predictions.
• The sample galaxies show a preference for NFW-like dark matter haloes, with an average inner slope of 1.00 ± 0.04 with a root mean square (rms) scatter of 0.15, and the results change little when including H i kinematics.The dark matter virial mass  200 and concentrated  are only weakly constrained with large uncertainties, generally following the  200 −  correlations from the TNG-100 and the EAGLE simulations, and the relation presented in Dutton & Macciò (2014).
• The total mass density profiles are proportional to  −2 , consistent with the previous dynamical results with data coverage to 4 e (Cappellari et al. 2015;Serra et al. 2016;Bellstedt et al. 2018).
where  mod, is the best-fitting modelled data,  data, and  data, is the 'real' data and its uncertainty, and  is the number of independent data points.We assume the observation uncertainty  obs, =  data, , then possible observational data follow a Gaussian distribution that  obs, ∼ N ( data, ,  obs, 2 ).This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure1.Mass-size relation of all sample galaxies, with the fast-rotators marked in blue and slow-rotators marked in red.The  * here is calculated from   according to Equation 2 ofCappellari et al. (2013).The grey dots are all the galaxies in the Atlas 3D Survey(Cappellari et al. 2011).The gold and violet solid curves are the mass-size relations of the early-type galaxies in the TNG-100 and EAGLE simulations, and the corresponding dark and light shaded areas cover the 16th to 84th fractiles and 0.15th to 99.85th fractiles, respectively.

Figure 2 .
Figure 2. Top: MGE fitting of the -band image of NGC 4203, with the surface brightness contours in black and its MGE model in red.The foreground star on the top right and its saturation were masked in the fitting.Bottom: -band mass-to-light ratio gradient of NGC 4203.The blue dots and the corresponding uncertainties are the median values and the scatters of the -band mass-to-light in each radial bin of Poci et al. (2017).The green dots and the corresponding uncertainties are the -band mass-to-light ratio from Boardman et al. (2017) and linearly rescaled to match the -band data.The orange dots are the extensions of the farthest green dot to 200 arcsec with twice the uncertainty to avoid nonphysical values or overfitting.The red curve is the best-fitting result.

Figure 3 .
Figure 3. Best-fitting model with combined stellar and H i kinematics of NGC 4203, with all maps rotated so that their major axis is aligned with the x-axis.The top and middle panels: SAURON and Mitchell stellar kinematics.In each of the top and middle panels: data (top), model (middle), and relative residual (bottom; defined as data-model/error) of the surface brightness, stellar velocity, velocity dispersion, and the third and fourth orders of Gauss-Hermite moments from left to right.The missing pixels in the Mitchell stellar are generated by gridding because of the gaps in the Mitchell fibres, which do not affect the dynamical modelling.The bottom panel: data (left), model (middle), and relative residual (right; defined as data-model/error) of the H i velocity.We have masked the regions that are possibly the H i blobs outside the disc, see Appendix A.

Figure 4 .
Figure 4. Top: Enclosed mass profile of NGC 4203.The black, red, and blue curves stand for the total, stellar, and dark matter mass obtained by our dynamical models; the solid curves are from the model with combined stellar and H i kinematics, and the dotted curves are from the model with stellar kinematics.The corresponding shaded areas and typical errorbars show the 1 −  uncertainties for each curve.Bottom: Dark matter fraction as a function of radius.The solid blue curve is from the model with combined stellar and H i kinematics, and the dotted blue curve is from the model with stellar kinematics.The corresponding shaded area and typical errorbars show the 1 −  uncertainties for each line.The figure extends to 5 e , the red dashed lines mark the coverage of stellar kinematics.The H i kinematics extend to ∼ 400 arcsec as labelled in the figure.

Figure 5 .
Figure 5. Relation of dark matter fraction and stellar mass for all sample galaxies.Top panel: from left to right are the dark matter fractions of the sample galaxies within 1 e , 2 e , and 5 e .Bottom panel: from left to right are the dark matter fractions of the sample galaxies within 1kpc, 5kpc, and 10kpc.The empty red circles represent the results with two-aperture stellar kinematics, while the full blue squares represent the results with stellar and H i kinematics, connecting to their counterparts in black dashed lines.For comparison, we plot the solid gold and violet curves (and the dark and light shadows) to represent the predictions (and the 1 −  and 3 −  uncertainties) of the TNG-100 simulation and the EAGLE simulation.

Table 2 .
Position angle (clock-wise from west to its major-axis) and inclination for the subsample. 1 : position angle measured from the Mitchell stellar velocity field;  2 and  2 : position angle and inclination measured from the H i velocity field, and we show  2 in the form of cos  2 as it is the direct output of KINEMETRY.

Table 3 .
Best-fitting parameters of the sample galaxies. indicates the stellar mass scale factor;  is the intrinsic flattening of the stellar components;   ,   , and  are the central density, the scale radius, and the inner slope of the dark matter halo profile.The asterisks represent the parameters obtained with combined stellar and H i misaligned discs for the sub-sample.

Table 5 .
Dark matter fractions of all sample galaxies.The asterisks represent the dark matter fractions within 5 e and 5kpc obtained with combined stellar and H i misaligned discs for the sub-sample.
The residual distribution for the possible observation data is obtained by  confidence interval due to observation uncertainties.