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

HAL
open science

# Dark matter measurements combining stellar and Hi kinematics: 30 per cent 1$\sigma$ outliers with low dark matter content at 5 R<sub>e</sub>

Meng Yang, Ling Zhu, Yu Lei, Nicholas Boardman, Anne-Marie Weijmans, Raffaella Morganti, Tom Oosterloo, Pierre-Alain Duc

► To cite this version:

Meng Yang, Ling Zhu, Yu Lei, Nicholas Boardman, Anne-Marie Weijmans, et al. Dark matter measurements combining stellar and Hi kinematics: 30 per cent 1$\sigma$ outliers with low dark matter content at 5 R<sub>e</sub>. Monthly Notices of the Royal Astronomical Society, 2024, 528 (3), pp.5295-5308. 10.1093/mnras/stae335. hal-04489978

HAL Id: hal-04489978

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

Submitted on 13 Mar 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Dark matter measurements combining stellar and H I kinematics: 30 per cent 1$\sigma$ outliers with low dark matter content at $5R_e$

Meng Yang<sup>1,★</sup> Ling Zhu<sup>1,★</sup> Yu Lei<sup>1</sup> Nicholas Boardman,<sup>2</sup> Anne-Marie Weijmans,<sup>2</sup> Raffaella Morganti<sup>1,3,4</sup> Tom Oosterloo<sup>3,4</sup> and Pierre-Alain Duc<sup>5</sup>

Accepted 2024 January 29. Received 2024 January 29; in original form 2023 October 23

## ABSTRACT

We construct the Schwarzschild dynamical models for 11 early-type galaxies with the SAURON and Mitchell stellar IFUs out to $2–4R_e$, and construct dynamical models with combined stellar and H I kinematics for a subsample of four galaxies with H I velocity fields out to $10R_e$ obtained from the Westerbork Synthesis Radio Telescope, thus robustly obtaining the dark matter content out to large radii for these galaxies. Adopting a generalized-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 \pm 0.04$ (rms scatter 0.15). The mean dark matter fraction for the sample is 0.2 within $1R_e$, and increases to 0.4 at $2R_e$, and 0.6 at $5R_e$. The dark matter fractions within $1R_e$ of these galaxies are systematically lower than the predictions of both the TNG-100 and EAGLE simulations. For the dark matter fractions within $2R_e$ and $5R_e$, 40 and 70 per cent galaxies are $1\sigma$ consistent with either the TNG-100 or the EAGLE predictions, while the remaining 60 and 30 per cent galaxies lie below the $1\sigma$ region. Combined with 36 galaxies with dark matter fractions measured out to $5R_e$ in the literature, about 10 per cent of these 47 galaxies lie below the $3\sigma$ region of the TNG-100 or EAGLE predictions.

**Key words:** galaxies: haloes – galaxies: kinematics and dynamics – galaxies: structure – dark matter.

## 1 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 haloes. The distribution of dark matter within galaxies is affected not only by the characteristics of dark matter particles (e.g. Kaplinghat, Tulin & Yu 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, Springel & Hernquist 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, Eke & Frenk 1996a; Read & Gilmore 2005; Cole, Dehnen & Wilkinson 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, Walker & Steger 2019). However, for massive galaxies, baryonic feedback might not be enough to

significantly modify the distribution of dark matter (Tollet et al. 2016; Bullock & Boylan-Kolchin 2017; Hopkins et al. 2018).

Various techniques are employed to investigate dark matter in galaxies, such as strong lensing (e.g. Bolton et al. 2008; Barnabè et al. 2012; Oguri, Rusu & Falco 2014; Leier et al. 2022), hot interstellar X-ray gas (e.g. Fukazawa et al. 2006; Harris et al. 2019, 2020), and dynamical methods. Rotation curves have been widely used to late-type galaxies in the local universe (e.g. Oh et al. 2008; Martinsson et al. 2013; Oh et al. 2015; Relatores et al. 2019; Tortora et al. 2019) and at higher redshift $z \sim 1–2$ (e.g. Tiley et al. 2019; Genzel et al. 2020; Sharma, Salucci & van de Ven 2021, 2022; Nestor Shachar et al. 2023), but their applications in early-type galaxies are limited. Dynamical modelling is a key method for accurately determining the dark matter distributions in early-type galaxies with complicated dynamical structures: spherical Jeans models (Binney & Mamon 1982; Tonry 1983) suffer from a mass-anisotropy degeneracy, which requires further assumptions to infer the mass distribution; anisotropic Jeans models (Cappellari 2008) reduce ad-hoc assumptions of underlining structures for axisymmetric galaxies; particle-based made-to-measure models (Syer & Tremaine 1996) and orbit-based Schwarzschild models (Schwarzschild 1979) sidestep this degeneracy and succeed in modelling triaxial galaxies.

Most studies that assess the dark matter mass fraction of large samples with dynamical modelling rely on stellar kinematic data

<footnote n="1"><sup>1</sup>Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China</footnote>

<footnote n="2"><sup>2</sup>School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews, KY16 9SS, UK</footnote>

<footnote n="3"><sup>3</sup>Netherlands Institute for Radio Astronomy (ASTRON), Postbus 2, NL-7990 AA Dwingeloo, the Netherlands</footnote>

<footnote n="4"><sup>4</sup>Kapteyn Astronomical Institute, University of Groningen, PO Box 800, NL-9700 AV Groningen, the Netherlands</footnote>

<footnote n="5"><sup>5</sup>Observatoire astronomique de Strasbourg, 11 rue de l'université, F-67000 Strasbourg, France</footnote>

<footnote>*E-mail: myang@shao.ac.cn (MY); lzhu@shao.ac.cn (LZ)</footnote>
<table><caption>**Table 1.** Basic information of sample galaxies. All values are taken from Cappellari et al. (2011), except that the K-band stellar mass $M_{*,K}$ is estimated from $M_K$ according to equation (2) in Cappellari (2013), and $R_e$ in kpc is calculated from $R_e$ in arcsec and distance. FR/SR stands for fast-rotators/slow-rotators as classified in Cappellari et al. (2011).</caption><thead><tr><th>Galaxy</th><th>RA(°)</th><th>DEC(°)</th><th>FR/SR</th><th>Distance (Mpc)</th><th>M<sub>K</sub></th><th>log (M<sub>*,K</sub>/M<sub>☉</sub>)</th><th>R<sub>e</sub> (arcsec)</th><th>R<sub>e</sub> (kpc)</th></tr></thead><tbody><tr><td>NGC0680</td><td>27.447035</td><td>21.970827</td><td>FR</td><td>37.5</td><td>-24.17</td><td>11.09</td><td>14.5</td><td>2.64</td></tr><tr><td>NGC1023</td><td>40.100052</td><td>39.063251</td><td>FR</td><td>11.1</td><td>-24.01</td><td>11.02</td><td>47.9</td><td>2.58</td></tr><tr><td>NGC2685</td><td>133.894791</td><td>58.734409</td><td>FR</td><td>16.7</td><td>-22.78</td><td>10.48</td><td>25.7</td><td>2.08</td></tr><tr><td>NGC2764</td><td>137.072983</td><td>21.443447</td><td>FR</td><td>39.6</td><td>-23.19</td><td>10.66</td><td>12.3</td><td>2.36</td></tr><tr><td>NGC3522</td><td>166.668549</td><td>20.085621</td><td>SR</td><td>25.5</td><td>-21.67</td><td>9.99</td><td>10.2</td><td>1.26</td></tr><tr><td>NGC3998</td><td>179.484039</td><td>55.453564</td><td>FR</td><td>13.7</td><td>-23.33</td><td>10.73</td><td>20.0</td><td>1.33</td></tr><tr><td>NGC4203</td><td>183.770935</td><td>33.197243</td><td>FR</td><td>14.7</td><td>-23.44</td><td>10.77</td><td>29.5</td><td>2.1</td></tr><tr><td>NGC5582</td><td>215.179703</td><td>39.693584</td><td>FR</td><td>27.7</td><td>-23.28</td><td>10.7</td><td>27.5</td><td>3.69</td></tr><tr><td>NGC5631</td><td>216.638687</td><td>56.582664</td><td>SR</td><td>27.0</td><td>-23.70</td><td>10.89</td><td>20.9</td><td>2.74</td></tr><tr><td>NGC6798</td><td>291.013306</td><td>53.624752</td><td>FR</td><td>37.5</td><td>-23.52</td><td>10.8</td><td>16.9</td><td>3.07</td></tr><tr><td>UGC03960</td><td>115.094856</td><td>23.275089</td><td>SR</td><td>33.2</td><td>-21.89</td><td>10.09</td><td>17.4</td><td>2.8</td></tr></tbody></table>

from spatially resolved integral field unit (IFU) spectrographs, such as Atlas³D (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 effective radii ($R_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 haloes (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 (Weijmans 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 (Côté et al. 2001, 2003; Romanowsky et al. 2003; Deason et al. 2012; Alabi et al. 2016, 2017; Dumont et al. 2023), or in combination with stellar kinematics (de Lorenzi et al. 2008, 2009; Napolitano et al. 2009; Das et al. 2011; Napolitano et al. 2011; Morganti et al. 2013; Napolitano et al. 2014; Zhu et al. 2014, 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 H1 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, McGaugh & Schombert 2016). Combining H1 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 $5R_e$; (2) Observations that reach out to $5R_e$ are inadequate and display a considerable scatter, with few measurements in galaxies with stellar masses below $10^{10} M_{☉}$; (3) Dark matter fractions measured using stellar kinematics only covering 1–2$R_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 summarizes our work.

## 2 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} M_{☉}$ and a range in effective radius of 1–4 kpc, and includes three slow-rotators dominated by random motions and eight fast-rotators dominated by circular motions, as classified in Cappellari et al. (2011). We show the mass-size relation of the sample in Fig. 1, and compare it with the full Atlas³D sample. We also plot the mass-size relations of the early-type galaxies in the TNG-100 simulation (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018) and the EAGLE simulation (Crain et al. 2015; Schaye et al. 2015), in which the $R_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^{-1}$ 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 $R_e$ to measure enclosed dark matter fractions.

We used the r-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 field-of-view of 1 deg², and these images have a limiting magnitude of ~29 mag arcsec⁻² in the r-band, which allows us to trace the surface brightness out to $5R_e$, beyond the range of our kinematic data sets. As the MATLAS r-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 r-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
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Mass-size relation of all sample galaxies, with the fast-rotators marked in blue and slow-rotators marked in red. The $M_*$ here is calculated from $M_K$ according to equation (2) of 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.</figcaption></figure>

images to reach a limiting magnitude of 23.2 in the *r*-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² with lenslets of 0.94 arcsec² and an instrumental velocity dispersion of 105 km s⁻<fnref n="1" />, and covers the central 1R⊙ of the sample galaxies. The spectra were Voronoi binned (Cappellari & Copin 2003) to a signal-to-noise ratio of 40, and fitted with the penalized pixel fitting method (pPXF; Cappellari & Emsellem 2004) to obtain the stellar kinematics (velocity, velocity dispersion, and Gaussian-Hermite moments $h_3, h_4$).<fnref n="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², a large fibre size of 2.08 arcsec and a maximum instrumental velocity resolution of 42 km s⁻<fnref n="1" />, and covers 2–4R⊙ 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 h for each galaxy. The raw data were processed and calibrated using the MIRIAD software package (Sault, Teuben & Wright 1995), resulting in a datacube with a spectral resolution of 16 km s⁻<fnref n="1" />. 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.

## 3 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 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.

### 3.1 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 *r*-band image, we adopt the MATLAS *r*-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 *r*-band image; we then replace the saturated central region in the MATLAS *r*-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 two-dimensional Gaussian profiles (MGE; Emsellem, Monnet & Bacon 1994; Cappellari 2002), and show the MGE fitting of NGC 4203 in the top panel of Fig. 2.

The sample galaxies display strong gradients in their *r*-band mass-to-light ($M_*/L_r$) ratio obtained from the SAURON spectra Poci, Cappellari & McDermid (2017) and the *V*-band mass-to-light ratio ($M_*/L_v$) obtained from the Mitchell spectra (Boardman et al. 2017), which are obtained with the stellar population synthesis models. Therefore, we take the mass-to-light ratio gradient into account instead of assuming a constant mass-to-light ratio (Yang et al. 2020). We first rescale the Mitchell $M_*/L_v$ linearly to match the SAURON $M_*/L_r$. 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-hand panel of Fig. 2.

The stellar mass and the corresponding mass-to-light ratio obtained 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 $\alpha$ to indicate the total mass scale due to the IMF: $\alpha = 1.0$ is set to the stellar mass produced by the

<footnote n="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</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: MGE fitting of the *r*-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: *r*-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 *r*-band mass-to-light in each radial bin of Poci, Cappellari & McDermid (2017). The green dots and the corresponding uncertainties are the *V*-band mass-to-light ratio from Boardman et al. (2017) and linearly rescaled to match the *r*-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.</figcaption></figure>

Salpeter IMF (Salpeter 1955), while the Kroupa IMF (Kroupa 2001) produces a stellar mass roughly with $\alpha = 0.6$ (Speagle et al. 2014).

We then deproject the two-dimensional MGE model to a triaxial MGE model with a set of viewing angles ($\theta$, $\psi$, $\phi$). The viewing angles are linked to the three ratios of the intrinsic shape of the galaxy ($p$, $q$, $u$) (van den Bosch et al. 2008), which are the free parameters of the gravitational potential. We set the minor-to-major axis ratio $q$ and the intermediate-to-major axis ratio $p$ as free parameters to allow triaxiality, and fix $u = 0.9999$ to reduce model freedom.

We choose a spherical generalized-NFW profile (Zhao 1996; Navarro, Frenk & White 1996b) for the dark matter halo,

$$ \rho(r) = \frac{\rho_s}{(r/r_s)^\gamma [1 + (r/r_s)^\eta]^{(3-\gamma)/\eta}} \quad (1) $$

which contains four parameters: $\rho_s$ is the scale density, $r_s$ is the scale radius, $\gamma$ is the inner slope (when $\gamma = 1$, it reduces to the NFW profile), while $\eta$ 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 M_{\odot}$ and its softening length at $10^{-3}$ pc.

In total, we have six free parameters in the dynamical model: the total mass scale $\alpha$, the minor-to-major axis ratio $q$, and the intermediate-to-major axis ratio $p$, the scale density $\rho_s$, the scale radius $r_s$, and the inner slope $\gamma$ of the dark matter halo.

## 3.2 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 (van den Bosch et al. 2008; Yang et al. 2020), and adopt a sampling number of each orbit library ($n_E \times n_\theta \times n_R$ for the general orbit library and $n_E \times n_\theta \times n_\phi$ 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 $n_{\text{bins},S}/n_{\text{bins},M}$, where $n_{\text{bins},S}$ and $n_{\text{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.

## 3.3 Model of misaligned gas kinematics

For each galaxy in the subsample, we measure the position angle ($\psi_1$) of the Mitchell stellar velocity field, and the position angle ($\psi_2$) and inclination ($i_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 $i_1$ is actually the viewing angle $\theta$ introduced in Section 3.1, which is bound to the stellar minor-to-major axis $q$. Considering that $q$ 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

<table><caption>**Table 2.** Position angle (clock-wise from west to its major-axis) and inclination for the subsample. $\psi_1$: position angle measured from the Mitchell stellar velocity field; $\psi_2$ and $i_2$: position angle and inclination measured from the H I velocity field, and we show $i_2$ in the form of $\cos i_2$ as it is the direct output of KINEMETRY.</caption><thead><tr><th>Galaxy</th><th>ψ<sub>1</sub></th><th>ψ<sub>2</sub></th><th>cos i<sub>2</sub></th></tr></thead><tbody><tr><td>NGC2685</td><td>53.5°</td><td>40.5°</td><td>0.267</td></tr><tr><td>NGC4203</td><td>-104.5°</td><td>64.0°</td><td>0.835</td></tr><tr><td>NGC5582</td><td>59.5°</td><td>58.0°</td><td>0.683</td></tr><tr><td>NGC6798</td><td>-49.0°</td><td>145.0°</td><td>0.525</td></tr></tbody></table>
its own short axis with a circular velocity $v_c$, and it reaches nearly dynamical equilibrium in the radial direction of the stellar disc. We project $v_c$ along the tangential direction of the stellar disc,

$$v_{\phi} = \frac{r_0 \cos \delta}{R_0} v_c, \qquad (2)$$

where $r_0$ is its distance to the centre of the gas disc, $R_0$ is the distance to the short axis of the stellar disc, and $\delta$ is the angle between the axes of the stellar and gas disc. $\delta$ is connected to the position angle and inclination of the stellar disc ($\psi_1, i_1$) and those of the gas disc ($\psi_2, i_2$) by

$$\cos \delta = \cos i_1 \cos i_2 \pm \sin i_1 \sin i_2 \cos (\psi_1 - \psi_2). \qquad (3)$$

The $\pm$ 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 ($\psi, i$), the $+$ sign represents the two discs fully aligned such that the model reduces to the H1 model described in Yang et al. (2020), and the "-" sign represents the rare case, where two discs rotate an angle of $i$ in opposite directions from the line of sight, so that $\delta = 2 \times i$.

The equation for the dynamical equilibrium is

$$v_{\phi} = \sqrt{-\text{AccR}_{R_0} \cdot R_0}, \qquad (4)$$

where $\text{AccR}_{R_0}$ is the gravitational acceleration along the radial direction of the stellar disc at the distance $R_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 $v_c$ of the gas disc.

The modelled light-of-sight velocity of the gas disc is given by

$$v_{\text{mod}} = v_c \sin i_2 \cos \phi', \qquad (5)$$

where $\phi'$ 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).

## 3.4 Confidence level

We have a total of six free parameters. There are two parameters for the stellar component sampled on linear grids:

(i) $q$ is the minor-to-major axis ratio, $q \in [0.06, q_{\text{max}}]$ with a step size of 0.02, where $q_{\text{max}}$ is the maximum intrinsic flattening allowed by the observed axis ratio of the surface brightness for each galaxy;

(ii) $p$ is the intermediate-to-major axis ratio, $p \in [0.90, 0.99]$ with a step size of 0.01 for a nearly oblate intrinsic shape with some triaxiality still allowed. $p$ 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) $\rho_s$ is the central density sampled on a logarithmic grid, $\log[\rho_s/(\text{M}_\odot \cdot \text{pc}^{-3})] \in [-5.00, -1.00]$ with a step of 0.25;

(ii) $r_s$ is the scale radius sampled on a logarithmic grid, $\log(r_s/\text{pc}) \in [3.75, 5.50]$ with a step of 0.25;

(iii) $\gamma$ is the inner slope, $\gamma \in [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) $\alpha$ represents the total mass produced by different IMFs, $\alpha \in [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 $\chi^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 $\chi^2$ in the kinematic residual maps as the best-fitting model, with

$$\chi^2 = \chi_{\text{star}}^2 + \chi_{\text{gas}}^2. \qquad (6)$$

The fluctuation in the observational data results in a fluctuation in $\chi^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\sigma$ uncertainties to the observational stellar kinematics and measure its corresponding $\chi^2_{\text{bts}}$. We repeat this process and take $\Delta\chi^2_{\text{star}} = \chi^2_{\text{bts}} + \sigma(\chi^2_{\text{bts}}) - \chi^2_{\text{star,min}}$ as the $1\sigma$ confidence level for stellar kinematics, where $\chi^2_{\text{bts}}$ and $\sigma(\chi^2_{\text{bts}})$ are the mean and standard deviation of $\chi^2_{\text{bts}}$ and $\chi^2_{\text{star,min}}$ is the minimum of $\chi^2_{\text{star}}$. We adopt $\Delta\chi^2_{\text{gas}} \sim \sqrt{4\chi^2_{\text{gas,min}} - 2N_{\text{gas}}}$ as the $1\sigma$ confidence level for gas kinematics (see Appendix B), where $\chi^2_{\text{gas,min}}$ is the minimum of $\chi^2_{\text{gas}}$ and $N_{\text{gas}}$ is the number of gas data points. In the general case, where the model is a simplification of the data, $\Delta\chi^2_{\text{gas}} \sim 2\sqrt{\chi^2_{\text{gas,min}}}$, consistent with the bootstrapping result in Yang et al. (2020). Therefore, the $1\sigma$ confidence level is

$$\Delta\chi^2 = \Delta\chi^2_{\text{star}} + \Delta\chi^2_{\text{gas}}. \qquad (7)$$

We show the best-fitting model with the combined stellar and H1 kinematics of NGC 4203 in Fig. 3 as an example of our fit, in which the major structures of the stellar kinematics and H1 velocity fields are recovered.

## 3.5 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\sigma$ 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\sigma$ confidence level and take the standard deviation as the $1\sigma$ 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\sigma$ confidence level of each galaxy. We then adopt the mean dark matter fraction profile, and take the standard deviation as the $1\sigma$ uncertainties.

Fig. 4 shows the enclosed mass profile and dark matter fraction profile of NGC 4203 from two sets of models: one with combined stellar and H1 kinematics and the other with stellar kinematics only. The figure extends to $5R_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\sigma$ 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} \times 10^{10} \text{M}_\odot$, while the model with combined stellar and H1 kinematics reports a total mass of $2.8^{+0.6}_{-0.6} \times 10^{10} \text{M}_\odot$, 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 H1
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Best-fitting model with combined stellar and H1 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 H1 velocity. We have masked the regions that are possibly the H1 blobs outside the disc, see Appendix A.</figcaption></figure>

kinematics. For NGC 4203, the stellar mass dominates the enclosed mass within $2R_e$ (approximately 5 kpc), while dark matter dominates outside of $5R_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.

## 4 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 H1 misaligned discs for the subsample in additional columns for comparison and provide the orientation of the H1 discs calculated using equation (3) in Table 4.
<figure><img src="image_6.png" /><figcaption>**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\sigma$ 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\sigma$ uncertainties for each line. The figure extends to $5R_\odot$, the red dashed lines mark the coverage of stellar kinematics. The H I kinematics extend to ~400 arcsec as labelled in the figure.</figcaption></figure>

<table><caption>**Table 4.** Orientations of misaligned H I discs. The projection direction represents $\pm$ in equation (3), and $\delta$ represents the angle between the rotating axes of the stellar and H I discs.</caption><thead><tr><th>Galaxy</th><th>Projection direction</th><th>δ</th></tr></thead><tbody><tr><td>NGC2685</td><td>–</td><td>151.4°</td></tr><tr><td>NGC4203</td><td>–</td><td>21.2°</td></tr><tr><td>NGC5582</td><td>–</td><td>116.1°</td></tr><tr><td>NGC6798</td><td>+</td><td>123.0°</td></tr></tbody></table>

## 4.1 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 $\alpha$ to indicate this mass scale. $\alpha$ 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 $\alpha$ into consideration, our full sample has a mean r-band stellar mass-to-light ratio of $3.5 M_\odot/L_\odot$ at the galaxy centre with a root mean square (rms) scatter of $1.2 M_\odot/L_\odot$, and $2.4 M_\odot/L_\odot$ at the galaxy outskirts with an rms scatter of $1.0 M_\odot/L_\odot$.

## 4.2 Dark matter fraction

We obtain the dark matter fraction within $1R_\odot$, $2R_\odot$, and $5R_\odot$ ($f_{DM}(R_\odot)$, $f_{DM}(2R_\odot)$, and $f_{DM}(5R_\odot)$) 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 \pm 0.03$ with an rms scatter of $0.10$ within $1R_\odot$, $0.35 \pm 0.04$ with an rms scatter of $0.15$ within $2R_\odot$, and $0.62 \pm 0.05$ with an rms scatter of $0.18$ within $5R_\odot$. The subsample modelled with combined stellar and H I kinematics has a mean dark matter fraction

<table><caption>**Table 3.** Best-fitting parameters of the sample galaxies. $\alpha$ indicates the stellar mass scale factor; $q$ is the intrinsic flattening of the stellar components; $\rho_s$, $r_s$, and $\gamma$ 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 subsample.</caption><thead><tr><th>Galaxy</th><th>α</th><th>α*</th><th>q</th><th>q*</th><th>ρ<sub>s</sub> (10<sup>-3</sup> M<sub>☉</sub> pc<sup>-3</sup>)</th><th>ρ<sub>s</sub><sup>*</sup></th><th>r<sub>s</sub> (kpc)</th><th>r<sub>s</sub><sup>*</sup></th><th>γ</th><th>γ*</th></tr></thead><tbody><tr><td>NGC0680</td><td>0.50 ± 0.07</td><td>–</td><td>0.66 ± 0.08</td><td>–</td><td>8.3<sup>+28.5</sup><sub>-6.4</sub></td><td>–</td><td>30.4<sup>+64.5</sup><sub>-20.7</sub></td><td>–</td><td>0.85 ± 0.33</td><td>–</td></tr><tr><td>NGC1023</td><td>0.60 ± 0.03</td><td>–</td><td>0.20 ± 0.05</td><td>–</td><td>1.6<sup>+1.1</sup><sub>-0.6</sub></td><td>–</td><td>87.9<sup>+83.0</sup><sub>-42.7</sub></td><td>–</td><td>0.97 ± 0.15</td><td>–</td></tr><tr><td>NGC2685</td><td>0.81 ± 0.03</td><td>0.81 ± 0.03</td><td>0.22 ± 0.03</td><td>0.23 ± 0.03</td><td>0.4<sup>+0.3</sup><sub>-0.2</sub></td><td>0.4<sup>+0.3</sup><sub>-0.2</sub></td><td>85.8<sup>+83.5</sup><sub>-42.3</sub></td><td>78.1<sup>+70.4</sup><sub>-37.0</sub></td><td>1.13 ± 0.12</td><td>1.14 ± 0.12</td></tr><tr><td>NGC2764</td><td>1.19 ± 0.20</td><td>–</td><td>0.10 ± 0.05</td><td>–</td><td>9.1<sup>+8.5</sup><sub>-4.4</sub></td><td>–</td><td>12.1<sup>+8.2</sup><sub>-4.9</sub></td><td>–</td><td>0.97 ± 0.22</td><td>–</td></tr><tr><td>NGC3522</td><td>0.59 ± 0.11</td><td>–</td><td>0.39 ± 0.10</td><td>–</td><td>3.4<sup>+27.8</sup><sub>-3.0</sub></td><td>–</td><td>17.1<sup>+89.1</sup><sub>-14.4</sub></td><td>–</td><td>0.94 ± 0.29</td><td>–</td></tr><tr><td>NGC3998</td><td>1.04 ± 0.10</td><td>–</td><td>0.31 ± 0.08</td><td>–</td><td>6.2<sup>+23.4</sup><sub>-4.9</sub></td><td>–</td><td>9.9<sup>+7.8</sup><sub>-4.4</sub></td><td>–</td><td>1.01 ± 0.05</td><td>–</td></tr><tr><td>NGC4203</td><td>0.47 ± 0.05</td><td>0.46 ± 0.05</td><td>0.63 ± 0.10</td><td>0.62 ± 0.09</td><td>0.8<sup>+4.3</sup><sub>-0.7</sub></td><td>4.0<sup>+9.5</sup><sub>-2.8</sub></td><td>14.9<sup>+23.9</sup><sub>-9.2</sub></td><td>21.0<sup>+26.7</sup><sub>-11.8</sub></td><td>0.92 ± 0.44</td><td>0.86 ± 0.35</td></tr><tr><td>NGC5582</td><td>0.78 ± 0.06</td><td>0.86 ± 0.05</td><td>0.47 ± 0.06</td><td>0.56 ± 0.01</td><td>1.5<sup>+0.6</sup><sub>-0.4</sub></td><td>1.5<sup>+0.1</sup><sub>-0.1</sub></td><td>25.4<sup>+8.6</sup><sub>-6.4</sub></td><td>17.8<sup>+13.8</sup><sub>-7.8</sub></td><td>1.36 ± 0.12</td><td>1.24 ± 0.20</td></tr><tr><td>NGC5631</td><td>0.72 ± 0.12</td><td>–</td><td>0.54 ± 0.15</td><td>–</td><td>11.7<sup>+24.8</sup><sub>-7.9</sub></td><td>–</td><td>16.2<sup>+26.3</sup><sub>-10.0</sub></td><td>–</td><td>0.76 ± 0.45</td><td>–</td></tr><tr><td>NGC6798</td><td>0.64 ± 0.06</td><td>0.90 ± 0.05</td><td>0.20 ± 0.05</td><td>0.34 ± 0.02</td><td>0.9<sup>+0.9</sup><sub>-0.5</sub></td><td>0.8<sup>+0.2</sup><sub>-0.2</sub></td><td>90.0<sup>+72.5</sup><sub>-40.2</sub></td><td>17.8<sup>+13.8</sup><sub>-7.8</sub></td><td>1.07 ± 0.12</td><td>1.24 ± 0.15</td></tr><tr><td>UGC03960</td><td>0.93 ± 0.20</td><td>–</td><td>0.67 ± 0.09</td><td>–</td><td>3.1<sup>+5.9</sup><sub>-2.0</sub></td><td>–</td><td>30.8<sup>+48.2</sup><sub>-18.8</sub></td><td>–</td><td>0.99 ± 0.16</td><td>–</td></tr></tbody></table>
<table><caption>**Table 5.** Dark matter fractions of all sample galaxies. The asterisks represent the dark matter fractions within $5R_e$ and 5 kpc obtained with combined stellar and H1 misaligned discs for the subsample.</caption><thead><tr><th>Galaxy</th><th>f<sub>DM</sub>(R<sub>e</sub>)</th><th>f<sub>DM</sub>(2R<sub>e</sub>)</th><th>f<sub>DM</sub>(5R<sub>e</sub>)</th><th>f<sup>*</sup><sub>DM</sub>(5R<sub>e</sub>)</th><th>f<sub>DM</sub>(1 kpc)</th><th>f<sub>DM</sub>(5 kpc)</th><th>f<sub>DM</sub>(10 kpc)</th><th>f<sup>*</sup><sub>DM</sub>(10 kpc)</th></tr></thead><tbody><tr><td>NGC0680</td><td>0.21 ± 0.10</td><td>0.44 ± 0.13</td><td>0.75 ± 0.10</td><td>—</td><td>0.06 ± 0.04</td><td>0.42 ± 0.12</td><td>0.67 ± 0.12</td><td>—</td></tr><tr><td>NGC1023</td><td>0.12 ± 0.04</td><td>0.29 ± 0.06</td><td>0.64 ± 0.06</td><td>—</td><td>0.04 ± 0.02</td><td>0.28 ± 0.06</td><td>0.53 ± 0.07</td><td>—</td></tr><tr><td>NGC2685</td><td>0.13 ± 0.04</td><td>0.28 ± 0.06</td><td>0.64 ± 0.06</td><td>0.62 ± 0.06</td><td>0.05 ± 0.02</td><td>0.34 ± 0.06</td><td>0.62 ± 0.06</td><td>0.60 ± 0.07</td></tr><tr><td>NGC2764</td><td>0.31 ± 0.10</td><td>0.48 ± 0.10</td><td>0.74 ± 0.06</td><td>—</td><td>0.22 ± 0.11</td><td>0.49 ± 0.09</td><td>0.70 ± 0.07</td><td>—</td></tr><tr><td>NGC3522</td><td>0.16 ± 0.13</td><td>0.30 ± 0.19</td><td>0.54 ± 0.25</td><td>—</td><td>0.12 ± 0.11</td><td>0.48 ± 0.24</td><td>0.65 ± 0.25</td><td>—</td></tr><tr><td>NGC3998</td><td>0.05 ± 0.06</td><td>0.12 ± 0.12</td><td>0.30 ± 0.23</td><td>—</td><td>0.03 ± 0.04</td><td>0.24 ± 0.19</td><td>0.40 ± 0.27</td><td>—</td></tr><tr><td>NGC4203</td><td>0.04 ± 0.05</td><td>0.09 ± 0.10</td><td>0.23 ± 0.21</td><td>0.58 ± 0.07</td><td>0.02 ± 0.03</td><td>0.11 ± 0.11</td><td>0.22 ± 0.20</td><td>0.56 ± 0.07</td></tr><tr><td>NGC5582</td><td>0.23 ± 0.08</td><td>0.39 ± 0.10</td><td>0.63 ± 0.11</td><td>0.40 ± 0.08</td><td>0.08 ± 0.03</td><td>0.29 ± 0.09</td><td>0.47 ± 0.11</td><td>0.26 ± 0.08</td></tr><tr><td>NGC5631</td><td>0.20 ± 0.15</td><td>0.39 ± 0.19</td><td>0.66 ± 0.22</td><td>—</td><td>0.07 ± 0.09</td><td>0.36 ± 0.18</td><td>0.58 ± 0.22</td><td>—</td></tr><tr><td>NGC6798</td><td>0.23 ± 0.06</td><td>0.45 ± 0.08</td><td>0.79 ± 0.05</td><td>0.29 ± 0.06</td><td>0.07 ± 0.03</td><td>0.38 ± 0.08</td><td>0.64 ± 0.07</td><td>0.19 ± 0.05</td></tr><tr><td>UGC03960</td><td>0.41 ± 0.15</td><td>0.64 ± 0.13</td><td>0.85 ± 0.08</td><td>—</td><td>0.16 ± 0.10</td><td>0.60 ± 0.14</td><td>0.78 ± 0.09</td><td>—</td></tr></tbody></table>

of $0.10 \pm 0.02$ with an rms scatter of 0.03 within $1R_e$, $0.21 \pm 0.03$ with an rms scatter of 0.06 within $2R_e$, and $0.47 \pm 0.03$ with an rms scatter of 0.13 within $5R_e$, which is lower than the full sample on average but consistent with it within a $1\sigma$ scatter.

We plot the relation between the dark matter fractions and stellar mass for the sample galaxies in the top panel of Fig. 5. The results with two-aperture stellar kinematics are shown in empty red circles, whereas the results with stellar and H1 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 H1 kinematics, while the dark matter fraction of NGC 4203 increases, and the dark matter fractions of NGC 5582 and NGC 6798 decrease, especially at $f_{DM}(5R_e)$. This is due to the coverage of stellar kinematics. As Fig. 4 shows, the total enclosed mass modelled with stellar kinematics is consistent with that modelled including H1 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 to that of all particles within a certain enclosed radius scaled by $R_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} M_\odot$ 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 per cent fractiles as the $1\sigma$ scatter, and the 0.15 and 99.85 per cent fractiles as the $3\sigma$ 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 per cent of the galaxies in the sample show $f_{DM}(R_e)$, $f_{DM}(2R_e)$, and $f_{DM}(5R_e)$ that is consistent within a $1\sigma$ scatter, while the remaining galaxies lie below the $1\sigma$ region, and 20 per cent galaxies have dark matter fractions lying below the $3\sigma$ region even at $5R_e$. Compared to the EAGLE

simulation, still 30, 40 per cent of the galaxies in the sample show $f_{DM}(R_e)$ and $f_{DM}(2R_e)$ that are consistent within a $1\sigma$ scatter. However, only 40 per cent galaxies have $f_{DM}(5R_e)$, which $1\sigma$ is consistent with the EAGLE simulation, and almost no galaxy has $f_{DM}(5R_e)$ lying below the $3\sigma$ 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 1, 5, and 10 kpc ($f_{DM}(1kpc)$, $f_{DM}(5kpc)$, and $f_{DM}(10kpc)$) to avoid the bias caused by $R_e$ measurements, as shown in the bottom panel of Fig. 5. There is a turnover around $10^{10.5} M_\odot$ in the relations from the simulations when the radius is scaled by $R_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} M_\odot$ 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 1, 5, and 10 kpc show similar discrepancies between observations and simulations as those within $1R_e$, $2R_e$, and $5R_e$. This similarity suggests that the difference of $R_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 $1R_e$, such as those of the Atlas$^3$D 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 $5R_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 H1 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 $5R_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 Fig. 6.

The observational dark matter fractions show a large scatter, and the $f_{DM}(R_e)$ of Ding et al. (2023) and this work are higher than
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Relation between 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 $1R_e$, $2R_e$, and $5R_e$. Bottom panel: from left to right are the dark matter fractions of the sample galaxies within 1, 5, and 10 kpc. 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\sigma$ and $3\sigma$ uncertainties) of the TNG-100 simulation and the EAGLE simulation.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** 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 $1R_e$, $2R_e$, and $5R_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 Atlas3D sample (Cappellari et al. 2013), the SAMI sample (Santucci et al. 2022), the early-type galaxies presented in Jin et al. (2020) from SDSS-IV MaNGA, early-type galaxies of the SDSS-IV MaNGA Dynpop sample (Zhu et al. 2024), the SLUGGS survey (Alabi et al. 2017), NGC2974 (Yang et al. 2020), NGC5846 (Zhu et al. 2016), M87 (Li et al. 2020), and M31 (Yuan et al. 2022). Simulations: TNG-100 simulation (in gold) and EAGLE simulation (in violet).</figcaption></figure>

the median of the results obtained from the central stellar IFU data (Atlas3D, 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} M_\odot$. 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 $f_{\text{DM}}(R_e)$, most of which have mean values below the $3\sigma$ scatter of the relations predicted by the simulations. The sample of galaxies with $f_{\text{DM}}(2-5 R_e)$ in the literature is relatively small, and the galaxy fractions with dark matter fractions that are $1\sigma$ consistent with simulation predictions are comparable to our sample. There are about 10 per cent of 47 galaxies with $f_{\text{DM}}(5R_e)$ lying below the $3\sigma$ region compared to the TNG-100 or the EAGLE simulations.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Normalized probability density distribution of dark matter inner slope $\gamma$ for all sample galaxies (subsample using stellar IFU + H1, remaining galaxies using stellar IFU).</figcaption></figure>

## 4.3 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 \pm 0.04$ with an rms scatter of 0.15. Including H1 kinematics does little change to the result. The subsample modelled with combined stellar and H1 kinematics has a mean value of $1.12 \pm 0.03$ with an rms scatter of 0.15. Fig. 7 displays the probability density distribution of the inner slopes of the dark matter profiles for our sample, which represents the normalized histogram of the inner slopes of all models within the $1\sigma$ 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 H1 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 ($\gamma \sim 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, Walker & Steger 2018; Hayashi, Chiba & Ishiyama 2020). Cuspy haloes ($\gamma \sim 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, McGaugh & Rubin 2001; Kuzio de Naray, McGaugh & de Blok 2008), but also in a few massive early-type galaxies (Forestell & Gebhardt 2010; Oldham & Auger 2018).

The virial mass $M_{200}$ (defined as the enclosed mass within radius $r_{200}$, where the average density within $r_{200}$ is 200 times the critical density) and the concentration $C$ (defined as the ratio of $r_{200}$ and the scale radius $r_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 $M_{200}$ and the concentration $C$ of dark matter haloes with that shown in the simulations, as shown in Fig. 8. The datasets are the same as those used in Fig. 6, except for the Fornax3D project which adopted the $M_{200}-C$ relation obtained with a simulation in Planck cosmology as presented in Dutton & Macciò (2014). Therefore, we compare this relation with the $M_{200}-C$ relations in the TNG-100 and EAGLE

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Relation of dark matter virial mass $M_{200}$ and concentration $C$ from observations and simulations. Observations: this work (subsample using stellar IFU + H1, remaining galaxies using stellar IFU), the SLUGGS survey (Alabi et al. 2017), NGC2974 (Yang et al. 2020), NGC5846 (Zhu et al. 2016), M87 (Li et al. 2020), and M31 (Yuan et al. 2022). Simulations: the TNG-100 simulation (in gold), the EAGLE simulation (in violet), and the relation presented in Dutton & Macciò (2014) (in black).</figcaption></figure>

simulations, as well as with the relation presented in Dutton & Macciò (2014). The observational $M_{200}-C$ relation overlaps those in the simulations but shows a much larger scatter.

## 4.4 Total mass density

Observations show that the total mass density of early-type galaxies follows an isothermal profile with a density slope of $\rho_{tot} \propto r^{-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\sigma$ 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.1R_e$ and $R_e$), and the outer slope (between $R_e$ and $4R_e$) for each profile. We then adopt the mean profile of all models within the $1\sigma$ confidence level and take their standard deviation as the $1\sigma$ 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\sigma$ confidence level, respectively. The total density profiles of all sample galaxies follow a nearly isothermal profile that $\rho_{tot} \propto r^{-2}$, and have a mean logarithmic slope $\gamma_{tot} = 2.19 \pm 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 \pm 0.04$ with an rms scatter of 0.31, and the mean outer slope is $1.84 \pm 0.09$ with an rms scatter of 0.24. Considering that only the galaxies in the subsample reach $4R_e$, we calculate the average outer slope for the subsample modelled with combined stellar and H1 kinematics, $2.15 \pm 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 Fig. 9.

Our results are consistent with previous dynamical results using SAURON IFU and GC kinematics extended to $4R_e$, which reports $\gamma_{tot} = 2.19 \pm 0.03$ with an rms scatter 0.11 in Cappellari et al. (2015).
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Profiles of total mass density distribution for all sample galaxies (subsample using stellar IFU + H1, remaining galaxies using stellar IFU). $\rho \propto r^{-1}$ (NFW inner slope), $\rho \propto r^{-2}$ (isothermal profile), and $\rho \propto r^{-3}$ are shown in dashed lines for comparison.</figcaption></figure>

and $\gamma_{\text{tot}} = 2.24 \pm 0.05$ in Bellstedt et al. (2018), as well as dynamical results using H1 kinematics extended to $6R_e$ in Serra et al. (2016), which reports $\gamma_{\text{tot}} = 2.18 \pm 0.03$ with an rms scatter 0.11.

## 5 DISCUSSION

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 $f_{\text{DM}}(R_e) = 0.07_{-0.07}^{+0.08}$. In this work, we measure $f_{\text{DM}}(R_e) = 0.05 \pm 0.06$ for NGC 3998 taking into account the mass-to-light gradient, consistent with the previous result. Since the stellar mass dominates the enclosed mass within $1R_e$ in NGC 3998, the measurements of $f_{\text{DM}}(R_e)$ are expected to be close to 0, regardless of the assumptions made.

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 H1 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 H1 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 H1 observations. Nevertheless, since this misaligned model reduces to the full aligned disc model when the angle between the stellar and H1 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 Fig. 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) $R_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 $R_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 $R_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 Fig. 5, and found a similar discrepancy between the observed and simulated galaxies, suggesting that the bias in $R_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 in 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 a 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 per cent, while this overestimation reduces to 18 per cent if a generalized-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 emphasize 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 generalized-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.

## 6 SUMMARY

In this paper, we map the mass distributions of 11 early-type galaxies with two-aperture stellar kinematics that extended to $2-4R_e$ using the Schwarzschild model, and we map the mass distributions of a subsample of four galaxies with H1 discs using dynamical models that combine two-aperture stellar and H1 gas kinematics extended to $10R_e$. Our main findings are summarized below:
(i) 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.

(ii) We robustly obtain the enclosed dark matter fraction within ~5R<sub>e</sub> for the 11 galaxies. The mean dark matter fraction within 1R<sub>e</sub> is 0.2 ± 0.1 (rms) and increases to 0.4 ± 0.2 (rms) within 2R<sub>e</sub> and 0.6 ± 0.2 (rms) within 5R<sub>e</sub> 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 5R<sub>e</sub>. The dark matter fractions of the sample galaxies are generally lower than the predictions of the TNG-100 and EAGLE simulations within 1R<sub>e</sub>, and there are 40 and 70 per cent sample galaxies with dark matter fractions that are 1σ consistent with either the TNG-100 or the EAGLE simulations at 2R<sub>e</sub> and 5R<sub>e</sub>, respectively, while the remaining 60 and 30 per cent galaxies lie below the 1σ region. Combined with dark matter fraction measurements out to 5R<sub>e</sub> in the literature, there are about 10 per cent of 47 galaxies lying below the 3σ region of the TNG-100 or the EAGLE predictions.

(iii) 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 M<sub>200</sub> and concentrated C are only weakly constrained with large uncertainties, generally following the M<sub>200</sub>–C correlations from the TNG-100 and the EAGLE simulations, and the relation presented in Dutton & Macciò (2014).

(iv) The total mass density profiles are proportional to r⁻², consistent with the previous dynamical results with data coverage to 4R<sub>e</sub> (Cappellari et al. 2015; Serra et al. 2016; Bellstedt et al. 2018).

## ACKNOWLEDGEMENTS

The authors thank the referee and the Scientific Editor for their helpful comments that helped to improve this manuscript.

The authors thank the referee for their helpful comments that helped to improve this manuscript. The research presented here is supported by the CAS Project for Young Scientists in Basic Research, Grant No. YSBR-062. MY acknowledges the China Postdoctoral Science Foundation, Grant No. 2021M703337. NFB acknowledges the Science and Technologies Facilities Council (STFC) grant ST/V000861/1.

The WRST is operated by the ASTRON (Netherlands Institute for Radio Astronomy) with support from the Netherlands Foundation for Scientific Research (NWO).

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No.

AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

## DATA AVAILABILITY

The data underlying in this article are available as follows: the MAT-LAS data at http://matlas.astro.unistra.fr/WP/, the PanSTARRS data at https://outerspace.stsci.edu/display/PANSTARRS/, the Atlas³D data at http://www-astro.physics.ox.ac.uk/atlas3d/, the Mitchell data at https://doi.org/10.17630/43b4e7fb-2f2c-4867-a4d5-21599299c9d9, the TNG-100 simulation at https://www.tng-project.org/ and the EAGLE simulation at https://eagle.strw.leidenuniv.nl/. The H I data underlying this article will be shared upon reasonable request to the corresponding authors.

## REFERENCES

Adams J. J. et al., 2014, ApJ, 789, 63

Alabi A. B. et al., 2016, MNRAS, 460, 3838

Alabi A. B. et al., 2017, MNRAS, 468, 3949

Bacon R. et al., 2001, MNRAS, 326, 23

Barnabè M. et al., 2012, MNRAS, 423, 1073

Battaglia G., Helmi A., Tolstoy E., Irwin M., Hill V., Jablonka P., 2008, ApJ, 681, L13

Bellstedt S. et al., 2018, MNRAS, 476, 4543

Binney J., Mamon G. A., 1982, MNRAS, 200, 361

Blumenthal G. R., Faber S. M., Primack J. R., Rees M. J., 1984, Nature, 311, 517

Boardman N. F. et al., 2016, MNRAS, 460, 3029

Boardman N. F. et al., 2017, MNRAS, 471, 4005

Bolton A. S., Treu T., Koopmans L. V. E., Gavazzi R., Moustakas L. A., Burles S., Schlegel D. J., Wayth R., 2008, ApJ, 684, 248

Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343

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

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

Cappellari M. et al., 2012, Nature, 484, 485

Cappellari M. et al., 2013, MNRAS, 432, 1709

Cappellari M. et al., 2015, ApJ, 804, L21

Cappellari M., 2002, MNRAS, 333, 400

Cappellari M., 2008, MNRAS, 390, 71

Cappellari M., 2013, ApJ, 778, L2

Cappellari M., 2017, MNRAS, 466, 798

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

Cappellari M., Emsellem E., 2004, PASP, 116, 138

Chambers K. C. et al., 2016, preprint (arXiv:1612.05560)

Cole D. R., Dehnen W., Wilkinson M. I., 2011, MNRAS, 416, 1118

Comerón S. et al., 2023, A&A, 675, A143

Cooke L. H. et al., 2022, MNRAS, 512, 1012

Côté P. et al., 2001, ApJ, 559, 828

Côté P., McLaughlin D. E., Cohen J. G., Blakeslee J. P., 2003, ApJ, 591, 850

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

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

Das P., Gerhard O., Mendez R. H., Teodorescu A. M., de Lorenzi F., 2011, MNRAS, 415, 1244

Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371

de Blok W. J. G., McGaugh S. S., Rubin V. C., 2001, AJ, 122, 2396

de Blok W. J. G., Walter F., Brinks E., Trachternach C., Oh S. H., Kennicutt R. C. J., 2008, AJ, 136, 2648

de Lorenzi F. et al., 2009, MNRAS, 395, 76

de Lorenzi F., Gerhard O., Saglia R. P., Sambhus N., Debattista V. P., Pannella M., Méndez R. H., 2008, MNRAS, 385, 1729

Deason A. J., Belokurov V., Evans N. W., McCarthy I. G., 2012, ApJ, 748, 2

Debuhr J., Quataert E., Ma C.-P., Hopkins P., 2010, MNRAS, 406, L55

Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604

Ding Y. et al., 2023, A&A, 672, A84
Duc P.-A. et al., 2015, MNRAS, 446, 120

Dumont A. et al., 2023, preprint (arXiv:2306.11786)

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

Emsellem E. et al., 2004, MNRAS, 352, 721

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

Forestell A. D., Gebhardt K., 2010, ApJ, 716, 370

Fukazawa Y., Botoya-Nonesa J. G., Pu J., Ohto A., Kawano N., 2006, ApJ, 636, 698

Genel S. et al., 2018, MNRAS, 474, 3976

Genzel R. et al., 2020, ApJ, 902, 98

Harris G. L. H., Babyk I. V., Harris W. E., McNamara B. R., 2019, ApJ, 887, 259

Harris W. E., Remus R.-S., Harris G. L. H., Babyk I. V., 2020, ApJ, 905, 28

Hayashi K., Chiba M., Ishiyama T., 2020, ApJ, 904, 45

Hill G. J. et al., 2008, in McLean I. S., Casali M. M. eds, SPIE Conf. Ser. Vol 7014, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, Bellingham, p.701470

Hopkins P. F. et al., 2018, MNRAS, 480, 800

Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., Robertson B., Springel V., 2006, ApJS, 163, 1

Jethwa P., Thater S., Maindl T., Van de Ven G., 2020, DYNAMITE: Dynamics, Age and Metallicity Indicators Tracing Evolution, Astrophysics Source Code Library, record ascl:2011.007

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

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

Józsa G. I. G., Oosterloo T. A., Morganti R., Klein U., Erben T., 2009, A&A, 494, 489

Kaplinghat M., Tulin S., Yu H.-B., 2016, Phys. Rev. Lett., 116, 041302

Krajnović D., Cappellari M., de Zeeuw P. T., Copin Y., 2006, MNRAS, 366, 787

Kroupa P., 2001, MNRAS, 322, 231

Kuzio de Naray R., McGaugh S. S., de Blok W. J. G., 2008, ApJ, 676, 920

Leier D., Ferreras I., Negri A., Saha P., 2022, MNRAS, 510, 24

Lelli F., McGaugh S. S., Schombert J. M., 2016, AJ, 152, 157

Li C. et al., 2020, MNRAS, 492, 2775

Li H. et al., 2017, ApJ, 838, 77

Lovell M. R. et al., 2018, MNRAS, 481, 1950

Marinacci F. et al., 2018, MNRAS, 480, 5113

Martinsson T. P. K., Verheijen M. A. W., Westfall K. B., Bershady M. A., Andersen D. R., Waters R. A., 2013, A&A, 557, A131

Morganti L., Gerhard O., Coccato L., Martinez-Valpuesta I., Arnaboldi M., 2013, MNRAS, 431, 3570

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

Napolitano N. R. et al., 2009, MNRAS, 393, 329

Napolitano N. R. et al., 2011, MNRAS, 411, 2035

Napolitano N. R., Pota V., Romanowsky A. J., Forbes D. A., Brodie J. P., Foster C., 2014, MNRAS, 439, 659

Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72

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

Nelson D. et al., 2018, MNRAS, 475, 624

Nestor Shachar A. et al., 2023, ApJ, 944, 78

Newman A. B., Treu T., Ellis R. S., Sand D. J., 2013, ApJ, 765, 25

Oguri M., Rusu C. E., Falco E. E., 2014, MNRAS, 439, 2494

Oh S.-H. et al., 2015, AJ, 149, 180

Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt Robert C. J., 2011, AJ, 141, 193

Oh S.-H., de Blok W. J. G., Walter F., Brinks E., Kennicutt Robert C. J., 2008, AJ, 136, 2761

Oldham L. J., Auger M. W., 2018, MNRAS, 476, 133

Ostriker E. C., Shetty R., 2011, ApJ, 731, 41

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

Poci A., Cappellari M., McDermid R. M., 2017, MNRAS, 467, 1397

Read J. I., Gilmore G., 2005, MNRAS, 356, 107

Read J. I., Walker M. G., Steger P., 2018, MNRAS, 481, 860

Read J. I., Walker M. G., Steger P., 2019, MNRAS, 484, 1401

Relatores N. C. et al., 2019, ApJ, 887, 94

Robertson A. et al., 2018, MNRAS, 476, L20

Romanowsky A. J., Douglas N. G., Arnaboldi M., Kuijken K., Merrifield M. R., Napolitano N. R., Capaccioli M., Freeman K. C., 2003, Science, 301, 1696

Salo H. et al., 2015, ApJS, 219, 4

Salpeter E. E., 1955, ApJ, 121, 161

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

Sand D. J., Treu T., Smith G. P., Ellis R. S., 2004, ApJ, 604, 88

Santini P. et al., 2012, A&A, 538, A33

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

Sault R. J., Teuben P. J., Wright M. C. H., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E. eds, ASP Conf. Ser. Vol. 77, Astronomical Data Analysis Software and Systems IV. Astron. Soc. Pac., San Francisco, p. 433

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

Schwarzschild M., 1979, ApJ, 232, 236

Serra P. et al., 2012, MNRAS, 422, 1835

Serra P., Oosterloo T., Cappellari M., den Heijer M., Józsa G. I. G., 2016, MNRAS, 460, 1382

Sharma G., Salucci P., van de Ven G., 2021, A&A, 653, A20

Sharma G., Salucci P., van de Ven G., 2022, A&A, 659, A40

Simon J. D., Bolatto A. D., Leroy A., Blitz L., 2003, ApJ, 596, 957

Sonnenfeld A., Treu T., Marshall P. J., Suyu S. H., Gavazzi R., Auger M. W., Nipoti C., 2015, ApJ, 800, 94

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

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

Strickland D. K., Heckman T. M., 2009, ApJ, 697, 2030

Syer D., Tremaine S., 1996, MNRAS, 282, 223

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

Tiley A. L. et al., 2019, MNRAS, 485, 934

Tollet E. et al., 2016, MNRAS, 456, 3542

Tonry J. L., 1983, ApJ, 266, 58

Tortora C., Posti L., Koopmans L. V. E., Napolitano N. R., 2019, MNRAS, 489, 5483

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

Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20

Wang Y. et al., 2020, MNRAS, 491, 5188

Wasserman A. et al., 2018, ApJ, 863, 130

Weijmans A.-M. et al., 2009, MNRAS, 398, 561

Yang M., Zhu L., Weijmans A.-M., van de Ven G., Boardman N., Morganti R., Oosterloo T., 2020, MNRAS, 491, 4221

Yıldız M. K., Serra P., Oosterloo T. A., Peletier R. F., Morganti R., Duc P.-A., Cuillandre J.-C., Karabal E., 2015, MNRAS, 451, 103

Yuan S., Zhu L., Liu C., Qu H., Fan Z., 2022, Res. Astron. Astrophys., 22, 085023

Zanisi L. et al., 2021, MNRAS, 501, 4359

Zhao H., 1996, MNRAS, 278, 488

Zhou S. et al., 2019, MNRAS, 485, 5256

Zhu K., Lu S., Cappellari M., Li R., Mao S., Gao L., 2023, MNRAS, 522, 6326

Zhu K., Lu S., Cappellari M., Li R., Mao S., Gao L., Ge J., 2024, MNRAS, 527, 706

Zhu L. et al., 2014, ApJ, 792, 59

Zhu L. et al., 2016, MNRAS, 462, 4001

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

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

### suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.
## APPENDIX A: HI VELOCITY FIELDS

We exhibit the H1 velocity fields of the subsample obtained with WSRT in Fig. A1.

<figure><img src="image_12.png" /><figcaption>**Figure A1.** H1 velocity fields of the subsample obtained with the WSRT. The median velocity of each galaxy is already subtracted.</figcaption></figure>

An tilted inner ring is located in the centre of NGC 2685 (Józsa et al. 2009), so we mask the central region within an ellipse with a radius of 120 arcsec.

There are some areas with discontinuous velocities in the velocity fields of NGC 4203, which are probably H1 blobs located outside the discs. Therefore, we mask these regions in our fitting of the misaligned gas models.

There are extended arms in the velocity fields of NGC 4203 and NGC 5582, which are located outside $10R_e$ of each galaxy. As these arms might be biased from the nearly dynamical equilibrium, we mask the regions outside $10R_e$ for these galaxies.

## APPENDIX B: CONFIDENCE LEVEL CALCULATION

We first define an 'ideal' residual of the best-fitting model by

$$ \chi_{\text{ideal}}^2 = \sum_{i}^{N} \left( \frac{\mu_{\text{mod},i} - \mu_{\text{data},i}}{\epsilon_{\text{data},i}} \right)^2, \qquad (\text{B1}) $$

where $f_{\text{mod},i}$ is the best-fitting modelled data, $\mu_{\text{data},i}$ and $\epsilon_{\text{data},i}$ is the 'real' data and its uncertainty, and $N$ is the number of independent data points. We assume the observation uncertainty $\epsilon_{\text{obs},i} = \epsilon_{\text{data},i}$, then possible observational data follow a Gaussian distribution that $f_{\text{obs},i} \sim N(\mu_{\text{data},i}, \epsilon_{\text{obs},i}^2)$.

The residual distribution for the possible observation data is obtained by

$$ \begin{aligned} \chi_{\text{f}}^2 &= \sum_{i}^{N} \left( \frac{\mu_{\text{mod},i} - f_{\text{obs},i}}{\epsilon_{\text{obs},i}} \right)^2 \\ &= \chi_{\text{ideal}}^2 + \sum_{i}^{N} x_i^2 - 2 \sum_{i}^{N} \left( \frac{f_{\text{mod},i} - \mu_{\text{data},i}}{\epsilon_{\text{obs},i}} \right) x_i, \end{aligned} \qquad (\text{B2}) $$

where each $x_i \sim N(0, 1)$. Its expectation $\mathbb{E}(\chi_{\text{f}}^2) = \chi_{\text{ideal}}^2 + N$, and its standard deviation $\sigma(\chi_{\text{f}}^2) = \sqrt{4\chi_{\text{ideal}}^2 + 2N} = \sqrt{4\mathbb{E}(\chi_{\text{f}}^2) - 2N}$.

For the given observational data, we expect its residual $\chi_{\text{obs}}^2 \sim \mathbb{E}(\chi_{\text{f}}^2)$. Therefore, we adopt $\Delta\chi_{\text{obs}}^2 = \sigma(\chi_{\text{f}}^2) \sim \sqrt{4\chi_{\text{obs}}^2 - 2N}$ as the $1\sigma$ confidence interval due to observation uncertainties.

This paper has been typeset from a TeX/Latex file prepared by the author.