Is the Core-cusp Problem a Matter of Perspective? Jeans Anisotropic Modeling against Numerical Simulations

Wenting Wang; Ling Zhu; Zhaozhou Li; Yang Chen; Jiaxin Han; Feihong He; Xiaohu Yang; Yipeng Jing; Carlos Frenk; Jialu Nie; Hao Tian; Chao Liu; Yanan Cao; Xiaoqing Qiu; John Helly; Robert J. J. Grand; Facundo A. Gomez

doi:10.3847/1538-4357/ac9b19

1. Introduction

The growth of large and intermediate scales of cosmic structures in our universe, such as cosmic filaments and the distribution of clumpy dark matter halos, can be modeled by the linear perturbation theory under the standard Lambda cold dark matter (ΛCDM) cosmological model, which has turned out to be remarkably successful (e.g., Yang et al. 2004; Cole et al. 2005; Henriques et al. 2012; Han et al. 2015; Wang et al. 2016; Springel et al. 2018). On small scales, galaxies form through the gas cooling and condensation within dark matter halos (e.g., White & Rees 1978). Smaller halos and galaxies can merge with larger halos, becoming the so-called substructures and satellite galaxies. On such small scales within dark matter halos, a series of challenges to the standard theory has been raised (e.g., Klypin et al. 1999; Boylan-Kolchin et al. 2011; Merritt et al. 2020), among which one famous and hotly debated issue is the core-cusp problem. Dark-matter-only simulations predict inner density slopes close to −1 (cusp), whereas the modeling of gas rotation curves or stellar kinematics in the central regions of low surface brightness galaxies, gas-rich dwarfs, and dwarf spheroids favor inner slopes close to 0 (core), which brings in tension with the theory (e.g., Flores & Primack 1994; Moore 1994; de Blok et al. 2001; Gentile et al. 2004; de Blok 2010; also see Bullock & Boylan-Kolchin 2017 for a review).

The core-cusp problem has invoked significant interest. It motivates alternative dark matter models such as the self-interacting dark matter (SIDM) model (Peter et al. 2013; Rocha et al. 2013; Oman et al. 2015; Foot & Vagnozzi 2015). Other promising solutions to the problem within the ΛCDM framework include, for example, the stellar feedback, that dark matter is heated by stellar winds or supernovae explosions, which drive repeated gravitational potential fluctuations and the dark matter particle orbits slowly expand (e.g., Read & Gilmore 2005; Read et al. 2019; Freundlich et al. 2020a; Li et al. 2022b; Boldrini 2021). Modern numerical simulations have shown that this is possible (e.g., Mashchenko et al. 2008; Pontzen & Governato 2012, 2014). However, stellar feedback is only shown to be efficient for dwarfs with stellar-to-halo mass ratios between 10⁻³ and 10⁻², whereas ultra-faint dwarfs with this ratio smaller than 10^−3.5 have formed very few stars, and thus stellar feedback seems very unlikely to be the main mechanism forming cores (e.g., Tollet et al. 2016; Freundlich et al. 2020b; Hayashi et al. 2020; Lazar et al. 2020). In addition to stellar feedback, tidal effects have been shown to be capable of forming cores in numerical simulations, even through galaxy interactions before infalling (e.g., Genina et al. 2022).

Note that core-like inner density profiles are only reported in some dwarf galaxies, and the constraints could still be limited by statistical error or by model extrapolations (e.g., Zhu et al. 2016a; Lazar & Bullock 2020; Shi et al. 2021). Observation also demonstrates a strong diversity in the rotation curves, even at fixed maximum circular velocity (e.g., Read et al. 2019; Hayashi et al. 2020; Kaplinghat et al. 2020; Santos-Santos et al. 2020), and the diversity is shown to be related to the stellar-to-halo mass ratio (Hayashi et al. 2020), to the star formation history (SFH; Read et al. 2019), and to the surface brightness (e.g., Santos-Santos et al. 2020). However, it seems that CDM hydrodynamical simulations cannot fully account for the diversity (Kaplinghat et al. 2020; Santos-Santos et al. 2020).

In addition to the purpose of serving as theoretical predictions under the ΛCDM framework and being compared with observational constraints, modern numerical simulations are perhaps the most useful source of data to help validate various assumptions behind dynamical modeling methods of constraining the inner dark matter profiles. It is impossible to directly infer the level of systematic uncertainties behind dynamical models based on purely observational data, but numerical simulations provide us the ground truth to compare with. For example, Walker & Peñarrubia (2011) used two stellar populations with different half-mass radii to infer the dark matter densities of Sculptor and Fornax. Both galaxies are claimed to have cored dark matter halos. However, based on the APOSTLE suite of hydrodynamic simulations of Local Group analogs, Genina et al. (2018) concluded that violations of the spherical symmetry can mistakenly result in best-fitting core profiles, when the truth is cuspy. In a later study by Genina et al. (2020), a spherical Jeans modeling method was further tested using both the CDM and SIDM models. For the CDM model, they reported 50% and 20% of scatter for the enclosed mass in inner regions or in the half-mass radius of tracers, respectively. For SIDM dwarfs, their recovered mass profiles are biased toward cuspy dark matter distributions. Campbell et al. (2017) also estimated the systematics of half-mass estimators developed in previous studies (Wolf et al. 2010; Walker & Peñarrubia 2011), and reported intrinsic scatters of 23%–25%.

In this study, we first construct mock observations for the member stars of 28 dwarf systems, selected from the cosmological auriga suite of simulations. We then validate the discrete axisymmetric Jeans Anisotropic Multi-Gaussian Expansion (jam) method based on the mock observational data. jam is widely used in estimating the underlying mass profiles for different types of galaxies (e.g., Cappellari et al. 2013; Li et al. 2017) and for dwarf galaxies in the Milky Way (e.g., Watkins et al. 2013). In a previous study, jam is well tested using mock IFU data based on an Illustris simulation (Li et al. 2016), and for a few mock dwarf galaxies with discrete data and in steady status (Zhu et al. 2016a). Here, our mock tracer stars and dwarf galaxy systems can more closely represent the dynamical status of real observed dwarf systems in the Milky Way, hence enabling us to investigate the performance of jam with realistic nonsteady tracers.

We not only investigate the performance of jam but also the diversity in the best fits. We find that the amount of bias is correlated with the current specific star formation rates (SFRs), which is due to the difference in the dynamical status of quiescent and star-forming dwarfs. While radial scales most relevant to the core-cusp problem are already below the resolution of the simulations, we extrapolate to the very inner radii to draw more general inferences. We check the performance of jam by considering both the error-free data and the data after incorporating realistic observational errors as well as contamination by unbound stars. One purpose of this study is to investigate whether current Gaia and future China Space Station Telescope (CSST) proper motions can be useful for such dynamical modelings.

In the following, we first introduce the auriga suite of simulations, our sample of dwarfs, mock stars, and the method of incorporating observational errors and modeling the contamination by foreground/background or unbound stars in Section 2. The jam modeling approach is introduced in Section 3. The results are presented in Section 4, which include the overall model performance and dependence of the systematic biases on the SFR and dynamical status of the systems. We also discuss the results after incorporating observational errors + unbound star contaminations, and we discuss the results with or without proper motions and with a smaller tracer sample size. We make connections to the core-cusp problem and discuss the detectability of contraction motions in member stars in Section 5. We conclude in Section 6.

2. Data

2.1. Sample of Dwarf Systems, Mock Stars, and Mock Images

2.1.1. auriga Suite of Simulations

The sample of dwarf galaxies is constructed from the auriga suite of simulations (Grand et al. 2017). Details of the auriga simulations can be found in Grand et al. (2017) and Grand et al. (2018). Here, we make a brief introduction.

The auriga simulations are a set of cosmological zoom-in simulations. The evolution of Milky Way mass systems are simulated and traced from redshift z = 127 to 0. They are identified as isolated halos from the parent dark-matter-only simulations of the EAGLE project (Schaye et al. 2015). The cosmological parameters adopted are from the third year of Planck data (Planck Collaboration et al. 2014) with Ω_m = 0.307, Ω_Λ = 0.693, Ω_b = 0.048, and H₀ = 67.77 km s⁻¹ Mpc⁻¹.

The simulations were conducted using the magnetohydrodynamic code arepo (Springel 2010) with full baryonic physics, which incorporates a comprehensive galaxy formation model and has higher resolutions than the parent simulation. The physical mechanisms of the galaxy formation model include atomic and metal line cooling (Vogelsberger et al. 2013), a uniform UV background (Faucher-Giguère et al. 2009), a subgrid model of the interstellar medium and star formation processes (Springel & Hernquist 2003), metal enrichment from supernovae and AGB stars (Vogelsberger et al. 2013), feedback from core-collapse supernovae (Okamoto et al. 2010), and the growth and feedback from supermassive black holes (Springel et al. 2005). A uniform magnetic field with comoving strength of 10⁻¹⁴ G is set at redshift z = 127, which quickly becomes subdominant in collapsed halos (Pakmor & Springel 2013; Pakmor et al. 2017).

In this study, we use six Milky Way–like systems from the level 3 set of simulations. The six systems are named Au6, Au16, Au21, Au23, Au24, and Au27. The virial masses¹² of their host dark matter halos are in the range of 1–2 × 10¹² M_⊙. The typical dark matter particle mass is about 4 × 10⁴ M_⊙, while the average baryonic particle mass is about 5 × 10³ M_⊙.

2.1.2. Dwarf Systems

Each of the six systems has dwarf satellite galaxies, and in our analysis, we only use those dwarf systems that are less massive than 10⁹ M_⊙ in stellar mass and also have more than 6000 star particles. Here, the upper limit of 10⁹ M_⊙ is adopted to avoid including dwarf galaxies that are significantly more massive than classical dwarf spheroid galaxies analyzed in previous studies (see, e.g., Hayashi et al. 2020). Besides, as discussed by Zhu et al. (2016a), a discrete data set with 6000 stars is required to distinguish between core and cusp inner profiles. Thus, we use such a large sample of stars as tracers, in order to control the statistical errors to be small enough, while we focus on discussing systematic errors behind the model, but for part of our analysis, we will try a smaller tracer sample of 2000 stars. Note, however, each star particle corresponds to a single stellar population, whose particle mass ranges on average from ∼4600 to ∼6500M_⊙ for different simulations. A lower limit of 6000 particles roughly corresponds to lower limits in the total stellar mass of ∼10^7.44 to ∼10^7.59 M_⊙. We also exclude those dwarfs whose dynamical spin axes and geometrical minor axes have strong misalignments (see explanations below in this section). In the end, we have 28 dwarf systems in total. For each dwarf system, we randomly pick 6000 bound star particles as tracers. Note that as we have explicitly checked, different random selections of the tracer sample lead to differences smaller than the symbol sizes in the relevant figures in this study.

Out of the 28 dwarfs, we will explicitly show the best-fitting and true density profiles for six systems similar to the Sagittarius dwarf spheroid (dSph)¹³ and other four representative star-forming systems with prominent galactic winds or gas outflows.¹⁴ Table 1 summarizes the information on these systems. Interestingly, only one out of the six Sagittarius dSph–like (hereafter Sag-like) systems is as old as 10^11.2 Gyr, whereas the other three systems still have star formations at ∼6–9 Gyr ago, indicating such massive satellites can persist star formations after falling to the current host halo. Compared with the nearby Sag-like systems, other star-forming systems with strong outflows are at relatively large distances. This is consistent with the conclusion of Simpson et al. (2018), who reported strong mass and distance-dependent quenching signals in satellites of Local Group–like systems from 30 zoom-in simulations.

Table 1. Information on the Six Sag-like Systems and the Four Star-forming Systems with Prominent Gas Outflows from the Level 3 Set of auriga Simulations

Host Name	Dwarf	${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })$	D [kpc]	Age (Gyr)
Au16	9	7.933	24.53	11.20
Au21	10	8.446	42.29	6.75
Au23	4	8.297	45.74	7.94
Au23	7	8.047	34.52	8.96
Au24	24	7.503	60.15	9.36
Au27	25	7.598	27.09	9.57

Au21	7	7.287	449.04	11.22
Au24	9	7.842	235.10	10.11
Au24	13	7.524	261.72	10.25
Au27	3	8.054	385.27	9.66

Note. The dwarf ID is provided in the second column, which is simply the position index of the dwarf in the corresponding subhalo catalog.

Download table as: ASCII Typeset image

2.1.3. Mock Stars

To create mock samples of observed stars in each dwarf system, we start from the original coordinates in the simulation boxes and subtract the stellar-mass weighted mean coordinates and velocities of all bound particles belonging to the dwarfs, to eliminate the perspective rotation (Feast et al. 1961). We place the observer on the disk plane, which is defined as the plane perpendicular to the minor axis of all bound star particles with Galactocentric distances smaller than 20 kpc. The observer is 8 kpc away from the Galactic Center, with a random position angle.

The coordinates and velocities are then transformed to the observing frame. The $z^{\prime}$ -axis of the observing frame is chosen as the line-of-sight direction. The $x^{\prime}$ -axis (major axis) is the cross product between the spin axis of the dwarf galaxy and the $z^{\prime}$ -axis, which is projected on the sky. The $y^{\prime}$ -axis (minor axis) is the cross product between the $z^{\prime}$ and $x^{\prime}$ vectors, taking a minus sign. We introduce the minus sign here because, according to Watkins et al. (2013), the observing frame of the jam model is a left-handed system.

2.1.4. Mock Dwarf Images and Multi-Gaussian Expansion

In jam modeling, the potential and density distribution of the luminous stellar component will be directly inferred from the projected optical images of the dwarfs, with the stellar mass-to-light ratios (M/L) being free parameters. The image will be deprojected to be in three dimensions based on the distance and inclination angle of the dwarf. The inclination angle is defined as the angle between the average spin axis of the dwarf and the line-of-sight $z^{\prime}$ direction, with the value between 0° and 180°.

Thus, to apply jam we need to create mock images for our sample of dwarfs. We simply adopt the projected stellar-mass density distribution to create the images, i.e., the read in each pixel is in units of stellar mass per square parsec based on all bound star particles associated with the dwarf galaxy, so in our case the true value of the M/L is unity.

Once the mock images are made, the luminous stellar-mass distribution, ${\rm{\Sigma }}(x^{\prime} ,y^{\prime} )$ , will be decomposed to a few different Gaussian components (multi-Gaussian expansion (MGE)), in order to enable the analytical deprojection for any arbitrary ${\rm{\Sigma }}(x^{\prime} ,y^{\prime} )$ and to bring analytical solutions for any arbitrary matter distribution (see Section 3 for more details). In fact, the process of creating MGEs requires the $x^{\prime}$ -axis to be defined as the longer axis of the galaxy image. For ideally axisymmetric systems, the $x^{\prime}$ -axis as we defined above is expected to be the projected longer axis as well. However, real triaxial systems might deviate from this because their minor and spin axes might misalign. As we have checked, most of the dwarfs in our analysis are not strongly triaxial, and their image $x^{\prime}$ -axis is indeed aligned with the long axis, with a misaligned angle of at most ∼15°. In a few cases, the $x^{\prime}$ -axis might be almost perpendicular to the actual long axis, but the galaxy images are very close to be spherical, with a very tiny difference between the major and minor axes. In a few very extreme cases, the dwarfs are strongly elongated, whose spin axes prominently deviate from their minor axes, and the $x^{\prime}$ -axes are very different from the image longer axis. These systems usually significantly deviate from axisymmetry and are out of equilibrium (not in steady states). The performance of jam is expected to be very poor once applied to such systems (e.g., Li et al. 2016), regardless of how the $x^{\prime}$ -axis is defined. Thus, we have excluded them from our analysis. Observationally, such extremely triaxial systems can also be excluded by comparing their image minor axis and the spin axis inferred from line-of-sight velocity maps.

2.2. Incorporating Observational Errors

In our analysis, we first use true positions and velocities of 6000 bound star particles without including realistic errors. In this way, we can better understand the intrinsic systematic uncertainties behind the model. We will then repeat our analysis for only 2000 bound star particles with only line-of-sight velocities plus typical errors. This is for more practical meanings because the number of observed member stars in their host dwarfs is at most ∼2000 for current observations, most of which do not have accurate proper motion measurements. Moreover, for the six nearby Sag-like systems, we will incorporate realistic observational errors to both line-of-sight velocities and proper motions, and model the contamination of foreground/background or unbound stars, with member stars selected based on the differences in their kinematics relative to the dwarf centers. This will be compared with the error-free and foreground-/background-free model performance. Lastly, we find most Sag-like systems are undergoing contractions due to the effect of tidal force. We will check whether the amount of contraction is observable after including the Gaia DR3 proper motion errors. We now introduce how the observational errors are incorporated and how the foreground/background is modeled.

2.2.1. Assigning Apparent Magnitudes of Individual Stars

Star particles in current state-of-the-art simulations represent simple stellar populations, and thus we cannot directly use their stellar properties, such as the magnitudes and colors, to represent single stars, though we can still use the position and velocity information of star particles in the simulations. To assign each star particle the appropriate magnitude information and hence observational errors, we use the trilegal code (Girardi et al. 2005; Vanhollebeke et al. 2009; Girardi et al. 2012; Girardi 2016, Y. Chen et al. 2022, in preparation, Dal Tio et al. 2022) to generate populations of different types of stars based on the input SFH, the age–metallicity relation (AMR), and the total stellar mass. The SFH of the Sagittarius dSph is taken from Weisz et al. (2014), and we calculate the AMR based on the closed-box model, which agrees well with real data (e.g., Layden & Sarajedini 2000). For the other dwarfs, the magnitude information is only required when we generate realistic errors for the line-of-sight velocities of mock stars (the case of using 2000 star particles as tracers). For them, we simply use the same SFH and AMR as that of the Sagittarius dSph. This is approximately a reasonable choice for dwarf galaxies because dwarfs are relatively old, and the differences in their SFH and AMR do not lead to significant variations in the magnitude distribution of member stars.

The filter response curve is for the CSST (Zhan 2011; Cao et al. 2018; Gong et al. 2019) g filter. The absolute magnitudes generated by trilegal are converted to apparent magnitudes according to the distance modulus of each dwarf galaxy with respect to the mock observer in the simulation (see Section 2.1 above). In our analysis, we will incorporate typical proper motion errors from either Gaia DR3 or from CSST.¹⁵ When adopting the Gaia errors, we convert CSST g magnitudes to Gaia G magnitudes based on an empirical and color-dependent relation linking the g filter of the Hyper Suprime-Cam survey (Miyazaki et al. 2018; Komiyama et al. 2018; Furusawa et al. 2018) to Gaia G (Wang et al. 2019), assuming that the CSST g filter is not very different from the HSC g filter. This is a reasonable approximation, according to the CSST and HSC filter response curves (Kawanomoto et al. 2018; Gong et al. 2019).

The total number of stars generated by trilegal is very large, while the number of star particles used as tracers and bound to each dwarf galaxy in the simulation is more limited (in our case we use a subset of 6000 or 2000 as mentioned above). We take those trilegal stars whose apparent magnitudes are in a certain range of apparent magnitudes at the corresponding distance of each dwarf, randomly select 6000 or 2000 of them, and assign their magnitudes to our tracer star particles. The magnitude range is chosen to be 10 < g < 19 (Gaia DR3 proper motion errors) or 18 < g < 19 (CSST proper motion errors) for the Sag-like systems. These nearby systems can all have more than 6000 stars with g < 19. The magnitude range is 18 < g < 19 when using CSST proper motion errors because stars brighter than g = 18 will be saturated in CSST observations. In particular, when we repeat our analysis with the smaller sample of 2000 star particles as tracers (line-of-sight velocities only and with typical errors), we only use dwarf galaxies that can have more than 2000 mock stars observed above the given apparent magnitude thresholds.¹⁶ For this case, we choose the magnitude range of 10 < g < 19 for nearby Sag-like systems, and the range is chosen as 10 < g < 21 for more distant systems.

2.2.2. Parallax, Line-of-sight Velocity, and Proper Motion Errors

According to the apparent magnitudes, we assign each star particle errors in parallax, line-of-sight velocity, and proper motion. The parallax errors are the median Gaia DR3 errors at the corresponding magnitude (Gaia Collaboration et al. 2016, 2021). The errors of the line-of-sight velocities are linear interpolations between 1 and 10 km s⁻¹, assuming stars with g = 10 have 1 km s⁻¹ of errors in the line-of-sight velocity, while stars with g = 21 have a 10 km s⁻¹ error. The line-of-sight velocity errors are typical for current or future spectroscopic surveys such as the Dark Energy Spectroscopic Instrument (Yano et al. 2016; Allende Prieto et al. 2020; Rehemtulla et al. 2022). The errors in proper motions are either the median Gaia DR3 errors or the typical CSST errors at the corresponding magnitudes.

For CSST, details about how its proper motion errors are modeled and predicted will be presented in a separate paper (J. Nie et al. 2022, in preparation). Nie et al. generated mock stars in the Galactic bulge region, and realistic star images are created at different epochs after including the point-spread function (PSF), noise due to correction of flat and bias, readout noise, sky background, parallax, and stellar motions. The process of source detection and centroid position determination at different epochs are then modeled. After calibrating the reference frame using mock reference stars, the final astrometric solutions and associated errors are obtained. The reference stars are assumed to have proper motions available from Gaia. Details of the distribution of reference stars, relevant CSST instrument parameters, and how the mock star samples, stellar motions, and PSF are modeled will be found in Nie et al. Moreover, their PSF model is simplified, which might be different from the real PSF during real operations. Based on about six astrometric measurements evenly distributed in 10 yr of baseline, the typical proper motion error is ∼0.2 mas yr⁻¹ at 18 < g < 19. Note that because the number density of member stars in dwarf galaxies is less dense than the Galactic bulge region, their error is likely an upper limit. For Gaia, the median errors are obtained by querying the Gaia database, with a typical DR3 proper motion error of ∼0.1 mas yr⁻¹ at 18 < g < 19.

After realistic observational errors are obtained, the proper motion, distance, and line-of-sight velocity of each star particle are displaced from their true values according to the assigned errors, by adding a displacement drawn from a Gaussian distribution with the standard deviation equaling the corresponding error.

2.2.3. Modeling of Foreground/background and Unbound Star Contamination

After observational errors are incorporated, we select as tracers those star particles whose differences in line-of-sight velocities, proper motions, and parallaxes with respect to the dwarf centers¹⁷ are all smaller than five times the corresponding errors (e.g., Vasiliev & Belokurov 2020). In this way not only those true bound star particles are selected as tracers, but also some foreground and background particles as well as unbound particles around the dwarfs are included. To assign magnitudes for foreground/background and unbound particles, we adopt the sample of stars created by trilegal for the Milky Way, and only use those stars in a rectangular region approximately centered on the observed Galactic latitude and longitude of the Sagittarius dSph (0° < l < 11° and −20° < b <−8°). We randomly select a subsample of trilegal stars and assign their magnitudes to foreground/background and unbound star particles in the simulation. The observational errors for foreground/background and unbound particles are generated in the same way as for bound member star particles based on their apparent magnitudes.

Notably, for Sag-like systems, almost all bound star particles will be selected as tracers in this way, and we ensure these bound tracers to be the same as the subset used in the error-free case, for fair comparisons. Other selected unbound star particles are randomly picked in proportion to the number of the selected fraction of bound star particles.

Figure 1 shows the differences between the proper motions after assigning errors and the true proper motions (Au16-9), as a function of the apparent magnitudes and for both the Gaia DR3 and CSST proper motion errors. The scatter is a factor of 2 smaller for Gaia DR3 errors than CSST errors at 18 < g < 19, and is smaller than 0.1 mas yr⁻¹ at g < 18 for Gaia DR3. The dense clump of stars in the right panels and at g ∼ 18.4 correspond to the red giant branch bump, whereas there are no such prominent features in the sample of foreground/background stars in the left panels.

$x^{\prime} $ — **Figure 1.** Differences between the proper motions after incorporating realistic observation errors and the true proper motions, reported as a function of the CSST g-band apparent magnitudes. We only show the $x^{\prime}$ component, as the $y^{\prime}$ component is very similar. Black (left) and red (right) points denote foreground/background stars and member stars of Au16-9 in the simulation, respectively. The proper motion errors are for either Gaia DR3 (top panels) or for CSST (bottom panels). The dense clumps at g ∼ 18.4 in the right panels correspond to the red giant branch bump. Gaia DR3 proper motion errors are about a factor of 2 smaller than the upper limits in CSST errors at 18 < g < 19.
Download figure:
Standard image High-resolution image

3. Methodology

jam is a publicly available source of code.¹⁸ It is a powerful tool to constrain both the underlying matter distribution and the internal dynamics of tracers (e.g., Zhu et al. 2016b), based on either line-of-sight velocities or proper motions of tracers. The version of jam we use for this paper is slightly different from the public version of the jam model for discrete data (Watkins et al. 2013), with an improved python interface and plotting tools. Details about jam can be found in Cappellari (2008) and Watkins et al. (2013), and here we only briefly introduce the method.

The method is based on solving the axisymmetric Jeans equation in an intrinsic frame defined on the dwarf galaxy with cylindrical coordinates, to solve for the first and second velocity moments:

$\begin{eqnarray}&&\displaystyle \frac{\nu (\overline{{v}_{R}^{2}}-\overline{{v}_{\phi }^{2}})}{R}+\displaystyle \frac{\partial \nu \overline{{v}_{R}^{2}}}{\partial R}+\displaystyle \frac{\partial \nu \overline{{v}_{R}{v}_{z}}}{\partial z}=-\nu \displaystyle \frac{\partial {\rm{\Phi }}}{\partial R},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\displaystyle \frac{\nu \overline{{v}_{R}{v}_{z}}}{R}+\displaystyle \frac{\partial \nu \overline{{v}_{R}{v}_{z}}}{\partial R}+\displaystyle \frac{\partial \nu \overline{{v}_{z}^{2}}}{\partial z}=-\nu \displaystyle \frac{\partial {\rm{\Phi }}}{\partial z},\end{eqnarray} \tag{ 2 }$

where ν is the tracer density distribution. Φ is the total potential. Upon solving the equation to obtain unique solutions, the cross-velocity terms are assumed to be zero, i.e., $\overline{{v}_{R}{v}_{z}}=0$ . In addition, the anisotropy parameter, b, is assumed to be constant and defined as $\overline{{v}_{R}^{2}}=b\overline{{v}_{z}^{2}}$ . A rotation parameter, κ, is introduced as $\overline{{v}_{\phi }}=\kappa {\left(\overline{{v}_{\phi }^{2}}-\overline{{v}_{R}^{2}}\right)}^{1/2}$ .

In our analysis, we define the z-axis of the intrinsic frame as the direction of the averaged spin of all bound star particles to the dwarf in the simulation, and the intrinsic frame is a right-handed system. The intrinsic frame is linked to the observing frame (see Section 2.1 above) through the inclination angle, i, of the dwarf galaxy:

$\begin{eqnarray}&&\left(\begin{array}{c}x^{\prime} \\ y^{\prime} \\ z^{\prime} \end{array}\right)=\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & -\cos i & \sin i\\ 0 & \sin i & \cos i\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right),\end{eqnarray} \tag{ 3 }$

and

$\begin{eqnarray}&&\left(\begin{array}{c}{v}_{x^{\prime} }\\ {v}_{y^{\prime} }\\ {v}_{z^{\prime} }\end{array}\right)=\left(\begin{array}{ccc}\cos i & -\sin i & 0\\ \sin i & \cos i & 0\\ 0 & 0 & 1\end{array}\right)\left(\begin{array}{c}{v}_{R}\\ {v}_{\phi }\\ {v}_{z}\end{array}\right),\end{eqnarray} \tag{ 4 }$

where $R=\sqrt{{x}^{2}+{y}^{2}}$ .

The total potential, Φ, on the right-hand side of Equations (1) and (2) is contributed by both luminous and dark matter. As we have mentioned, the luminous matter distribution is directly inferred from the surface brightness of the dwarf galaxy (see Section 2.1 above). To model the density profile of dark matter, we adopt in our analysis either the Navarro–Frenk–White (NFW) profile or a double power-law functional form of

$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{s}}{{\left(r/{r}_{s}\right)}^{\gamma }{\left(1+r/{r}_{s}\right)}^{\alpha }}\end{eqnarray} \tag{ 5 }$

with the model parameters (ρ_s, r_s, γ, and α) to be constrained. Note in our analysis throughout this paper, the outer power-law index, α, will be fixed to 3.

In order to have analytical solutions for any given potential model and tracer distribution, MGE is not only applied to the two-dimensional surface density distribution of the luminous stellar component (see Section 2.1.4 above), but also to the underlying model for the dark matter distribution and to the density distribution of tracers (ν) as well.¹⁹ Each MGE component would have analytical solutions to Equations (1) and (2). In principle, each MGE component of the tracer population can have its own rotation parameter, κ_k, and velocity anisotropy parameter, b_k. The M/L for each MGE component can also differ, but in our analysis, we treat κ, b, and M/L to be the same for different MGEs.

For an observed star with position ${{\boldsymbol{x}}}_{i}^{{\prime} }=(x{{\prime} }_{i},y{{\prime} }_{i})$ on the image plane, which has observed velocity ${{\boldsymbol{v}}}_{i}=\left({v}_{x^{\prime} ,i},{v}_{y^{\prime} ,i},{v}_{z^{\prime} ,i}\right)$ and an error matrix of

$\begin{eqnarray}{{\boldsymbol{S}}}_{i}=\left(\begin{array}{ccc}{\sigma }_{{v}_{x^{\prime} },i}^{2} & 0 & 0\\ 0 & {\sigma }_{{v}_{y^{\prime} },i}^{2} & 0\\ 0 & 0 & {\sigma }_{{v}_{z^{\prime} },i}^{2}\end{array}\right),\end{eqnarray} \tag{ 6 }$

its position, ${{\boldsymbol{x}}}_{i}^{{\prime} }$ , can be transformed to the intrinsic frame to solve the corresponding velocities and velocity dispersions, based on a set of model parameters, Θ. A solution for each MGE is sought, and the solutions of different MGEs are added together in the end. The solutions are then transformed back to the observing frame. The mean velocity predicted by the model in the observing frame is denoted as ${{\boldsymbol{\mu }}}_{i}=\left({v}_{x^{\prime} ,i},{v}_{y^{\prime} ,i},{v}_{z^{\prime} ,i}\right)$ , and the covariance matrix is defined through both the first and the second velocity moments,

$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{C}}}_{i}\\ =\,\left(\begin{array}{ccc}\overline{{v}_{{\rm{x}}^{\prime} ,{\rm{i}}}^{2}}-{\overline{{v}_{{\rm{x}}^{\prime} ,{\rm{i}}}}}^{2} & \overline{{v}_{{\rm{x}}^{\prime} {\rm{y}}^{\prime} ,{\rm{i}}}^{2}}-\overline{{v}_{{\rm{x}}^{\prime} ,{\rm{i}}}}\,\overline{{v}_{{\rm{y}}^{\prime} ,{\rm{i}}}} & \overline{{v}_{{\rm{x}}^{\prime} {\rm{z}}^{\prime} ,{\rm{i}}}^{2}}-\overline{{v}_{{\rm{x}}^{\prime} ,{\rm{i}}}}\,\overline{{v}_{{\rm{z}}^{\prime} ,{\rm{i}}}}\\ \overline{{v}_{{\rm{x}}^{\prime} {\rm{y}}^{\prime} ,{\rm{i}}}^{2}}-\overline{{v}_{{\rm{x}}^{\prime} ,{\rm{i}}}}\,\overline{{v}_{{\rm{y}}^{\prime} ,{\rm{i}}}} & \overline{{v}_{{\rm{y}}^{\prime} ,{\rm{i}}}^{2}}-{\overline{{v}_{{\rm{y}}^{\prime} ,{\rm{i}}}}}^{2} & \overline{{v}_{{\rm{y}}^{\prime} {\rm{z}}^{\prime} ,{\rm{i}}}^{2}}-\overline{{v}_{{\rm{y}}^{\prime} ,{\rm{i}}}}\,\overline{{v}_{{\rm{z}}^{\prime} ,{\rm{i}}}}\\ \overline{{v}_{{\rm{x}}^{\prime} {\rm{z}}^{\prime} ,{\rm{i}}}^{2}}-\overline{{v}_{{\rm{x}}^{\prime} ,{\rm{i}}}}\,\overline{{v}_{{\rm{z}}^{\prime} ,{\rm{i}}}} & \overline{{v}_{{\rm{y}}^{\prime} {\rm{z}}^{\prime} ,{\rm{i}}}^{2}}-\overline{{v}_{{\rm{y}}^{\prime} ,{\rm{i}}}}\,\overline{{v}_{{\rm{z}}^{\prime} ,{\rm{i}}}} & \overline{{v}_{{\rm{z}}^{\prime} ,{\rm{i}}}^{2}}-{\overline{{v}_{{\rm{z}}^{\prime} ,{\rm{i}}}}}^{2}\end{array}\right).\end{array}\end{eqnarray} \tag{ 7 }$

By assuming the velocity distribution predicted by the model is a trivariate Gaussian with mean velocity μ _i and covariance C _i at ${{\boldsymbol{x}}}_{i}^{{\prime} }$ , the likelihood can be written as

$\begin{eqnarray}\begin{array}{rcl}{L}_{i}^{\mathrm{dwarf}} & = & p\left({{\boldsymbol{v}}}_{i}| {{\boldsymbol{x}}}_{i}^{{\prime} },{{\boldsymbol{S}}}_{i},{\rm{\Theta }}\right)\\ & = & p\left({{\boldsymbol{v}}}_{i}| {{\boldsymbol{x}}}_{i}^{{\prime} },{{\boldsymbol{S}}}_{i},{{\boldsymbol{\mu }}}_{i},{{\boldsymbol{C}}}_{i}\right)\\ & = & \displaystyle \frac{\exp \left[-\tfrac{1}{2}{\left({{\boldsymbol{v}}}_{i}-{{\boldsymbol{\mu }}}_{i}\right)}^{{\rm{T}}}{\left({{\boldsymbol{C}}}_{i}+{{\boldsymbol{S}}}_{i}\right)}^{-1}\left({{\boldsymbol{v}}}_{i}-{{\boldsymbol{\mu }}}_{i}\right)\right]}{\sqrt{{\left(2\pi \right)}^{3}\left|\left({{\boldsymbol{C}}}_{i}+{{\boldsymbol{S}}}_{i}\right)\right|}}.\end{array}\end{eqnarray} \tag{ 8 }$

In addition to the above model for the dwarf galaxy itself, the discrete jam model can also model a population of foreground/background stars. The surface density of foreground/background stars is assumed to be constant throughout the field of the dwarf, and is modeled to be a given fraction, , of the central surface density of the dwarf, Σ(0, 0). With the assumptions, the prior on dwarf membership can be written as

$\begin{eqnarray}&&{m}_{i}({{\boldsymbol{x}}}_{i})=\displaystyle \frac{{\rm{\Sigma }}({{\boldsymbol{x}}}_{i})}{{\rm{\Sigma }}({{\boldsymbol{x}}}_{i})+\epsilon {\rm{\Sigma }}(0,0)}.\end{eqnarray} \tag{ 9 }$

The velocity distribution of the foreground/background stars can also be modeled as a trivariate Gaussian with the given mean velocity and velocity dispersions. Then the likelihood of the dwarf and foreground/background can be combined through

$\begin{eqnarray}&&{L}_{i}={m}_{i}({{\boldsymbol{x}}}_{i}){L}_{i}^{\mathrm{dwarf}}+[1-{m}_{i}({{\boldsymbol{x}}}_{i})]{L}_{i}^{\mathrm{bkgd}}.\end{eqnarray} \tag{ 10 }$

The total likelihood is the product of the likelihood for each star

$\begin{eqnarray}&&L=\prod _{i=1}^{{N}_{\mathrm{star}}}{L}_{i}.\end{eqnarray} \tag{ 11 }$

The parameters used in our modeling are the following:

1.
Rotation parameter, κ;
2.
Velocity anisotropy, b;
3.
Stellar-mass-to-light ratio, M/L;
4.
Dark matter halo scale density, ρ_s;
5.
Dark matter halo scale radius, r_s;
6.
Inner density slope of the host dark matter halo, γ;
7.
The background fraction, .

In particular, instead of directly fitting ρ_s and r_s, we fit ${d}_{1}={\mathrm{log}}_{10}{\rho }_{s}^{2}{r}_{s}^{3}$ and ${d}_{2}={\mathrm{log}}_{10}{\rho }_{s}$ . This is to alleviate the degeneracy between the halo parameters and also make the parameter space to be more efficiently explored in log space (Zhu et al. 2016a). In our analysis, the M/L will mostly be fixed to unity, i.e., its true value, except for some cases which will be otherwise specified. The inner density slopes are also sometimes fixed to 1, i.e., for the NFW model profile. For foreground/background contamination-free cases, we fix to 0. Throughout this paper, we fix the distance and the inclination angle to their true values. However, we explicitly tested the results when the distance is allowed to be a free parameter. In this case, the best-constrained mass profiles remain very similar, while the best-fitting second velocity moments in the very inner regions can be slightly improved, though the best-fitting distance might slightly differ from the true distance by up to ∼10%. Moreover, the outer density slope, α, will be fixed to 3, but we have also tried to vary the outer slopes and add a constant background density, and our conclusions are not affected.

4. Results

4.1. Overall Performance

The left plot of Figure 2 shows the best-fitting versus true masses for all 28 dwarfs. We show this for the mass enclosed within the half-mass radius of tracers, M(< r_half), and the mass between 200 and 300 pc, M(200–300 pc). The fitting is based on the double power-law model profile for the underlying dark matter (Equation (5)) with free inner slopes. Note that the matter density within 150 pc is often used as the proxy to the inner slopes (e.g., Read et al. 2019). Unfortunately, 150 pc is below the resolution limit of auriga. So we focus on the radial ranges that are above the resolution limit for now, and we will investigate the constraint within 150 pc later in Section 5.1, subject to extrapolations to the very center.

In general, M(< r_half) is better constrained, with a scatter of ∼0.067 dex. This is in good agreement with the commonly accepted experience that the mass enclosed within the half-mass (or half-light) radius of tracers can be more robustly constrained upon dynamical modeling (e.g., Wolf et al. 2010; Walker & Peñarrubia 2011; Wang et al. 2015, 2020). This is due to the degeneracy between the two halo parameters, ρ_s and r_s, reflecting that dynamical models mostly constrain the gravity or mass at the median radius of the tracers, but are less sensitive to the shape of the mass profile (e.g., Han et al. 2016a; Li et al. 2021, 2022a). There are a few cases for which M(< r_half) is constrained slightly worse than M(200–300 pc), which is quite rare though.

M(200–300 pc) is recovered reasonably, with a scatter of 0.167 dex, which is larger than the scatter of M(< r_half). We further classify systems into different groups according to their star formation activities. At first, we notice that there are 11 star-forming systems with prominent galactic winds or gas outflows in the simulation, and on top of the corresponding red circles, we overplot green stars. An example of such outflows is shown in the left plot of Figure 15, which we will discuss later. There are another five systems with larger than zero SFRs, but they do not show obvious galactic winds, which we denote with blue diamonds. The six Sag-like systems are denoted by cyan squares. Interestingly, systems classified in these groups show different trends. Quiescent Sag-like systems are more likely to have M(200–300 pc) underestimated, whereas the best-fitting M(200–300 pc) for star-forming systems with strong outflows are biased to be more above the black dashed line by ∼0.05 dex on average. The other systems are more symmetrically distributed on both sides of the black dashed line.

In the right plot of Figure 2, it is shown that the bias in best-fitting M(200–300 pc) from the truth correlates with the sSFR. The correlation is weak or absent for M(< r_half). For the five systems with zero SFRs but without prominent galactic winds, their sSFRs are weak, while the biases in their best fits are smaller on average. The six Sag-like systems are all quiescent with zero sSFR. Those systems with zero sSFR but are not too close to the Galaxy Center and are not yet undergoing strong tidal effects, i.e., not Sag-like, tend to show smaller biases in their best fits as well (red circles without star, diamond, or square symbols).

In order to understand the cause for such a correlation between the bias in best fits and the sSFR, in Figure 3 we provide the best-fitting versus true masses between 200 and 300 pc in log space, versus the dynamical status of the systems (x-axis) between 200 and 400 pc. In the x-axis of Figure 3, θ is the radial action angle, or referred to as the phase angle by Han et al. (2016a, 2016b). For a star at radius r, θ is defined as

$\begin{eqnarray}&&\theta (r)=\displaystyle \frac{1}{T}{\int }_{{r}_{p}}^{r}\displaystyle \frac{{dr}^{\prime} }{{v}_{r}},\end{eqnarray} \tag{ 12 }$

where T is the radial period, r_p is the radius at pericenter, and v_r is the radial velocity. Note that when radial cuts are adopted in the tracer sample, the integral boundaries have to be properly changed accordingly. Equation (12) can be easily evaluated for a spherical potential. Here, we first calculate potential profiles for each dwarf based on the true radial distribution of all particles (star, dark matter, gas, and black hole) under the spherical approximation, and we evaluate θ for each tracer star particle using the potential profile. Since our systems are not strongly deviating from spherical symmetry, the spherical assumption should approximately hold for most of our dwarfs, though not all. For a system in steady state, θ evaluated from the true potential is expected to distribute uniformly between 0 and 1, with the median value close to 0.5. The deviation of the medians from 0.5 thus indicates the number of deviations from steady states. This is why the x-axis quantity in Figure 3 is chosen as the median value of the θ for each system, which we denote as θ_med. The larger θ_med deviates from 0.5 (the black vertical dashed line), and the stronger the system is deviating from the steady state.

We adopt particles between 200 and 400 pc to calculate the phase angles. The bias of best-fitting versus true M(200–300 pc) correlates with the dynamical status over this range.²⁰ The scatter is large, but there exists the trend that systems more strongly deviating from steady states between 200 and 400 pc are on average more likely to have larger biases in M(200–300 pc). Explicitly, we find systems with underestimated M(200–300 pc) mostly tend to have θ_med larger than 0.5. In other words, their phase angle distribution is biased to large values. On the contrary, many systems with overestimated M(200–300 pc) tend to have θ_med smaller than 0.5, i.e., phase angle distribution is biased to small values, though some of them still have θ_med > ∼0.5. Note for these systems with overestimated M(200–300 pc) and θ_med > 0.5, we find part of them look a bit flattened, such as the blue diamond in the upper right corner (Au27-2), so the calculation of their θ_med might be slightly affected by the fact that we adopted spherical potential profiles to calculate their phase angles, but the overall trend is reasonable.

Nevertheless, Figure 3 unambiguously indicates that the amount of bias in best fits is related to the dynamical status of the systems, and thus the correlation between the bias in best fits and the sSFR reflects the dynamical status behind.²¹ For systems with strong galactic winds, despite the fact that wind particles are not used as tracers in our analysis, the galactic winds have caused deviations from steady states, resulting in larger biases (mostly overestimates in M(200–300 pc) in their best fits). For the other five systems with larger than zero but weak star formations, perhaps because they do not present prominent outflows to disturb the system, they tend to be closer to steady states on average and also have smaller biases in the best fits.

On the other hand and very interestingly, we find many of the Sag-like systems are undergoing some level of contractions along the longer image axis²² , which perhaps results in the underestimated M(200–300 pc). These include Au16-9, Au23-7, and Au24-24, as marked in the figure with corresponding numbers. An example of the contraction is shown in the left panel of Figure 14 for Au16-9, which we will discuss later. In addition to the three systems, the other three Sag-like systems in Table 1 (Au21-10, Au23-4, and Au27-25) show less significant underestimates in M(200–300 pc). We only see some very weak contractions along the $x^{\prime}$ -axis for Au21-10 and Au23-4, and this is perhaps why the best fits for the two dwarfs agree better with the truth. Au27-25 undergoes strong expansions along the $y^{\prime}$ - and $z^{\prime}$ -axes, and along the $x^{\prime}$ -axis, it shows contractions. Note the expansion is not due to galactic winds as all the Sag-like systems have zero SFR and no wind particles.

The contractions and expansions in these Sag-like systems are due to the strong tidal effects at small distances. In fact, all six Sag-like systems are undergoing some expansions along the line of sight. This is expected because the tidal disruption occurs along the radial direction pointing to the Galaxy Center, and although the artificial observer is 8 kpc away from the Galaxy Center, the line-of-sight direction still aligns with the radial direction. As shown by Ogiya et al. (2022), when the dwarf is close to its pericentric passage, the effective radius first undergoes a small decrease, followed by a large increase, sometimes associated with oscillations. On the contrary, if without tidal forces, Ogiya et al. (2022) showed that the effective radius stays almost constant.

Our results seem to indicate strong tidal effects, which can cause both contractions and expansions, tend to cause underestimates in M(200–300 pc), and perhaps it is the contraction that is more related to the underestimated M(200–300 pc). Perhaps the coherent contractions have caused smaller velocity dispersions in the inner regions and hence underestimated the inner profiles. Besides, due to the global contraction, the apocenters of tracer orbits likely move inward with time, resulting in θ_med greater than 0.5 if θ is evaluated based on the current underlying potential. This is because the potential changes with time, but the apocenters will be overestimated if they are calculated from the current potential by assuming that the potential does not evolve (steady state). On the other hand, outward motions or expansions, such as gas outflows or galactic winds, are more likely to cause overestimates in M(200–300 pc). The galactic winds cause lowered density in central regions, driving the whole system out of equilibrium. When θ is evaluated based on the current potential, assuming the system is in a steady state, the distribution of θ is biased to be smaller, with θ_med less than 0.5.

In the next subsection, we perform more detailed investigations on the best fits for the six quiescent Sag-like systems and another four star-forming systems with strong outflows.

4.2. Examples of the Best-fitting Matter Density Profiles

In this subsection, we further investigate the best-fitting and true density profiles for six Sag-like systems and another four systems with prominent galactic winds.

4.2.1. Best-fitting Density Profiles and Velocity Maps for Sag-like Systems

We start with the six Sag-like systems. In Figures 2 and 3, we have shown that the best-fitting masses between 200 and 300 pc, M(200–300 pc), are underestimated for the six systems and to different levels, when using the double power-law model profile for the underlying dark matter. In Figure 4, the best fits (magenta curves) based on the double power-law potential model and true total mass profiles (green dots with error bars) are shown for the six systems. For Au16-9, Au21-10, Au23-7, and Au24-24, their best-fitting inner profiles are lower than the truth, i.e., the magenta solid curves are below the green dots within ∼0.5 kpc, and the differences are significant compared with the small error bars. M(200–300 pc) of Au23-4 and Au27-25 are mildly underestimated, so the differences between the magenta curves and the green dots are not very obvious in the corresponding panels.

The M/L is fixed to unity for the magenta curves, but the underestimates in the inner profiles remain the same, even if we allow the M/L to be free. Besides, we have fixed the outer slopes to be α = 3. However, the poor constraints in the outskirts are not improved if we allow the outer slopes to be free. This is probably due to the significantly smaller number of tracer star particles in the outskirts. These systems mostly have overestimated outer profiles and underestimated inner profiles. The true and best-fitting profiles cross each other at the radii close to the half-mass radii. Note although only bound particles are used as tracers, the true profiles in Figure 4 are calculated based on all particles at the corresponding radial ranges.²³ This is because though particles belonging to the host halo might be approximated to have a constant density over the scale of the dwarf, they still contribute to the potential gradient. Since background particles belonging to the host halo could be more dominant for such nearby Sag-like systems, we have also tried to model the underlying matter distribution using the double power-law model plus a constant-density model. However, the best-constrained model still tends to underestimate the inner profiles and overestimate the outer profiles.

We also show in Figure 4 the best-fitting star and dark matter profiles separately. The red, black, and green curves denote the best-fitting star, dark matter, and total profiles based on the NFW model profile. Here, the M/L is allowed to be free. It is clearly shown that the profiles of the stellar component, and hence the M/L, are poorly constrained by jam, with the best fits (red curves) significantly deviating from the truths (red diamonds) in most cases. The dark matter profiles, on the other hand, can be constrained better, but still show some significant deviations from the truth in inner regions, especially for Au23-4 and Au23-7. Fortunately, despite the relatively worse fits to the stellar and dark matter components individually, the total profiles can be fit significantly better.

Due to the poorly constrained M/L, we have fixed the M/L to unity in Figures 2 and 3, when the double power-law potential model is used. However, with the double power-law model, we have tried to allow the M/L to be free as well, though we choose not to explicitly show the results to avoid redundancy. Our conclusions remain very similar. The stellar and dark matter profiles are poorly constrained individually. Actually, with the free inner slopes in the double power-law model, the fitting to the dark matter component can become significantly lower than the truth in inner regions for a few cases, whereas the stellar component is significantly overestimated, with the total profiles (stellar + dark matter) still more reasonably constrained. Hereafter, unless otherwise specified, the M/L will be fixed to unity, i.e., its true value. This is a reasonable operation because observationally, the M/L can be alternatively constrained through stellar population synthesis modeling and fixed upon dynamical modeling.

The poor fits to stellar and dark matter profiles individually are mainly due to the degeneracy between the M/L and halo parameters. First of all, the true stellar and dark matter profiles share similar shapes as shown in Figure 4. Note that all particles have been included in these true profiles, and if only using bound particles, the stellar and dark matter profiles would appear more similar in the outskirts. This is due to the strong tidal stripping that happened to these Sag-like systems, making them have similar outer dark matter and stellar halo profiles (e.g., Errani et al. 2022). As a result, the model has difficulties distinguishing the stellar and dark matter components, as they have similar contributions to the shape of the underlying potential profiles. This explains why the best fits to stellar and dark matter components are poor individually. Moreover, we show in Figure 5 the likelihood contours of different parameter combinations for Au16-9. The anticorrelations or degeneracies among the M/L and the two halo parameters (d₁ and d₂) are prominent.

**Figure 5.** The likelihood contours of all possible combinations of every two out of the five model parameters for Au16-9 (NFW potential model with free M/L and 6000 star particles as tracers). The yellow contour includes 10% of the MCMC sample centered on the best fits. The light blue and two dark blue contours represent the 1σ, 2σ, and 3σ confidence levels. The degeneracies between the tracer parameters (κ and b) and the potential parameters (M/L, ${d}_{1}={\mathrm{log}}_{10}{\rho }_{s}^{2}{r}_{s}^{3}$ , and ${d}_{2}={\mathrm{log}}_{10}{\rho }_{s}$ ) are very weak, but there are stronger degeneracies among M/L, d₁, and d₂.
Download figure:
Standard image High-resolution image

**Figure 5.** The likelihood contours of all possible combinations of every two out of the five model parameters for Au16-9 (NFW potential model with free M/L and 6000 star particles as tracers). The yellow contour includes 10% of the MCMC sample centered on the best fits. The light blue and two dark blue contours represent the 1σ, 2σ, and 3σ confidence levels. The degeneracies between the tracer parameters (κ and b) and the potential parameters (M/L, ${d}_{1}={\mathrm{log}}_{10}{\rho }_{s}^{2}{r}_{s}^{3}$ , and ${d}_{2}={\mathrm{log}}_{10}{\rho }_{s}$ ) are very weak, but there are stronger degeneracies among M/L, d₁, and d₂.
Download figure:
Standard image High-resolution image

Comparing the magenta and green curves in Figure 4, we find that the inner profiles are more poorly constrained when the double power-law model profile is used for Au16-9, Au21-10, Au23-7, and Au24-24 than when the NFW profile is used. To explore why the inner profiles are more poorly constrained when the inner slopes are allowed to free, we check the best-fitting and true velocity maps in Figure 6, which shows the mean velocity and velocity dispersion profiles for the $x^{\prime}$ component and for Au16-9 as an example. Both the velocity and velocity dispersion profiles have large scatters, and the double power-law model tends to go through the middle of the velocity dispersion profiles in the inner regions, while the NFW model leads to a curve that is higher in the center. Both are reasonable fits, but in fact, the double power-law model achieves a better fit with a log-likelihood larger than that of the NFW model by ∼33, which is significantly larger than the expected 1σ uncertainty of the five free parameters. This in fact benefited from our usage of a large sample of 6000 member stars. While the velocity profiles are better fit with the double power-law model profile, the deviation of the best-fitting double power-law model from the true density profile is larger, which reflects the deviation from the steady state due to strong tidal effects, and thus a better fit to the velocity maps lead to a worse fit to the density profile. Explicitly, as we have mentioned, the strong contraction of member stars in Au16-9 (see Figure 14) has likely caused smaller velocity dispersions in the inner regions, and thus the double power-law model with underestimated inner profiles better fits the velocity map. By using the NFW model, more masses are forced to be in the inner region, leading to a better match to the true density profile but a worse fit to the velocity map.

**Figure 6.** The $x^{\prime}$ component of the mean velocity (left) and velocity dispersion (right) profiles of star particles, binned along the major (green) and minor (black) axes for Au16-9. Note that along the major and minor axes, only stars within sectors of ±45° to the corresponding axes are used to calculate the mean velocity and velocity dispersion in each bin. Each bin contains 100 stars. The x and y errors indicate the bin width and the 1σ scatters, respectively. Green and black solid curves denote the model predictions along the major and minor axes, based on the NFW model profile. Blue and magenta solid curves denote the model predictions along the major and minor axes, based on the double power-law model profile for the dark matter with free inner slopes. Dashed curves around the solid ones with corresponding colors indicate the 1σ uncertainties in the model. The green, black, blue, and magenta curves in the left panel are all very close to the horizontal lines at y = 0, and hence are overlapping each other.
Download figure:
Standard image High-resolution image

4.2.2. Best-fitting Density Profiles and Velocity Maps for Star-forming Systems with Outflows

Now we discuss a few examples of star-forming systems with prominent outflows. Figure 7 shows the true and best-fitting total profiles for four representative star-forming systems. They all have strong outflows and large biases in their best-fitting M_fit(200–300 pc). The four systems are marked with corresponding numbers in Figure 3. For simplicity, we do not show the profiles for the stellar and dark matter components, but only show the total profiles. These systems all have overestimated inner density profiles and underestimated outer profiles, with the true and best-fitting profiles crossing at roughly the half-mass radii.

The mean velocity and velocity dispersion profiles of star particles for Au21-9 are shown in Figure 8, and we only show the $x^{\prime}$ component as an example. The apparent sizes of such distant star-forming dwarfs are smaller than the nearby Sag-like objects shown in Figure 6. Overall, the agreement between the truth and the best fit is reasonable. The model tends to fit a rotation, i.e., ${v}_{x}^{\prime}$ positive along the positive $y^{\prime}$ -axis (minor axis) and negative along the negative $y^{\prime}$ -axis. These systems have strong galactic winds (wind particles), but the star particles do not show obvious signatures of outward motions. However, there are some indications of inflows of the stars on large scales along the major axis, that is, ${v}_{x}^{\prime}$ is negative along the positive $x^{\prime}$ -axis (major axis) and positive along the negative $x^{\prime}$ -axis. Such inflows are perhaps due to tidal effects because this dwarf is on its way to approaching the host Galaxy Center. Besides, the stellar system can expand in response to the lowered potential due to gas outflows at the beginning, which undergoes some overshooting and as a result would contract and get stabilized later (e.g., Navarro et al. 1996; Li et al. 2022b). Maybe we happened to observe this dwarf at the contraction stage. The best-fitting velocity dispersion along the major axis is slightly lower than the truth beyond 40', while the best-fitting dispersion along the minor axis seems to be a bit higher in amplitude than the truth beyond 10', which are perhaps related to the deviation from steady states in the outskirts.

Comparing Figures 4 and 7, it seems both figures have the best-fitting and true total profiles crossing at approximately the half-mass radii. As we have mentioned, due to the degeneracy between the halo parameters, most dynamical models mainly constrain the mass at the median radius of the tracers, but are less sensitive to the shape of the mass profile. To maintain the relatively robust constraints on the mass at the half-mass radius of tracers, the biases in the best-fitting inner and outer density profiles for the same system often happen in opposite ways. In addition, quiescent Sag-like systems and star-forming systems also tend to have opposite trends. Sag-like systems tend to have underestimated inner profiles and overestimated outer profiles, whereas star-forming systems with prominent outflows tend to have overestimated inner profiles and underestimated outer profiles. This is because, as we have explained in the previous subsection, the contractions/expansions (perhaps mainly the contractions) due to strong tidal effects for nearby Sag-like systems and galactic winds/gas outflows for more distant star-forming systems drive the systems out of equilibrium, but in opposite ways.

4.3. Constraints Based on Fewer Tracers and Only Line-of-sight Velocities with Observational Errors

Our results in the previous sections are based on a large sample of 6000 star particles as tracers, with both line-of-sight velocities and proper motions, but without errors. The error-free and large sample of tracers has enabled us to control the size of statistical errors while focusing on investigating the intrinsic systematics behind jam. However, current observations of nearby dwarf galaxies at most provide kinematical information on ∼2000 member stars (e.g., Zhu et al. 2016a), most of which do not have accurate proper motion measurements. Thus, for more practical meanings, in this subsection, we repeat our analysis by using a smaller sample of 2000 star particles as tracers, and we only use their line-of-sight velocities after including observational errors (see Section 2.1 for details). After excluding distant systems that cannot have more than 2000 stars brighter than g = 21, we end up with 16 systems.

The best-constrained M_fit(200–300 pc) versus their true values is shown in Figure 9. We overplot the error-free constraints based on 6000 star particles as black empty triangles. Approximately, it seems M(200–300 pc) can still be an unbiasedly constrained ensemble, but compared with the measurements based on 6000 star particles and both line-of-sight velocities and proper motions, the measurements are now slightly biased more above the diagonal reference line at M(200–300 pc) > 10⁷ M_⊙. The scatter is increased by ∼50%, with the black triangles denoting a scatter of 0.169 dex²⁴ and the red circles denoting a scatter of 0.260 dex. Our results thus indicate that systematic uncertainties caused by deviations from the steady states of the dwarf systems are dominating the final errors of mass estimation. If without accurate proper motion measurements, line-of-sight velocities of ∼2000 tracers with typical observational errors can still achieve reasonable and approximately unbiased ensemble measurements for M(200–300 pc) on average, but with slightly larger scatters.

4.4. Constraints on Sag-like Systems after Including Observational Errors and Unbound Star Contamination

In the previous subsections, our results have been based on bound star particles as tracers. In this section, we discuss the results after incorporating line-of-sight velocity errors, Gaia DR3 parallax errors, and either Gaia DR3 or CSST (Zhan 2011; Cao et al. 2018; Gong et al. 2019) proper motions errors. We also select member stars according to their kinematics. We first present results based on Gaia DR3 errors for the six Sag-like systems. We then show results based on the CSST proper motion errors for the dwarf with the closest distance to the mock observer (Au16-9).

After including typical line-of-sight velocity errors, Gaia DR3 proper motions, and parallax errors, we select tracer star particles based on their differences in proper motions, line-of-sight velocities, and parallaxes with respect to the dwarf centers (see Section 2.2 for details), and the best fits are shown as blue dashed curves in Figure 10 for the six Sag-like systems. Stars selected in this way not only include those true bound member stars, but may also include foreground/background contamination and some unbound star particles in the vicinity of the dwarf. As we have explicitly checked, because the associated observational errors are relatively small for such nearby Sag-like systems, the fraction of actual foreground or background stars is very small. In addition to bound member stars, most of the other tracer stars selected in this way are those unbound star particles surrounding the dwarf system. Part of these star particles might have once been bound to the dwarf, and have been stripped by the host halo recently.

**Figure 10.** Similar to Figure 4. Green dots with error bars denote the true total density profiles of the six Sag-like systems, and magenta solid curves denote the best fits based on the double power-law functional form for the underlying potential (errors not modeled). These are exactly the same as in Figure 4. We further include realistic observational errors (1–10 km s⁻¹ of errors in line-of-sight velocities and Gaia DR3 proper motion and parallax errors), and the tracer stars are selected according to their difference in kinematics with respect to the dwarf centers (see Section 2.2 for details). Blue dashed curves denote the best fits based on tracer star particles selected after including the observational errors. Cyan and yellow dashed curves denote best fits based only on line-of-sight velocities and only on proper motions, respectively, with observational errors included.
Download figure:
Standard image High-resolution image

After including observational errors and contamination by unbound stars, the uncertainties as indicated by the blue-shaded regions are now significantly larger than the original magenta-shaded regions, but the best fits stay very similar to the original error-free best fits. On the basis of the blue dashed curves, the yellow and cyan dashed curves denote the best fits using only proper motion data or only using line-of-sight velocities, respectively. For Au16-9, Au21-10, Au23-4 Au24-24, and Au27-25, the results based on pure proper motions and pure line-of-sight velocities are consistent with each other within 1σ. However, for Au23-7, pure proper motions and pure line-of-sight velocities show more significant differences. Interestingly, the yellow curves, i.e., with pure proper motions, show much better agreement with the true profiles. The larger difference between the yellow and cyan curves for Au23-7 is mainly related to how the tracer stars are selected. We found those unbound star particles selected as tracers and mainly in the outskirts of these systems show larger velocity dispersions in their line-of-sight velocities than those of bound particles. The line-of-sight velocity distribution of these particles also deviates from the assumed Gaussian distribution by jam, and jam does not model them as background. As a result, the increased velocity dispersions lead to the much higher amplitudes in the outskirts of the best fits, when both line-of-sight velocities and proper motions are used and when only line-of-sight velocities are used, but this is not the case when only proper motions are used.

Our results indicate that, for such nearby Sag-like systems, we can observe a sufficient number of bright stars, and with the accuracy of Gaia DR3 errors, proper motions can be as useful as or even sometimes better than line-of-sight velocities. For example, for Au16-9 at slightly beyond 20 kpc, the typical Gaia DR3 proper motion error of 0.1 mas yr⁻¹ corresponds to the error in tangential velocity of ∼5 km s⁻¹, which is comparable to the typical error of 1–10 km s⁻¹ for the line-of-sight velocity. We have repeated our analysis by incorporating observational errors for more distant dwarf galaxies. Unfortunately, for systems close to or beyond 100 kpc, the number of bright stars is limited, and Gaia DR3 proper motion errors are not good enough to achieve reasonable constraints, with the fitting based on pure proper motion data very difficult to converge. This is expected because most stars with 18 < g < 19 have typical Gaia DR3 errors of 0.1 mas yr⁻¹, which corresponds to ∼40 km s⁻¹ of errors in tangential velocities at ∼80 kpc, and thus are useless.

However, the precision in future Gaia proper motion data is expected to be improved in a way proportional to t^1.5 (t is the time), indicating more than a factor of 2 reduction in error for DR4 compared with DR3, and additionally, a factor of ∼2.8 in improvement for the extended mission compared with the end of the 5 yr mission at fixed magnitude, pushing to a limit that can be even better than the typical line-of-sight velocity errors (e.g., Mateu et al. 2017; Brown 2021). This means that based on pure proper motion data, dynamical modelings, which presently can only be achieved for nearby Sag-like systems at the farthest distance of ∼50–60 kpc, now hold the promise of being applied to a sufficient number of fainter tracers in dwarf systems at or beyond ∼100 kpc in the future.

We now investigate the expected performance of using future CSST proper motions. The expected upper limit in the CSST proper motion error (∼0.2 mas yr⁻¹) is about a factor of 2 larger than that of Gaia DR3 at 18 < g < 19. For very nearby systems such as Au16-9 at a distance a bit more than 20 kpc, the corresponding tangential velocity error is ∼10 km s⁻¹, which quickly increases to ∼40 km s⁻¹ beyond 40 kpc. Thus, it seems dynamical modeling with pure proper motion data is difficult with 0.2 mas yr⁻¹ of error for systems more distant than ∼30−40 kpc. In Figure 11, we show the constraint based on the typical CSST proper motion error for Au16-9 as an example. The uncertainties are significantly larger, but the usage of pure proper motion data still leads to a reasonable fit.

**Figure 11.** Similar to Figure 10, but shows the results after incorporating CSST proper motion errors for only Au16-9.
Download figure:
Standard image High-resolution image

The typical CSST proper motion error is based on about six astrometric measurements evenly distributed in 10 yr of baseline. Our results thus indicate that this particular observational strategy would provide useful proper motions for dynamical modeling of very nearby systems at distances < ∼20 kpc. However, due to the fact that the CSST proper motion error estimates are based on the star number density in the Galactic bulge region, the actual errors based on the same operational mode can be smaller, so with pure CSST proper motion data, it could be possible to push slightly further. Notably, as will be shown in Nie et al., the associated proper motion errors increase slowly toward fainter magnitudes, reaching ∼0.4 mas yr⁻¹ at the CSST g-band apparent magnitude of 24, whereas Gaia proper motion errors increase very quickly beyond G ∼ 20. This is because the aperture and exposure time of CSST is designed to have higher signal-to-noise ratios at fainter magnitudes than Gaia, with the g-band limiting magnitude of about 26–27 (Zhan 2011). Thus, future CSST proper motion measurements hold the promise of compensating Gaia measurements at fainter magnitudes.

5. Discussion

5.1. Connection to the Core-cusp Problem

The softening scale of the level 3 auriga set of simulations is 185 pc at z = 0. For radial range below this scale, the density profiles and dynamics of tracer particles are not realistic, so it is difficult for us to directly reach the scale that is the most relevant to the core-cusp problem (150 pc; e.g., Read et al. 2019). However, based on our results in the previous section, the best-fitting inner density profiles show different amounts of biases for individual systems, which also show a dependence on the star formation activities. It seems similar biases would exist on scales below the resolution limit as well, and thus we should be cautious when interpreting the results for individual systems.

In Section 4.2, the best-fitting inner profiles of a few quiescent Sag-like systems are shown to be below the true profiles when the inner slopes are allowed to be free. An interesting case is Au16-9. For this system, both the NFW and the double power-law model profile lead to reasonable fits to the velocity map, but the qualitative goodness of fit is significantly better for the double power-law model. Despite of the better fit to the tracer kinematics, the actual inner density profile is significantly underestimated. Such inconsistency between the kinematical prediction and the actual density distribution indicates strong deviations from steady states. Violation of the steady-state assumption could result in more core-ish inner profiles when the truth is more cuspy for individual systems. Au16-9 is close to the Galactic Center and is undergoing some prominent contractions due to strong tidal effects, and we should avoid using systems that are strongly out of equilibrium, otherwise, we may end up drawing the wrong conclusions regarding the core-cusp issue, especially when looking at individual or a small sample of systems.

In real observation, instead of directly constraining the inner slopes, the best-fitting total mass within 150 pc, M(<150 pc), is often adopted to infer the steepness of the inner slope (e.g., Read et al. 2019). We now compare the best-fitting and true mass within 150 pc. To achieve the comparison, we extrapolate the actual density profiles in the simulation down to scales below the resolution limit, based on the double power-law model profile. The extrapolation is used to calculate M_true(<150 pc), which is regarded as the truth. On the other hand, the best-fitting mass within 150 pc, M_fit(<150 pc), is measured through tracer dynamics beyond 185 pc, extrapolated to the central region. Here we have excluded tracers within 185 pc to avoid wrong dynamics. Note in previous sections, our discussions are focused on the radial range not affected by star particles within 185 pc, and thus including or excluding tracers within 185 pc barely affects our results.

The results are shown as red circles in Figure 12. The associated uncertainties in the true masses are computed according to the boundaries which enclose 68% of the MCMC chains.²⁵ Because there are no reliable tracers within 185 pc, now the jam-constrained M_fit(<150 pc) through stellar dynamics has significantly larger error bars than those of M_fit(200–300 pc) and M_fit(< r_half) in previous plots. For convenience, we also overplot the black empty triangles which are the measurements for M(200–300 pc).

The scatter in M(< 150 pc) is 0.255 dex, which is larger than the scatter of M(200–300 pc) (0.167 dex), but the overall fitting is reasonable. Star-forming and quiescent systems are well separated, with quiescent systems having smaller M(<150 pc). We believe the amount of scatter of the red circles in Figure 12 is an upper limit because the scatter can be further decreased if tracers within 150 pc are available, and indeed in real observations stars within 150 pc can be observed and used for dynamical constraints.

We now check whether M(<150 pc) can be used as a good proxy to the inner slopes. Here, we calculate the true inner slopes, γ_true, based on the best-fitting parameter of the double power-law model profile to the actual total density profiles beyond 185 pc in the simulation, so γ_true depends on extrapolations to the very center. Figure 13 shows γ_true versus both M_fit(<150 pc) and M_true(<150 pc). The scatters and error bars are large, but there are correlations between γ_true and M_true(<150 pc) or M_fit(<150 pc). In the left plot, M_true(<150 pc) increases for steeper γ_true, indicating M(<150 pc) has the potential of being used as the proxy to infer the inner slopes. In the right plot, the correlation between γ_true and M_fit(<150 pc) seems to be weakened due to the scatter in the recovery of M(<150 pc). However, we emphasize that since we do not have tracers within 150 pc, the scatter in the right plot of Figure 13 and in Figure 12 is likely an upper limit. If tracers within 150 pc are available, the scatter is likely to be decreased.

**Figure 13.** The true inner slopes vs. the true (left) and best-fitting (right) masses enclosed within 150 pc. The true inner slopes and true masses within 150 pc are computed by extrapolating the actual density profiles down to the very center, using a double power-law function. The best-fitting mass within 150 pc is based on tracer dynamics beyond 185 pc. The meanings of the different symbols are the same as those in Figure 12.
Download figure:
Standard image High-resolution image

5.2. Detectability of Contractions and Outflows

In Section 4, we discussed the contractions of member stars due to strong tidal effects for nearby systems, and star-forming systems are also mentioned to have outflows. Such systems are strongly out of equilibrium, and thus caution has to be taken if using them for dynamical modeling. Hence, observational identification of such contractions and outflows is important. In this subsection, we briefly discuss the detectability of such motions.

To detect isotropic contractions or expansions, the usage of line-of-sight velocities of member stars or IFU observation might not be very useful because the effect would be averaged and canceled out along the line-of-sight direction. In addition, if the contraction mainly happens along the image longer axis but not along the line-of-sight direction as in our case, line-of-sight velocities will not be helpful either. Proper motions thus become more important.

In the left plot of Figure 14, we show the contraction of member stars of Au16-9, without including any observational errors. It is clearly shown that ${v}_{y}^{\prime}$ is negative/positive along the positive/negative $y^{\prime}$ -axis (minor axis), which demonstrates the contraction.²⁶ Note the minor axis of Au16-9 is in fact the longer image axis because our image major axis is defined according to the spin direction (see Section 3 for details), but the difference between the image major and minor axes of Au16-9 is very small anyway. Besides, ${v}_{y}^{\prime}$ along the image major axis tends to show some rotations, but the rotation is absent in Figure 6, indicating Au16-9 could be triaxial.

${v}_{y}^{\prime} $ — **Figure 14.** A demonstration of the detectability of contraction in nearby Sag-like systems. The figure is based on Au16-9, which is similar to the left plot of Figure 6 but shows the ${v}_{y}^{\prime}$ velocity component. Each bin contains 300 star particles. The left panel is based on bound member stars without including observational errors. Observational errors are incorporated in the middle panel. In the right panel, observational errors are included and tracers are selected according to their kinematics, which includes some unbound particles in the outskirts. We only show the $y^{\prime}$ component of the tangential velocities along the major and minor axes of the dwarf. The black horizontal dashed lines of y = 0 are to guide the eyes.
Download figure:
Standard image High-resolution image

In the middle plot of Figure 14, we include Gaia DR3 proper motion errors (see Section 2.2 for details), and the signal is noisier but still clearly detectable. In the right plot, not only the proper motion errors are included, but also the stellar tracers are selected according to their difference in kinematics with respect to the dwarf center. The inclusion of some unbound particles in the outskirts has caused some offsets toward negative y. This is because the global motion of the dwarf is defined through bound star particles and subtracted from all selected tracer particles. However, those unbound particles are not expected to have zero mean velocity after the subtraction. Despite the offset, the contraction is still detectable.

Our results indicate, promisingly, that current Gaia DR3 proper motion errors can enable us to detect such contraction motions for member stars of nearby Sag-like systems. The detection also benefited from the usage of a large sample of 6000 member stars.

In contrast to the detection of global contractions in nearby systems, the detection of gas outflows can be more challenging. The 11 dwarfs with gas outflows are at fairly large distances. The most nearby dwarf with such outflow is at 130 kpc. Figure 15 shows the major and minor axis proper motions of wind particles around this dwarf. In the left plot, no observational errors are incorporated, and we can clearly see the expansions, i.e., the $y^{\prime}$ tangential velocities are positive at the positive $y^{\prime}$ -axis and negative at the negative $y^{\prime}$ -axis (black dots with error bars). We now simply incorporate the Gaia DR3 proper motion errors to these wind particles in the right plot, treating them as individual stars. The trend becomes significantly noisy, but the signal is still marginally detectable for this system at 130 kpc. In a real observation, the detection of such gas outflows perhaps has to depend on observations of the hot gas distribution around these dwarfs at different epochs.

**Figure 15.** A demonstration of the detectability of outflows in a dwarf at a distance of ∼130 kpc. The figure is based on ∼4000 wind particles in Au23-2. Each bin contains 300 particles. Observational errors are included in the right plot, but not in the left plot. We only show the $y^{\prime}$ -component tangential velocities along the major and minor axes of the dwarf, and the trend is the same for the $x^{\prime}$ -component tangential velocities. The black dots show prominent expansions, while the green dots show some rotations. The black horizontal dashed line at y = 0 is to guide the eye.
Download figure:
Standard image High-resolution image

$y^{\prime} $ — **Figure 15.** A demonstration of the detectability of outflows in a dwarf at a distance of ∼130 kpc. The figure is based on ∼4000 wind particles in Au23-2. Each bin contains 300 particles. Observational errors are included in the right plot, but not in the left plot. We only show the $y^{\prime}$ -component tangential velocities along the major and minor axes of the dwarf, and the trend is the same for the $x^{\prime}$ -component tangential velocities. The black dots show prominent expansions, while the green dots show some rotations. The black horizontal dashed line at y = 0 is to guide the eye.
Download figure:
Standard image High-resolution image

6. Conclusions

In this study, we construct mock samples of member stars for 28 simulated dwarf galaxies from the cosmological level 3 set of auriga simulations. We applied the discrete axisymmetric jam modeling method to the 28 dwarf systems, to recover the underlying mass distribution through the stellar dynamics under the steady-state assumption. Among our sample of dwarf galaxies, there are six nearby Sag-like systems and 11 star-forming systems with strong galactic winds or gas outflows.

We first use 6000 bound member star particles without errors as stellar tracers for each system, which is to ensure good statistics. After assigning star particles apparent magnitudes using the multi-population synthesis code trilegal, we incorporate typical line-of-sight velocity errors (1−10 km s⁻¹), Gaia DR3 parallax errors, and Gaia DR3 or CSST proper motion errors. We also try a smaller sample of 2000 star particles as tracers and try to select stellar tracers according to their differences in kinematics to the dwarf centers.

Compared with the true matter distribution in the simulation, we find the recovered density profiles of the stellar and dark matter components individually are poor due to the degeneracy between stellar mass and dark matter distribution, and the constraints on the M/L are very weak. Fortunately, the density profiles of the total mass distribution can be constrained more reasonably.

The total mass within the half-mass radius of tracers can be constrained the best, with a scatter of ∼0.067 dex. The best-fitting masses between 200 and 300 pc, M_fit(200–300 pc), are an unbiasedly recovered ensemble, with a scatter of 0.167 dex. However, the amount of bias shows a dependence on the current sSFR. This is a result of the deviation from steady states. Most quiescent Sag-like systems undergoing strong tidal effects have their M_fit(200–300 pc) more underestimated than the truth, whereas most star-forming systems with outflows have M_fit(200–300 pc) overestimated.

The dependence of the best fits on sSFR reflects how tidal effects and galactic winds (gas outflows) perturb the system, driving it away from the steady state. We discover that most of the member stars in Sag-like systems have contractions along the image longer axis, which are the outcomes of strong tidal effects (e.g., Ogiya et al. 2022) and cannot be properly modeled by equilibrium models. On the other hand, the strong galactic winds or gas outflows in star-forming systems have driven the system out of equilibrium as well, despite the fact that wind particles themselves are not directly used as tracers. Tidal effects and galactic winds both cause the deviation from steady states, but in opposite ways, resulting in underestimated and overestimated inner density profiles, respectively. When talking about individual systems, one might end up with cored constraints on the inner profiles, whereas the truth is, on the contrary, more cuspy. On the other hand, cored systems might also be constrained as cuspy in the best fits for individual systems. Thus, one should avoid drawing strong conclusions based on the individual or limited number of systems. Instead, the averaged result based on a relatively large sample of dwarf galaxies is an unbiased ensemble and more robust.

Interestingly, systems that still have some amount of star formation, but without prominent outflows, mostly show smaller amounts of biases in their best fits. Besides, massive quiescent systems at large distances and are not yet undergoing strong tidal effects, i.e., not Sag-like, also show smaller amounts of biases on average. When using dynamical models based on the steady-state assumption, we should be cautious when using systems undergoing prominent contractions or expansions, which are strongly out of equilibrium. Promisingly, with Gaia DR3 proper motion errors, the contraction in nearby Sag-like systems is shown to be detectable.

Using a smaller sample of 2000 star particles with only line-of-sight velocities (typical errors from 1–10 km s⁻¹) as tracers, M(200–300 pc) can still be an unbiased constrained ensemble, but the scatter is increased by ∼50%.

After incorporating realistic observational errors (proper motion errors correspond to Gaia DR3) and selecting tracer stars based on the differences in their kinematics with respect to the dwarf center, we find their best-constrained density profiles agree with the error-free constraints within 1σ. For such nearby dwarf systems, dynamical constraints based on pure proper motions with Gaia DR3 errors perform equally well or even better than the case using pure line-of-sight velocities, which is very encouraging. On the other hand, the expected upper limit in proper motion errors for the future CSST survey is about twice as large as that of Gaia DR3, at CSST g-band apparent magnitudes of 18 < g < 19, which can already lead to reasonable constraints for very nearby systems at distances < ∼20 kpc.

By extrapolating the true density profiles in the simulation down to the radial range below the resolution limit, and by extrapolating the best-fitting density profiles based on tracer dynamics at r > 185 pc down to the center, we find the mass within 150 pc can be an unbiasedly recovered ensemble as well, but with a larger scatter of 0.255 dex. We also see correlations between the best-fitting or true mass within 150 pc and the true inner slopes, but with large scatters. The scatter can be improved if tracers within 150 pc are available.

We acknowledge the science research grants from the China Manned Space (CSST) Project with No. CMS-CSST-2021-B03. This work is also supported in part by the CSST Project grant No. CMS-CSST-2021-A02, NSFC (12022307, 12273021), the National Key Basic Research and Development Program of China (No. 2018YFA0404504), 111 project (No. B20019) and Shanghai Natural Science Foundation (No. 19ZR1466800). We thank the sponsorship from Yangyang Development Fund. W.W. is grateful for useful discussions with Xiaoting Fu, Leo Girardi, Zhiyuan Li, Lu Li, Fabo Feng, and Zhen Yuan. The computation of this work is carried out on the Gravity supercomputer at the Department of Astronomy, Shanghai Jiao Tong University, and is partly supported by the STFC DiRAC HPC Facility, at the Institute of Computational Cosmology (ICC), Durham University. L.Z. acknowledges funding from the National Key R&D Program of China under grant No. 2018YFA0404501, and the National Natural Science Foundation of China under grant No. Y945271001. R.G. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through the Spanish State Research Agency, under the Severo Ochoa Program 2020-2023 (CEX2019-000920-S). F.A.G. acknowledges support from ANID FONDECYT Regular 1211370 and by the ANID BASAL project FB210003. F.A.G. acknowledges funding from the Max Planck Society through a Partner Group grant.

Is the Core-cusp Problem a Matter of Perspective? Jeans Anisotropic Modeling against Numerical Simulations

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Data

2.1. Sample of Dwarf Systems, Mock Stars, and Mock Images

2.1.1. auriga Suite of Simulations

2.1.2. Dwarf Systems

2.1.3. Mock Stars

2.1.4. Mock Dwarf Images and Multi-Gaussian Expansion

2.2. Incorporating Observational Errors

2.2.1. Assigning Apparent Magnitudes of Individual Stars

2.2.2. Parallax, Line-of-sight Velocity, and Proper Motion Errors

2.2.3. Modeling of Foreground/background and Unbound Star Contamination

3. Methodology

4. Results

4.1. Overall Performance

4.2. Examples of the Best-fitting Matter Density Profiles

4.2.1. Best-fitting Density Profiles and Velocity Maps for Sag-like Systems

4.2.2. Best-fitting Density Profiles and Velocity Maps for Star-forming Systems with Outflows

4.3. Constraints Based on Fewer Tracers and Only Line-of-sight Velocities with Observational Errors

4.4. Constraints on Sag-like Systems after Including Observational Errors and Unbound Star Contamination

5. Discussion

5.1. Connection to the Core-cusp Problem

5.2. Detectability of Contractions and Outflows

6. Conclusions

Footnotes

Is the Core-cusp Problem a Matter of Perspective? Jeans Anisotropic Modeling against Numerical Simulations

Article metrics

Share this article

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Data

2.1. Sample of Dwarf Systems, Mock Stars, and Mock Images

2.1.1. auriga Suite of Simulations

2.1.2. Dwarf Systems

2.1.3. Mock Stars

2.1.4. Mock Dwarf Images and Multi-Gaussian Expansion

2.2. Incorporating Observational Errors

2.2.1. Assigning Apparent Magnitudes of Individual Stars

2.2.2. Parallax, Line-of-sight Velocity, and Proper Motion Errors

2.2.3. Modeling of Foreground/background and Unbound Star Contamination

3. Methodology

4. Results

4.1. Overall Performance

4.2. Examples of the Best-fitting Matter Density Profiles

4.2.1. Best-fitting Density Profiles and Velocity Maps for Sag-like Systems

4.2.2. Best-fitting Density Profiles and Velocity Maps for Star-forming Systems with Outflows

4.3. Constraints Based on Fewer Tracers and Only Line-of-sight Velocities with Observational Errors

4.4. Constraints on Sag-like Systems after Including Observational Errors and Unbound Star Contamination

5. Discussion

5.1. Connection to the Core-cusp Problem

5.2. Detectability of Contractions and Outflows

6. Conclusions

Footnotes