HSTPROMO Internal Proper Motion Kinematics of Dwarf Spheroidal Galaxies: II. Velocity Anisotropy and Dark Matter Cusp Slope of Sculptor

TL;DR

This study uses 20 years of HST data to measure 3D velocities in Sculptor, constraining its dark matter profile and indicating a core rather than a cusp, with results sensitive to galaxy inclination.

Contribution

First radially-resolved 3D velocity dispersion profiles for Sculptor, revealing dark matter density slope and inclination effects using combined proper motion and line-of-sight data.

Findings

Dark matter density slope $ ightarrow$ core-like profile

Inclination degeneracy affects dark matter profile constraints

Higher-order velocity moments help constrain galaxy inclination

Abstract

We analyze three epochs of HST imaging over 20 years for the Sculptor dwarf spheroidal galaxy, measuring precise proper motions for 119 stars and combining them with 1760 existing line-of-sight velocities. This catalog yields the first radially-resolved 3D velocity dispersion profiles for Sculptor. We confirm mild oblate rotation, with major-axis velocities reaching $\sim 2$ km s$^{-1}$ beyond 20.0 arcmin. Using a methodology similar to that in the first paper in this series, we solve the Jeans equations in oblate axisymmetric geometry to infer the galaxy's mass profile. Our modeling reveals a significant degeneracy due to the unknown galaxy inclination, which is overlooked under spherical symmetry assumptions. This degeneracy allows acceptable fits across a range of dark matter profiles, from cuspy to cored. While we do not directly constrain the inclination with our Jeans models, higher-order line-of-sight velocity moments provide useful additional constraints: comparisons with scalefree models from de Bruijne et al. (1996) favor highly flattened (more face-on) configurations. Adopting an inclination well consistent with these comparisons ($i = 57.1$ degrees), we find, alongside radial velocity anisotropy, a dark matter density slope of $\Gamma_{\rm dark} = 0.29^{+0.31}_{-0.41}$ within the radial extent of the 3D velocity data, ruling out a cusp with $\Gamma_{\rm dark} \leq -1$ at 99.8% confidence. This confidence increases for lower inclinations and decreases drastically for nearly edge-on configurations. The results qualitatively agree with $\Lambda$CDM, SIDM, and Fuzzy DM scenarios that predict core formation, while our specific measurements provide quantitative constraints on the prescriptions of feedback, cross sections, or particle masses required by these models, respectively.

Tables

Table 1. Table 1: Overview of Sculptor structural parameters.

Reference	Data	${\alpha }_0$	${\delta }_0$	$\theta$	$ϵ$	$R_h$
(1)	(2)	(3)	(4)	(5)	(6)	(7)
This work	Gaia EDR3	$15\stackrel{\circ }{.}02328\pm 8\stackrel{\prime \prime }{.}2$	$-33\stackrel{\circ }{.}71504\pm 4\stackrel{\prime \prime }{.}8$	$95\stackrel{\circ }{.}4\pm 1\stackrel{\circ }{.}6$	$0.272\pm 0.013$	$11\stackrel{\prime }{.}12\pm 0\stackrel{\prime }{.}17$
Muñoz et al. (2018)	MegaCam	$15\stackrel{\circ }{.}01830\pm 0\stackrel{\prime \prime }{.}3$	$-33\stackrel{\circ }{.}71860\pm 2\stackrel{\prime \prime }{.}6$	$92\stackrel{\circ }{.}0\pm 1\stackrel{\circ }{.}0$	$0.330\pm 0.010$	$11\stackrel{\prime }{.}17\pm 0\stackrel{\prime }{.}05$

Table 2. Table 2: Main results of the JamPy axisymmetric Jeans modeling.

$i$	$D$	$\Omega$	${\beta }_J$	$r_{\mathrm{dark}}$	$M_{\mathrm{dark}}^{R_{\max }}$	${\gamma }_{\mathrm{dark}}$	${\Gamma }_{\mathrm{dark}}$	$\overline{{\beta }_B}$	$v_{\mathrm{circ}}^{R_{\max }}$
[∘]	[kpc]			$\left[{10}^2\mathrm{pc}\right]$	$\left[{10}^8{M}_{\odot }\right]$				$\left[\mathrm{km}{s}^{-1}\right]$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
$43.7$	${87.69}_{-1.96}^{+1.92}$	$-{2.20}_{-0.37}^{+0.32}$	${0.98}_{-0.00}^{+0.00}$	${3.32}_{-0.32}^{+0.55}$	${1.68}_{-0.14}^{+0.16}$	${1.35}_{-0.62}^{+0.45}$	$-{0.16}_{-0.20}^{+0.15}$	${\mathrm{\hspace{0.17em}0.94}}_{-0.00}^{+0.00}$	${20.79}_{-0.79}^{+0.87}$
$47.7$	${85.98}_{-1.94}^{+1.94}$	$-{0.53}_{-0.20}^{+0.19}$	${0.92}_{-0.01}^{+0.01}$	${4.26}_{-0.32}^{+0.42}$	${2.11}_{-0.21}^{+0.25}$	${1.64}_{-0.45}^{+0.23}$	${0.52}_{-0.22}^{+0.19}$	${\mathrm{\hspace{0.17em}0.82}}_{-0.03}^{+0.02}$	${23.49}_{-1.26}^{+1.13}$
$52.4$	${85.79}_{-1.88}^{+2.40}$	$-{0.42}_{-0.24}^{+0.26}$	${0.86}_{-0.03}^{+0.03}$	${4.49}_{-0.42}^{+0.69}$	${2.04}_{-0.27}^{+0.24}$	${1.58}_{-0.37}^{+0.28}$	${0.57}_{-0.20}^{+0.19}$	${\mathrm{\hspace{0.17em}0.68}}_{-0.07}^{+0.06}$	${22.80}_{-1.26}^{+1.49}$
$57.1$	${85.46}_{-1.81}^{+2.07}$	$-{0.91}_{-0.75}^{+0.45}$	${0.72}_{-0.12}^{+0.07}$	${4.06}_{-0.51}^{+0.65}$	${1.50}_{-0.30}^{+0.32}$	${1.45}_{-0.68}^{+0.41}$	${0.29}_{-0.41}^{+0.31}$	${\mathrm{\hspace{0.17em}0.35}}_{-0.39}^{+0.17}$	${19.88}_{-1.93}^{+1.97}$
$61.8$	${85.78}_{-1.92}^{+1.79}$	$-{4.76}_{-4.75}^{+3.08}$	$-{0.10}_{-1.06}^{+0.60}$	${6.37}_{-2.64}^{+10.14}$	${0.95}_{-0.17}^{+0.25}$	$-{1.12}_{-0.75}^{+2.09}$	$-{1.39}_{-0.51}^{+1.18}$	$-{1.04}_{-2.15}^{+0.99}$	${16.05}_{-1.42}^{+1.89}$
$66.5$	${85.73}_{-1.73}^{+1.75}$	$-{6.35}_{-2.89}^{+2.64}$	$-{0.79}_{-0.67}^{+0.71}$	${12.17}_{-5.82}^{+60.81}$	${0.90}_{-0.18}^{+0.23}$	$-{1.78}_{-0.15}^{+0.65}$	$-{1.83}_{-0.13}^{+0.44}$	$-{1.91}_{-1.39}^{+1.36}$	${15.67}_{-1.52}^{+1.77}$
$71.2$	${85.69}_{-2.01}^{+2.02}$	$-{5.42}_{-2.65}^{+2.83}$	$-{0.60}_{-0.65}^{+0.72}$	${9.64}_{-4.84}^{+27.90}$	${0.88}_{-0.16}^{+0.18}$	$-{1.65}_{-0.28}^{+1.10}$	$-{1.74}_{-0.21}^{+0.69}$	$-{1.92}_{-1.10}^{+1.34}$	${15.47}_{-1.37}^{+1.46}$
$75.9$	${85.28}_{-1.86}^{+1.87}$	$-{5.48}_{-2.40}^{+1.91}$	$-{0.77}_{-0.57}^{+0.55}$	${13.05}_{-6.06}^{+51.39}$	${0.88}_{-0.18}^{+0.18}$	$-{1.79}_{-0.16}^{+0.50}$	$-{1.85}_{-0.12}^{+0.34}$	$-{2.22}_{-0.93}^{+1.04}$	${15.56}_{-1.57}^{+1.35}$
$80.6$	${85.48}_{-1.83}^{+1.82}$	$-{5.71}_{-5.48}^{+1.83}$	$-{0.83}_{-0.50}^{+0.54}$	${15.02}_{-7.24}^{+104.45}$	${0.88}_{-0.18}^{+0.18}$	$-{1.84}_{-0.11}^{+0.42}$	$-{1.87}_{-0.09}^{+0.29}$	$-{2.18}_{-0.88}^{+1.00}$	${15.48}_{-1.53}^{+1.45}$
$85.3$	${85.45}_{-1.76}^{+1.84}$	$-{5.83}_{-11.00}^{+1.81}$	$-{0.88}_{-0.47}^{+0.48}$	${16.18}_{-7.63}^{+93.35}$	${0.89}_{-0.19}^{+0.16}$	$-{1.87}_{-0.09}^{+0.33}$	$-{1.90}_{-0.08}^{+0.24}$	$-{2.25}_{-0.91}^{+1.14}$	${15.62}_{-1.62}^{+1.21}$
$90.0$	${85.53}_{-1.81}^{+1.90}$	$-{5.47}_{-4.00}^{+1.79}$	$-{0.86}_{-0.51}^{+0.50}$	${16.64}_{-8.33}^{+173.46}$	${0.90}_{-0.19}^{+0.18}$	$-{1.86}_{-0.10}^{+0.36}$	$-{1.88}_{-0.09}^{+0.25}$	$-{2.16}_{-0.98}^{+1.23}$	${15.65}_{-1.66}^{+1.40}$
	${85.83}_{-1.91}^{+2.02}$	$-{3.16}_{-4.96}^{+1.75}$	${0.13}_{-1.15}^{+0.78}$	${7.86}_{-7.86}^{+8.51}$	${1.39}_{-0.58}^{+0.68}$	$-{0.02}_{-1.82}^{+1.76}$	$-{0.67}_{-1.24}^{+1.23}$	$-{0.56}_{-2.00}^{+1.33}$	${18.88}_{-3.94}^{+4.23}$

Table 3. Table 3: Proper motion catalog.

$\alpha$	$\delta$	${\mu }_{\alpha ,\ast }$	${ϵ}_{{\mu }_{\alpha ,\ast }}$	${\mu }_{\delta }$	${ϵ}_{{\mu }_{\delta }}$	F775W	F606W	ID
[∘]	[∘]	[$\mathrm{mas}{\mathrm{yr}}^{-1}$]	[$\mathrm{mas}{\mathrm{yr}}^{-1}$]	[$\mathrm{mas}{\mathrm{yr}}^{-1}$]	[$\mathrm{mas}{\mathrm{yr}}^{-1}$]	[dex]	[dex]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
$14.98182$	$-33.84174$	$0.0259$	$0.0189$	$-0.0052$	$0.0196$	$-11.3723$	$-11.4352$	1
$14.99856$	$-33.83554$	$0.0238$	$0.0225$	$-0.0191$	$0.0227$	$-11.5420$	$-11.0647$	2
$15.00169$	$-33.84224$	$0.0128$	$0.0189$	$0.0222$	$0.0155$	$-11.3253$	$-10.8545$	3
⋮

Table 4. Table 4: Line-of-sight velocity catalog.

$\alpha$	$\delta$	$v_{\mathrm{LOS}}$	${ϵ}_{v_{\mathrm{LOS}}}$	$N_{\mathrm{spec}}$	IDTolstoy	IDWalker
[∘]	[∘]	[$\mathrm{km}{s}^{-1}$]	[$\mathrm{km}{s}^{-1}$]
(1)	(2)	(3)	(4)	(5)	(6)	(7)
$13.66484$	$-33.52300$	$120.5981$	$0.4182$	$2$	scl011_09_	–
$13.87452$	$-33.76831$	$119.8343$	$0.9640$	$3$	scl_11_040	–
$13.93013$	$-34.02570$	$123.5945$	$0.4500$	$1$	scl020-12_	–
⋮

Equations

$$x_i\equiv \Delta {\mu }_i/{ϵ}_{\Delta {\mu }_i},$$

$$\mathrm{MAD}=\frac{{\eta }_{50}(|x-{\eta }_{50}(x)|)}{0.6745}$$

$$\frac{\partial \left(\nu ⟨v_r^2⟩\right)}{\nu \partial r}+\frac{(1+{\beta }_{\mathrm{J}})⟨v_r^2⟩-⟨v_{\varphi }^2⟩}{r}=-\frac{\partial \Phi }{\partial r},$$

$$\frac{(1-{\beta }_{\mathrm{J}})}{\nu }\left[\frac{\partial \left(\nu ⟨v_r^2⟩\right)}{\partial \theta }+\frac{\nu ⟨v_r^2⟩}{\tan \theta }\right]-\frac{⟨v_{\varphi }^2⟩}{\tan \theta }=-\frac{\partial \Phi }{\partial \theta },$$

$${\beta }_{\mathrm{J}}\equiv 1-\frac{⟨v_{\theta }^2⟩}{⟨v_r^2⟩}=1-\frac{{\sigma }_{\theta }^2}{{\sigma }_r^2},$$

$${\beta }_{\mathrm{B}}\equiv 1-\frac{⟨v_{\theta }^2⟩+⟨v_{\varphi }^2⟩}{2⟨v_r^2⟩},$$

$${⟨v_{\varphi }⟩}^2$$

$${\sigma }_{\varphi }^2$$

$$\mathcal{L}$$

$$\times \prod _{i,\mathrm{POS}}{G}_{\mathrm{POSr}}({\mu }_{\mathrm{POSr},i})G_{\mathrm{POSt}}({\mu }_{\mathrm{POSt},i}).$$

$${\Gamma }_{\mathrm{dark}}\equiv \frac{{\int }_{r_{\min }}^{r_{\max }}\frac{\mathrm{d}\log \rho }{\mathrm{d}\log r}\rho (r)\mathrm{d}r}{{\int }_{r_{\min }}^{r_{\max }}\rho (r)\mathrm{d}r},$$

$$v_{\mathrm{circ}}(r)=\sqrt{\frac{GM(r)}{r}},$$

Full text

dark matter — galaxies: dwarf — galaxies: structure — methods: data analysis — proper motions — stars: kinematics and dynamics

††facilities: HST; Gaia; Magellan:Clay; VLT:Kueyen. ††software: Python (Van Rossum & Drake, 2009), BALRoGO (Vitral, 2021), JamPy (Cappellari, 2020), emcee (Foreman-Mackey et al., 2013), Scipy (Jones et al., 2001–), Numpy (van der Walt et al., 2011), Matplotlib (Hunter, 2007), Agama (Vasiliev, 2019), aplpy (Robitaille & Bressert, 2012), scalefree (de Bruijne et al., 1996), and the Scientific color maps from Crameri, Shephard, & Heron (2020). ††thanks: Royal Society Newton International Fellow††thanks: NASA Hubble Fellow

I Introduction

The nature of dark matter (DM) remains one of the most profound open questions in modern astrophysics, with far-reaching implications for both cosmology and our understanding of galaxy formation (Bertone, Hooper, & Silk, 2005; Silk & Mamon, 2012). Among the environments in which DM can be studied, dwarf spheroidal galaxies (dSphs) stand out as particularly valuable laboratories (Strigari, 2013; Collins & Read, 2022). These small, faint systems are among the most DM-dominated objects known (Pryor & Kormendy, 1990), making them highly sensitive to the properties of DM on small scales (Bullock & Boylan-Kolchin, 2017). Their relatively simple stellar populations and quiescent evolutionary histories (Tolstoy, Hill, & Tosi 2009; Savino et al. 2025) also provide a cleaner view of internal dynamical processes than the more complex environments found in larger galaxies (Grebel, 2009).

Accurately inferring the internal kinematics of dSphs remains a long-standing challenge. Analyses based solely on line-of-sight (LOS) velocities111Motions toward or away from the observer, measured through spectroscopic data typically obtained from ground-based observatories. are fundamentally limited by the well-known degeneracy between the mass distribution and the orbital anisotropy of the stellar component (Binney & Mamon, 1982). This degeneracy hampers efforts to robustly distinguish between different DM density profiles, complicating tests of predictions from $\Lambda$ Cold Dark Matter ($\Lambda$CDM) and alternative cosmological models (see Salucci 2019 and Boldrini 2021 for reviews and references). Overcoming this limitation typically requires either modeling techniques that incorporate higher-order velocity moments (Chanamé, Kleyna, & van der Marel, 2008; Read & Steger, 2017) – which remain difficult to constrain with current LOS data – or the inclusion of transverse stellar motions in the plane of the sky (van der Marel & Anderson, 2010). The latter has historically been inaccessible for galaxies beyond the Milky Way (Gaia Collaboration et al., 2018), owing to the simultaneous need for long time baselines (typically $\gtrsim 10$ yrs) and high-precision astrometry to reach sufficiently small uncertainties.

The advent of space-based astrometry has now opened new avenues for addressing this limitation. Transverse velocity measurements, or proper motions (PMs), can be obtained with high precision and, when combined with LOS velocities, enable the reconstruction of full three-dimensional velocity vectors for stars in nearby dSphs (Massari et al., 2017, 2020; del Pino et al., 2022; Libralato et al., 2023; Vitral et al., 2024; McKinnon et al., 2024). These additional kinematic dimensions significantly enhance the constraining power of dynamical models and allow for a more complete and direct view of a system’s internal dynamics than can be obtained from LOS velocities alone (Read et al., 2021).

In a series of papers beginning with Vitral et al. (2024, hereafter Paper I), we aim to construct comprehensive PM catalogs for nearby dSphs and use them, alongside LOS velocities, to build detailed dynamical models of these systems. By systematically applying these techniques, we seek to refine our understanding of their underlying DM distributions and evaluate the influence of baryonic processes in shaping both their internal kinematics and DM density profiles. In the first paper of this series, we analyzed the Draco dSph and found that its internal dynamics are consistent with a cuspy DM profile, as predicted by DM-only $\Lambda$CDM simulations (Navarro, Frenk, & White, 1997). In this second study, we turn our attention to the Sculptor dSph (also known as PGC3589, MCG-06-03-015, or ESO351-30), a well-studied system (Battaglia et al., 2008; Łokas, 2009; Walker & Peñarrubia, 2011; Richardson & Fairbairn, 2014; Massari et al., 2017; Arroyo-Polonio et al., 2025) which, until now, lacked a large and homogeneous PM dataset with uncertainties smaller than its intrinsic velocity dispersion. Such precision is essential to mitigate the impact of measurement errors and systematic biases in dynamical mass estimates (Watkins et al., 2013; Vitral et al., 2022), as discussed in Paper I. Among Sculptor’s notable features is its higher stellar mass relative to Draco (Martin et al., 2008; de Boer et al., 2012), placing the two galaxies on opposite sides of the threshold where baryonic processes are thought to flatten DM cusps (Fitts et al., 2017), despite their broadly similar structural properties (Muñoz et al., 2018).

With this background in mind, the remainder of the paper is structured as follows. Section II outlines the general properties of Sculptor relevant to our modeling and describes the kinematic datasets used throughout the study. Section III details the modeling techniques employed to analyze the galaxy. Sections IV and V present our main results and assess their robustness, respectively. Section VI discusses the implications of our findings for mass modeling and broader cosmological questions. Finally, Section VII summarizes our key conclusions.

II Sculptor Data and General Characteristics

II.1 Structural Parameters

We begin by defining and presenting the structural parameters used as input throughout this work – namely, the shape of the stellar mass density profile, characterized by the total stellar mass, its scale radius, and the position of the galaxy’s center of mass (hereafter simply center). The position of Sculptor’s center, as listed in the McConnachie (2012) catalog and the Local Volume Database (Pace, 2024), originates from its initial discovery by Shapley (1938). Recently, Muñoz et al. (2018) refined Sculptor’s center and structural parameters using Canada France Hawaii Telescope’s MegaCam imaging. Their analysis provided more reliable estimates based on improved stellar measurements. In this study, we adopt a similar fitting methodology to Muñoz et al. (2018), also employed in Paper I (see the latter’s section 2.1.1 and eq. A2). Specifically, we fit an axisymmetric Plummer (1911) model to Gaia EDR3 star counts, incorporating a constant-density interloper background to derive Sculptor’s structural parameters and central position.

Our goodness-of-fit results are summarized in Figure 1, which compares the Gaia EDR3 data with our best-fitting structural model. The top panel shows a centered sky view of Sculptor, where individual Gaia sources are color-coded by the local stellar count within $0\farcm 5$. The fitted center (white cross) lies close to the observed density maximum, and the model’s 2D half-number ellipse (white dashed line) traces the projected stellar distribution well. The lower panel presents the azimuthally averaged surface number density as a function of projected radius $R$. The observed profile (blue) is well matched by the sum (solid red) of two components: the dSph itself (dashed red) and a uniform interloper background (dotted red). The good agreement in both the spatial distribution and radial profile indicates that our model captures Sculptor’s sky orientation, central position, and overall surface density structure. The parameters derived from this fit are presented in Table 1, alongside those reported by Muñoz et al. (2018).

The overall shape of Sculptor – quantified here by its position angle $\theta$, ellipticity $\epsilon$, and projected half-number radius $R_{h}$ – agrees closely with the results of Muñoz et al. (2018). In our convention, $\epsilon\equiv 1-b/a$, where $a$ and $b$ are the major and minor axes of the best-fitting photometric ellipse, and $\theta$ is the position angle of the major axis measured from north toward east. The half-number radius $R_{h}$ corresponds to the major axis length of the ellipse enclosing half the stars in a Plummer model fit. Our fitted center lies in between the Shapley (1938) and Muñoz et al. (2018) values, with respective distances of $52\farcs 5$ and $19\farcs 7$ from each.

The flattening of Sculptor closely resembles that of Draco, resulting in a probability distribution function (PDF) of possible inclinations that follows a similarly shaped curve. This PDF is calculated using the formalism from Paper I (see section 2.1.3), which is based on the statistical distribution of elliptical galaxy shapes inferred by Lambas et al. (1992). The resulting distribution is shown in Figure 2. The $[16,\,50,\,84]$ percentiles for the inclination and 3D flattening distributions are $[48\fdg 34,\,59\fdg 50,\,77\fdg 82]$ and $[0.38,\,0.60,\,0.71]$, respectively.

Finally, we consider Sculptor’s stellar mass as estimated by de Boer et al. (2012) based on its star formation history and color-magnitude diagram (CMD). They report a stellar mass of $7.8\times 10^{6}~\rm M_{\odot}$ within a $1\fdg 0$ elliptical radius.222Their definition of an elliptical radius, clarified through private communication with T. J. L. de Boer, involves projecting the ($\alpha$, $\delta$) data onto ($x$, $y$) Cartesian coordinates, rotating them by Sculptor’s position angle, and dividing the $y$ coordinate by $(1-\epsilon)$. The elliptical radius is then calculated as $\sqrt{x^{\prime 2}+y^{\prime 2}}$, where $x^{\prime}$ and $y^{\prime}$ are the rotated and scaled Cartesian coordinates. Extrapolating this measurement to infinity, assuming an axisymmetric Plummer density profile, yields a total stellar mass of $8.08\times 10^{6}~\rm M_{\odot}$, which is adopted throughout this work.333For context, this value relates to a DM halo mass of $M_{\rm dark}=4.85\times 10^{9}~\rm M_{\odot}$ in $\Lambda$CDM cosmology, using the relation for low-mass dwarfs from Read et al. (2017).

II.2 Line-of-Sight Velocities

Similarly to Paper I, we use the spectroscopic catalog from Walker et al. (2023) to obtain LOS velocities for our mass modeling. Following the same data cleaning procedures detailed in Paper I (see section 2.2.1), we identify 1151 Sculptor members from this catalog. To augment this dataset, we incorporate the recently published spectroscopic measurements by Tolstoy et al. (2023), which are derived from Gaia and VLT/FLAMES observations. We clean the Tolstoy et al. (2023) catalog following their recommended criteria: setting the $z$ coefficient444Tolstoy et al. (2023) define $z$ as a membership score based on Gaia’s PMs and parallaxes. $\leq 14.2$, enforcing ${\rm S/N}>20$, and selecting only stars classified as members in their analysis. This procedure yields an additional 1123 LOS velocities. For both catalogs, we apply corrections for perspective motions using the formalism from van der Marel et al. (2002, eq. 13).

To leverage the strengths and coverage of both datasets, we merge them using the following approach:

We first identify common stars in the two catalogs by performing a symmetric sky match with a maximum separation of $1\farcs 0$. This matched subset serves as the basis for comparison. For each catalog, we compute zero-points for LOS velocities using the following steps:

1. (a)

Compute the weighted mean of the LOS velocities for the matched stars in each catalog. 2. (b)

For each catalog, calculate the distribution of velocity offsets as the difference between the LOS velocities of individual stars and the weighted mean of the matched subset. 3. (c)

Derive the zero-point correction for each catalog by taking the weighted mean of these velocity offset distributions. 2. 2.

The computed zero-point corrections are applied to the LOS velocities in each catalog,555This common treatment is done since we lack an independent means of determining which catalog’s zero-point is intrinsically more accurate. Instead, we align both to a common scale using uncertainty-weighted averages, which incorporate each catalog’s reported precision. The resulting shifts are negligible compared to the velocity scales of interest in this work (e.g., $\sim 2~\rm km\,s^{-1}$ rotation and $\sim 10~\rm km\,s^{-1}$ dispersion): $0.03\pm 0.01~\rm km\,s^{-1}$ for the Walker et al. (2023) data and $-0.32\pm 0.03~\rm km\,s^{-1}$ for the Tolstoy et al. (2023) data. ensuring consistency between the datasets. 3. 3.

Finally, we construct the merged catalog by treating matched stars as single entries. For these matched stars, already corrected for the zero-point, we calculate a common LOS velocity and its associated uncertainty using the weighted mean of their individual measurements.666We do not simply adopt the entry from the catalog with the smallest nominal error, as uncertainty-weighted averaging makes optimal use of both measurements and yields the most precise estimate. Stars without matches are retained in the catalog unchanged, aside from the applied zero-point adjustment.

The resulting combined catalog contains 1,760 Sculptor members with LOS velocity uncertainties smaller than the galaxy’s intrinsic velocity dispersion.777This threshold is critical for robust mass modeling, as highlighted in Paper I and other studies (e.g., Vitral et al. 2022, 2023a). This merged dataset represents the most comprehensive spectroscopic catalog assembled for any dwarf galaxy to date.

II.2.1 Binaries

In recent years, unresolved binaries in stellar systems have garnered increased attention due to their potential role in inflating velocity dispersion profiles, which can result in overestimating the mass budget (e.g., Rastello et al., 2020; Pianta et al., 2022). However, our previous analyses of Draco showed no significant impact from binaries, despite its high binary fraction (Spencer et al., 2018). This raises the question of how such inflation patterns might depend on various data characteristics, such as the velocity dispersion and the number of measured epochs (Minor et al., 2010; Wang et al., 2023; Arroyo-Polonio et al., 2023).

Following the approach in Paper I, we assess the influence of unresolved binaries on Sculptor’s dynamics by comparing LOS velocity dispersion profiles for stars grouped by their number of measured epochs. As discussed in our earlier work, if binaries contribute to inflated velocity dispersions, we would expect stars with single-epoch LOS velocity measurements to exhibit higher dispersions than those with multi-epoch measurements, where the effects of binary motion are largely averaged out (Wang et al., 2023).

Figure 3 shows that, despite Sculptor’s high binary fraction,888Arroyo-Polonio et al. (2023) recently estimated a binary fraction of $0.55^{+0.17}_{-0.19}$ for this dwarf galaxy. the velocity dispersion profiles derived from stars with different numbers of epochs remain consistent within 1-$\sigma$ uncertainties. This finding aligns with the previous results for Draco. We therefore conclude that Sculptor’s binary population does not introduce any bias into our mass-anisotropy estimates beyond the level of statistical uncertainties.

II.2.2 Rotation

Beyond the perspective rotation due to Sculptor’s motion in the sky, which we have removed in Section II.2, this dSph might still have residual internal rotation that can be observed through the variation of its LOS first-order velocity moments. Indeed, recent studies have suggested that Sculptor carries mild rotation in the LOS (e.g., $\left<v/\sigma\right>=0.15\pm 0.15$, Zhu et al., 2016), with a spin axis aligned in between typical oblate/prolate configurations (Arroyo-Polonio et al., 2024, figure 6).

While the Jeans modeling described in Section III.1 directly accounts for the LOS first-order velocity moments, we do opt to model Sculptor as an oblate spheroid.999The Jeans modeling code we use, JamPy, is currently designed to handle either oblate or prolate configurations, but not intermediate arrangements. To ensure this assumption is reasonable, we first analyze Sculptor’s rotation profile. This analysis includes determining the direction of the rotation axis and the LOS velocity amplitude along the galaxy’s photometric major axis, allowing for further comparison with our Jeans fits.

We model rotation following a similar approach to that in Paper I, by defining eight concentric annuli with linearly-spaced projected radii. For each annulus, we fit a sinusoidal function to the LOS velocity data. While each annulus is assigned an independent amplitude, we enforce a common offset and phase across all annuli to recover the bulk LOS velocity and the global spin axis, respectively. In our formalism, we first project the galaxy onto Cartesian coordinates on the sky (as in eq. 2 from Gaia Collaboration et al. 2018), and then align it with Sculptor’s major and minor axes, measuring angles clockwise from the positive $x$-axis (associated with an angle of $0^{\circ}$) to the positive $y$-axis. 101010For Sculptor, these directions approximately align with the East and North directions, respectively. To better estimate the uncertainties in our fits, we employ a bootstrap methodology with 10,000 realizations.

Figure 4 shows the fitted rotation axis direction alongside a smoothed map of LOS velocities for the stars included in our modeling. We determine a rotation axis direction of $\phi_{0}=302\fdg 0\pm 52\fdg 0$, consistent with an oblate configuration (i.e., $\phi_{0}=270\fdg 0$). This result is in good agreement with recent estimates by Arroyo-Polonio et al. (2024), whose spin axis direction, expressed within our formalism, corresponds to $331\fdg 4\pm 18\fdg 0$. The rotation amplitude along the photometric major axis is presented later in Figure 8 together with our PM data, but a brief look reveals a modest but non-negligible rotation signal of $\sim 2~\rm km\,s^{-1}$, which amounts to approximately 20% of Sculptor’s internal velocity dispersion. This low rotation amplitude also explains why the rotation axis direction remains poorly constrained in our fits.

II.3 Proper Motions

II.3.1 Observations and Astrometric Catalogs

For our new PM measurements of Sculptor stars, we utilized multi-epoch HST ACS/WFC imaging data. Details about the field locations and earlier epoch observations for our target fields, F1 and F2, are provided in Sohn et al. (2017). The locations of these fields are also shown in Figure 5. In brief, both fields had two epochs of imaging data obtained in 2002 and 2013. Both 2002 and 2013 observations were performed using the F775W filter, and for 2013 observations, an additional filter F606W was used to create CMDs. These fields were revisited in October 2022 and October 2023 through our HST program GO-16737 (Sohn et al., 2021), using the same telescope pointing and orientation as in the previous epochs. During this latest epoch, we acquired 12 individual exposures using the F775W filter for each field, with each exposure lasting 495 seconds.

Data analysis largely followed the procedures outlined in Bellini et al. (2018), Libralato et al. (2018), and Paper I. Here, we provide a high-level summary of the PM derivation process and refer the reader to those works for methodological details. We downloaded the flat-fielded _flt.fits images for all target fields and epochs from the Mikulski Archive for Space Telescopes (MAST). These were processed using the hst1pass program (Anderson, 2022) to derive positions and fluxes for each star in every exposure. Instead of using the _flc.fits images, which are corrected for charge transfer efficiency (CTE) losses, we employed the table-based CTE correction option in hst1pass. This is an improved version of the corrections used in earlier works (Anderson 2022; Anderson, in prep.). Positional corrections were applied using the ACS/WFC geometric distortion solutions from Kozhurina-Platais et al. (2015), further extended to account for time-dependent distortion variations beyond 2020 (V. Kozhurina-Platais, private communication).

For each field, a “master frame” was constructed using the average positions of stars from the repeated first-epoch exposures. The ($X,Y$) axes of these master frames were aligned to ($\alpha,\delta$) coordinates by registering stellar positions to the Gaia DR3 astrometric system. Positions of stars from subsequent epochs were aligned to these master frames using a six-parameter linear transformation, and average positions were computed for each epoch. This procedure inherently aligns the star fields across epochs, resulting in zero PM on average for Sculptor dSph stars. Since our focus is on measuring internal velocity dispersion in the plane of the sky, this alignment does not affect our results. Uncertainties on average positions were determined as the root mean square of repeated measurements divided by the square root of the number of exposures. Finally, for each field and epoch, we prepared a catalog containing star positions as derived above, along with average instrumental magnitudes111111The instrumental magnitude in a given filter is defined as ${\rm mag}=-2.5\log c$, where $c$ is the electron count per exposure for a source. in F606W and F775W (from previous epochs), as output by hst1pass.

II.3.2 Photometric Cleaning & Local Corrections

With the catalogs catalogs of star positions, we compute raw PMs through a least-squares line fit of the master frame ($X,Y$) positions as a function of the epoch time. We use the polyfit routine from NumPy (van der Walt et al., 2011), assuming measured ($X,Y$) uncertainties, and no $\chi^{2}$ re-scaling. Next, the procedures for cleaning the photometry, removing CMD outliers, and applying local PM corrections to account for HST systematics followed the methodology outlined in Paper I (sections 2.3.2 and 2.3.3). Below, we briefly summarize the key steps, which are applied independently to each field.

We begin by identifying and removing clear outliers in the CMD and in the magnitude-QFIT121212QFIT quantifies a combined measure of goodness-of-fit and S/N; see Anderson et al. (2006); Libralato et al. (2014) for details. diagrams, employing a friends-of-friends algorithm. 2. 2.

Next, we address HST instrument systematics, which often appear as spatially-correlated mean motions of stars across different regions of the telescope detectors (e.g. Bellini et al., 2018; Libralato et al., 2018). To correct for this, we calculate and apply median-based shifts, with uncertainties propagated to both PM directions, using a local network of neighboring stars.131313For clarity, this correction is not applied iteratively but only once, during the final stage of the PM construction. As described in Paper I, the correction introduces an additional spread in the data, which is accounted for in the final PM uncertainty through a quadrature sum. To assess the reliability of this correction, we performed a dedicated Monte Carlo test, where we generated mock stellar fields with realistic uncertainties and injected spatial distortions mimicking detector systematics.141414The injected spatial distortions were constructed using a sum of Fourier modes across a mock stellar field. The procedure was repeated thousands of times with new realizations of velocity fields, star positions, and distortion phases to ensure statistical significance. We then applied our standard local correction procedure and compared the recovered velocity dispersions to the known input values. The test showed that the propagated uncertainties in our locally-corrected estimates are realistic, despite any unmodeled correlations that the corrections may introduce. Moreover, we find that the locally-corrected dispersions are unbiased to with a fraction of these propagated uncertainties. Hence, the final results of our work are insensitive to the details of our local corrections. 3. 3.

The effectiveness of these corrections is validated through diagnostic checks, such as reproducing the equivalent of figure 8 from Paper I. This process effectively removes first-order velocity moments, such as rotation or global streaming motions, while preserving the second-order moments necessary for accurate mass-anisotropy modeling. We also test the impact of this data manipulation later on Section V.2.1.

A key feature of our new data is that, while each field is initially processed independently, they include overlapping regions, as shown in Figure 5. This intersection allows for a direct comparison of matched stars. If the local corrections are effective and properly applied, the overall PM differences between matched stars should be consistent with their associated PM uncertainties and should not exhibit significant systematic shifts. Specifically, the distribution of PM differences, normalized by their combined error budget, is expected to approximate a Gaussian distribution with a mean of zero and a standard deviation of one.

To test this, we account for the non-Gaussian nature of the datasets, arising primarily from two factors: the absence of error cleaning in this initial subset of matched stars and the low-number statistics. These factors can influence the distribution of PM differences, potentially introducing asymmetries or heavier tails. To address these issues, we employ median-based statistics to constrain the center and spread of the data. Specifically, we define

[TABLE]

where the suffix $i$ stands for $\alpha,*$ or $\delta$151515Hereafter, $\alpha,*$ refers to the component $\mu_{\alpha,*}=({\rm d}\alpha/{\rm d}t)\cos{\delta}$ from the PM vector. and $\epsilon$ is used to refer to uncertainties, rather than the $\sigma$ symbol used for the standard deviation. Next, we use the median and the median absolute deviation (MAD) of the $x_{i}$ distribution to construct a corresponding Gaussian curve. The MAD, corrected to account for Gaussian standard deviations as defined in Beers, Flynn, & Gebhardt (1990), is

[TABLE]

where $\eta_{50}(x)$ refers to the 50th percentile (i.e., the median) of the distribution of $x$.

Figure 6 compares the Gaussian curve computed using the median and MAD (shown in red) with the expected normal distribution (shown in blue). The agreement between these two curves supports the validity of the local corrections applied to our PM dataset. Not only are the PM differences not significantly shifted toward positive or negative values, but the propagated uncertainties also adequately encompass the observed spread in the data. As a further check, we also performed this analysis using a skewed Gaussian fit,161616This test was motivated by the prominent tail towards ($x_{\alpha,*}$, $x_{\delta}$) $\sim(2.5,-2.5)$ visible in Figure 6. finding that the resulting mean and dispersion remained consistent with those of the expected normal distribution. Following this validation, we proceed to the next stage, which involves cleaning stars with high PM errors and removing outliers in PM.

II.3.3 Outliers and Underestimated Errors

Before constructing PM profiles, we perform a final round of data cleaning for our PM measurements, guided by the analyses presented in Paper I. Specifically, we compute the global velocity dispersion171717Here and throughout this work, the dispersion of a discrete dataset with uneven uncertainties is computed as described in appendix A of van der Marel & Anderson (2010), a procedure thoroughly validated in multiple works (e.g., Watkins et al. 2015; Vitral et al. 2023a and Paper I). of our PM dataset for each field, applying various thresholds for the maximum allowable PM uncertainty. We select the highest threshold that meets two criteria: (i) it remains within the same order of magnitude as the expected LOS velocity dispersion, for a given distance, and (ii) it lies in a plateau of velocity dispersions before increasing due to the inclusion of measurements with large formal errors, which – if underestimated – can bias the results. These criteria can be better visualized in Figure 7. As discussed in detail in Paper I and prior works (e.g. Watkins et al. 2015; Vitral et al. 2022), this step is critical for accurate mass modeling. Large formal errors, particularly when underestimated and exceeding the intrinsic velocity dispersion, can introduce significant biases into equilibrium-based mass measurements.

Based on this, we adopt a threshold of $\epsilon_{\mu}<0.023~\rm mas\,yr^{-1}$ on the PM error;181818This translates to a transverse velocity of $9.4~\rm km\,s^{-1}$, for a given distance of $85.8$ kpc (cf. Table 2). stars with larger errors are excluded from the subsequent analysis. Additionally, we apply a 3-$\sigma$ outlier cut, as in Paper I. We once again use median-based statistics (median and MAD) rather than the nominal mean and standard deviation to ensure robustness against outliers.191919Note that this step mainly serves to exclude stars too far from the galaxy’s bulk motion; in practice, only one of the stars removed lay within 5-$\sigma$ from the median. Finally, we treat matched stars from overlapping fields as single entries, following a similar approach to the merging of our LOS datasets (i.e. using the weighted mean of the individual measurements). Hence, after applying the rigorous selection criteria described throughout this section, we retained a high-quality subset of 119 stars from the original HST observations,202020For context, the original subset contained 8,317 and 5,087 stars in Fields 1 and 2, respectively, while the field-merged PM subset with no error cleaning whatsoever contained 2,070 stars. all of which having proper-motion uncertainties below the intrinsic velocity dispersion. This yields the most precise PM dataset ever assembled for this dwarf galaxy.

II.3.4 Proper Motion Dispersion Profiles

With our new PM dataset, we present for the first time in Figure 8 the radially-resolved transverse velocity dispersion profiles of the Sculptor dSph.212121For context, the error bars on the proper motion dispersion bins reflect comparable contributions from shot noise and measurement uncertainties. Although these profiles are less spatially extended than those shown for Draco in Paper I, they remain valuable for breaking degeneracies between the galaxy’s mass distribution and orbital structure – that is, the velocity anisotropy, which describes whether the stellar orbits are preferentially radial or tangential.

Quantitatively, we find a ratio of tangential to radial PM dispersions of $\left<\sigma_{\rm POSt}\right>/\left<\sigma_{\rm POSr}\right>=1.19\pm 0.19$ across the spatial extent of the data, where POSr and POSt refer to the plane-of-sky radial and tangential directions, respectively. The ratio of LOS to PM dispersions is $\left<\sigma_{\rm LOS}\right>/\left<\sigma_{\rm POS}\right>=0.93\pm 0.08$, with $\sigma_{\rm POS}$ denoting the average dispersion over both transverse directions. The former ratio is independent of the galaxy’s distance, while the latter scales inversely with it.222222The quoted value assumes $D=85.83$ kpc, as predicted by our Jeans fits in the final row of Table 2.

These measurements place important constraints on Sculptor’s internal velocity structure, particularly on its velocity anisotropy. As discussed in Sections III and IV, the relative values of the velocity dispersions in different directions – captured here by ratios such as $\left<\sigma_{\rm POSt}\right>/\left<\sigma_{\rm POSr}\right>$ and $\left<\sigma_{\rm LOS}\right>/\left<\sigma_{\rm POS}\right>$ – provide strong leverage for disentangling the orbital structure and intrinsic shape of the galaxy.

III Methods

III.1 Axisymmetric Jeans Modeling: JamPy

As in Paper I, we use the mass anisotropy modeling code JamPy (Cappellari, 2020), configured such that the velocity ellipsoid aligns with spherical coordinates. This choice follows our assumption of a spherical DM potential, which is expected to dominate the system’s dynamics when compared to Sculptor’s luminous axisymmetric component. In this configuration, JamPy applies the Jeans equations for rotating axisymmetric systems (e.g., Bacon, Simien, & Monnet 1983),

[TABLE]

where $\Phi$ is the gravitational potential, and $\left<.\right>$ denotes the distribution function-averaged quantity. The anisotropy parameter $\beta_{\rm J}$, assumed to be independent of $\theta$, is given by

[TABLE]

where the second equality assumes $\left<v_{\theta}\right>=\left<v_{r}\right>=0$. Due to symmetry and continuity, axisymmetric models always satisfy $\left<v_{\phi}\right>=0$ and $\left<v_{\phi}^{2}\right>=\left<v_{\theta}^{2}\right>$ along the symmetry axis. Consequently, along this axis, the anisotropy parameter $\beta_{\rm B}$, as defined by Binney (1980),

[TABLE]

equals $\beta_{\rm J}$. Models with $\beta_{\rm J}=0$ yield the same predicted second velocity moments as models in which the distribution function $f(E,Lz)$ does not depend on a third integral. Such models have been widely used for fitting data of axisymmetric systems (e.g. van der Marel, 1991). Away from the symmetry axis, models with $\beta_{\rm J}=0$ do not have an isotropic velocity dispersion tensor.

By not enforcing $\left<v_{\phi}\right>=0$ throughout the entire system, JamPy allows users to model rotation in the fitted galaxy – a feature often neglected in previous kinematic models of Sculptor (e.g., Massari et al., 2017; Hayashi et al., 2020). The first moment in the $\phi$ direction is related to the second-order moment in the radial direction through

[TABLE]

where $\Omega$ is the rotation parameter, also assumed to be independent of $\theta$. This parameter was denoted as $\gamma$ in Cappellari (2020), but we adopt a different notation here to avoid confusion with the inner slope of the DM mass density, which is commonly represented by the same symbol.

From the solutions to these equations, JamPy calculates projected velocity moments, which we then use to compute the respective quantities in the LOS and POS directions. We sample the parameter space using the Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al., 2013), by using a likelihood $\mathcal{L}$ defined as

[TABLE]

With respect to equation III.1, a few key considerations must be clarified:

• •

$G(x)$ represents a Gaussian function of argument $x$.

• •

Since the stars with PMs in our dataset do not have a corresponding $v_{\rm LOS}$ measurement,232323This is because our HST exposures target the more numerous faint stars, while the brighter stars for which ground-based spectroscopic measurements are typically possible are often saturated in the respective HST images. the product of the POS Gaussian components can be effectively treated separately from the LOS Gaussian component, with no risk of neglecting cross-terms in spherical or axisymmetric symmetry.

• •

For all Gaussian components, the respective standard deviation is derived from the output of JamPy along the corresponding dimension, using the relation $\sigma^{2}=\left<v^{2}\right>-\left<v\right>^{2}$. To this, we add the individual stellar errors quadratically, assuming Gaussian-distributed uncertainties.

• •

For the LOS component, the Gaussian mean is taken as the bulk LOS motion of Sculptor, plus the respective first-order moment estimated by JamPy, and adjusted to the pre-defined sense of rotation in Sculptor. For the POS components, however, we adopt a null mean as a result of our local cleaning process explained in Section II.3.2, which effectively removes streaming motions in our PMs.

• •

Although the assumption of Gaussian-distributed projected velocities is formally incorrect in the presence of anisotropy (Merritt, 1987, figure 8), we retain it due to the demonstrated convergence of this approach when applied to mock datasets generated by self-consistent distribution functions (Read et al., 2021), as well as its consistency with results from well-studied dynamical systems (see Watkins et al., 2013, for an application to $\omega$ Centauri).

Throughout our posterior sampling, we employ priors as defined in Section III.1.1. Our MCMC routine is configured to run a minimum of 20 000 steps per fit,242424If proper convergence is not achieved, we add an additional 10 000 steps, which was sufficient in our case. and we execute the runs in parallel on the local cluster at the University of Edinburgh. In these configurations, each run takes approximately 3 days to complete. We repeat the process for a different set of possible inclinations, drawn from the PDF derived in Section II.1. To discard the burn-in phase, we remove the first 75% of the steps from each chain and visually check that the chains remain stable thereafter.

III.1.1 Parametrizations & Priors of JamPy

Our parametrization follows the approach of Paper I. The baryonic content is fixed and modeled with a flattened Plummer profile, with parameters determined in Section II.1. Owing to the lack of external constraints on the halo shape, the DM halo is modeled as a spherical generalized Plummer profile,252525A special case of the $\alpha\beta\gamma$ profile (Zhao, 1996), with $\alpha=2$, $\beta=5$, and a free inner slope $\gamma$ (see eq. A3 from Paper I). with a steeper outer slope ($\beta=5$) motivated by tidal interactions inferred from Sculptor’s orbital history (Sohn et al., 2017), consistent with the findings of Peñarrubia et al. (2010). This DM parametrization was also adopted in Paper I, Hayashi et al. (2020), and Goldstein & Strigari (2025).

This density model depends on three parameters: a scale radius $r_{\rm dark}$, a total DM mass ($M_{\rm dark}$), and an inner slope $\gamma_{\rm dark}$. We assume flat/log-flat priors for all of these variables,

• •

$\gamma_{\rm dark}\in[-2,2]$, which includes both cuspy ($\gamma=-1$) and cored ($\gamma=0$) profiles. We allow for positive slopes to account for potential physical mechanisms not necessarily anticipated by $\Lambda$CDM.

• •

$\log{(M_{\rm dark}/[\rm M_{\odot}])}\in[6,12]$. Read et al. (2017) extrapolated classical $M_{\star}-M_{200}$ relations to lower-mass dSphs, such that Sculptor, with a luminous mass of $\sim 8.08\times 10^{6}~\rm M_{\odot}$, should have $M_{200}\sim 10^{9}~\rm M_{\odot}$. Therefore, our priors largely cover this range.

• •

$\log{(r_{\rm dark}/[{\rm kpc}])}\in[-1,2]$. Given the expected $M_{200}$ mass from Read et al. (2017),262626We thank Justin Read for sharing his algorithm to compute the precise $M_{200}$ value expected for Sculptor. the concentration relation from Dutton & Macciò (2014) yields a Navarro et al. (1997) profile scale radius of nearly $2$ kpc. Our prior thus covers this with more than an order of magnitude range.

We assume that the anisotropy parameters, $\beta_{\rm J}$ and $\Omega$, are constant and fit their symmetrized quantity,272727$x_{\rm sym}$ ranges monotonically from $-2$ to $+2$, given that $x$ spans similarly from $-\infty$ to $+1$. defined as $x_{\rm sym}=x/(1-x/2)$, using a flat prior between $-1.99$ and $1.99$. To assess the impact of this assumption, we have separately tested variations with a radial dependence, such as generalized parametrizations (Osipkov, 1979; Merritt, 1985), and found no significant effect on our results and conclusions presented further on.

Finally, we apply Gaussian priors for the bulk $v_{\rm LOS}$ of Sculptor, setting $v_{\rm LOS}=111.2\pm 0.3~\rm km\,s^{-1}$ (Arroyo-Polonio et al., 2024), and for the distance modulus, defined as $\mu_{0}=5\,\log{(D/[{\rm kpc}])}+10$. We adopt $\mu_{0}=19.67\pm 0.05$, which is the weighted mean of measurements reported in (Rizzi, 2002, RGB and TRGB), (Pietrzyński et al., 2008, RR-Lyrae on V and I bands), and (Górski et al., 2011, TRGB on J and K bands).

III.2 Practical Quantities

As in Paper I, we compute derived quantities to better interpret the results of our fits, as some quantities can be difficult to assess in their original form. For instance, due to the mass-anisotropy degeneracy, constraining the asymptotic DM density slope is particularly challenging in the absence of 3D velocities near the galaxy center. Instead, a more meaningful approach is to evaluate the DM density slope within the radial range where both PMs and LOS velocities are available. We define this as an effective DM slope:

[TABLE]

where $\rho(r)$ is the DM density. The lower bound, $r_{\rm min}$, corresponds to the smallest projected radius in the data where PM information is available. We set the upper bound as $r_{\rm max}={\rm min}(R_{\rm max,PM},r_{\Lambda{\rm CDM}}/3)$, where $R_{\rm max,PM}$ is the largest projected radius with available PM data, and $r_{\Lambda{\rm CDM}}$ is the NFW scale radius expected for DM halos in $\Lambda$CDM simulations.282828We adopt this reference value to facilitate comparison with theoretical predictions. The factor $(1/3)$ is introduced to exclude the region where the density profile steepens due to the transition from the inner to outer profile. (computed with the relations from Read et al. 2017).

Additionally, as in our previous work, we compute the circular velocity:

[TABLE]

which depends on the total292929The total mass includes both luminous and dark components. cumulative mass enclosed within radius $r$ and the gravitational constant $G$.

IV Results

We ran JamPy over a set of eleven inclinations, linearly spaced from $43\degree.7$ to $90\degree$ (edge-on), covering the range of Sculptor’s inclination PDF computed in Section II.1, to mimic the analyses from Paper I. The results are presented in Table 2. In addition to the anisotropy parameter $\beta_{\rm J}$, defined in Eq.(4), we also report the globally-averaged Binney anisotropy, $\overline{\beta_{\rm B}}$. The latter is computed by first determining the stellar-mass-weighted second velocity moments, averaged over the entire system, and then substituting them into Eq.(5).

IV.1 Inclination Dependence

Regarding the dependence of model parameters on the adopted galaxy inclination, Paper I found that, while correlations were present, Draco’s overall parameters remained relatively consistent. That is, the parameters largely retained their signs, preserving the general regime they indicated (e.g., velocity anisotropy remained tangential, and the DM slope differed from the cored case across all inclinations).

In this study, we observe similar correlations with inclination; however, their impact on our modeling conclusions is much more pronounced. Specifically, the best-fitting flattened systems exhibit extreme radial anisotropy with a corresponding steep density drop, whereas the best-fitting edge-on systems exhibit extreme tangential anisotropy with a corresponding highly-cuspy density profile. This diversity of DM slopes is showcased in Figure 9.

The inclination-averaged best-fit parameters of our model were computed using the inclination PDF shown in Figure 2, following the formalism presented in Paper I, Section 4.2.1. As expected from the strong degeneracy discussed above, the DM slope and velocity anisotropy are well constrained only for specific inclination values, but not in the averaged case. This result highlights the critical role of properly accounting for inclination in the mass modeling of flattened systems.

IV.1.1 Higher-order velocity moments

When first-order velocity moments are available in all three dimensions, the inclination of an axisymmetric system can be directly constrained using relations like those in Evans & de Zeeuw (1994, eq. [A4]) or via Jeans modeling with the inclination as a free parameter. However, correcting systematic effects in HST-based proper motions (as detailed in Section II.3.2) removes transverse first-order moments by design. Our Jeans analysis thus relies solely on the LOS first-order moment, which is insufficient to determine inclination on its own.

To address this degeneracy, we turn to the scalefree models of de Bruijne et al. (1996), further detailed in Paper I (appendix B). These models assume an axisymmetric baryonic component in a spherical potential303030This ensures consistency with our current Jeans setup. and, by adopting fiducial distribution functions and potential forms, offer a reasonable and generalized approximation of systems like dSphs with flattened dispersion profiles. While these models do not capture the detailed radial variation of the logarithmic density and potential slope, they can, unlike Jeans modeling, use higher-order moments of a single velocity dimension (here, LOS) to break degeneracies between mass, anisotropy, and crucially, inclination. This offers educated guesses of the inclination via features such as the kurtosis of the LOS velocity distribution (e.g., Paper I, fig. 18).

Following this approach, we compute the observed LOS kurtosis to be $\overline{K_{\rm LOS}}=0.63\pm 0.05$ and compare it to scalefree predictions across a range of inclinations, using our inferred velocity anisotropies. Figure 10 shows two panels, one per scalefree distribution function type, with the observed kurtosis distribution color-coded as in the leftmost color bar, and its mean marked by a dashed black line. For each inclination, we compute the expected kurtosis from scalefree using three representative values of the average velocity anisotropy $\overline{\beta_{\rm B}}$: the 16th percentile, the median, and the 84th percentile (from Table 2, column 9).

In both distribution function cases, the models favor lower inclinations (i.e., more flattened systems). The first panel yields the best match313131Defined by the most evident overlap between the three anisotropy curves and the dashed black line. at $i=57\fdg 1$, with the second at $i=52\fdg 4$. Although these generalized scalefree-based models are not designed for precise inclination fitting, the consistent need for radial velocity anisotropy to match the observed kurtosis supports the idea that Sculptor lies toward the more flattened end of the viable axial ratio range. This result remains robust even under conservative 3- and 4-$\sigma$ clipping of the LOS velocity distribution.

IV.2 Data-Model Comparison

While having in mind the degeneracies in our dataset and the inclination-averaged estimates summarized in the last row of Table 2, adopting a Jeans model tied to a fiducial inclination offers a clearer benchmark for comparing data and interpreting results. For this purpose, we adopt the $i=57\fdg 1$ Jeans model as our fiducial setting, as it corresponds both to the best scalefree match under its distribution function of type 1 and approximates the median inclination of the posterior shown in Figure 2.

Figure 8 illustrates the goodness-of-fit of the radial kinematic profiles from our Jeans models. Black points denote binned data: LOS kinematics on the left – annularly-averaged velocity dispersion (top) and major-axis rotation (bottom); and PM kinematics on the right – velocity dispersions in the radial (top) and tangential (bottom) directions. The solid black lines show the median model, while shaded regions represent MCMC percentiles, using the same fiducial inclination of $57\fdg 1$ to convert PMs to physical velocities. Model curves are interpolated over projected radius $R$ for display, although they also vary with projected angle. The model provides a satisfactory overall agreement, consistent with fits at other inclinations.

IV.3 Inferred Quantities

Degeneracies between inclination and parameters such as the DM inner slope and the velocity anisotropy profile render the inclination-averaged estimates of these quantities inconclusive, spanning a range of scenarios consistent with various cosmological models and DM candidates. However, incorporating higher-order moments, as outlined in Section IV.1.1, allows us to favor more flattened or less inclined configurations in the Jeans modeling. In this context, confidence in shallower DM slopes becomes notable: for our fiducial inclination of $i=57\fdg 1$, we rule out $\Gamma_{\rm dark}\leq-1$ at 99.8% confidence, while the inferred velocity anisotropy is mildly radial but still statistically consistent with isotropy. These trends are clearly reflected in the corner plots shown in Figure 11. Respectively, the confidence in such general patterns decreases or increases for more or less inclined systems.

The same trend is illustrated in the upper panel of Figure 12. Over the radial range spanned by the 3D-velocity data (i.e. the limits of the respective $x$ axis) for the $i=57\fdg 1$ model, the DM density profile clearly favors a shallow inner slope. As noted earlier, edge-on configurations yield more cuspy profiles, though they are disfavored due to the higher LOS velocity kurtosis in the data. The same figure also presents Sculptor’s circular velocity profile for the fiducial inclination. Here, the velocity peaks at $\mathrm{max}(v_{\rm circ})=24.18^{+1.29}_{-1.19}~\rm km\,s^{-1}$, followed by a decline within the data extent – a feature not captured for Draco in Paper I, and in good agreement with the similarly non-spherical analysis of Hayashi et al. (2020, figure 9).

The DM mass enclosed within the data extent for the inclination-averaged fit – similar to values at other inclinations – is $M_{\rm dark}^{R_{\rm max}}=1.39^{+0.68}_{-0.58}\times 10^{8}~\rm M_{\odot}$, with an extrapolated total mass of $M_{\rm dark}^{\rm total}=1.58^{+1.01}_{-0.62}\times 10^{8}~\rm M_{\odot}$. This is far below the $\sim 10^{10}~\rm M_{\odot}$ inferred by Errani, Peñarrubia, & Walker (2018) for cored profiles, together with their $r_{\rm max}\sim 10$ kpc and $v_{\rm max}\sim 65~\rm km\,s^{-1}$. The difference is likely explained as follows: their method matches $M(<1.8\,R_{\rm h})$ – a quantity consistent with ours323232Errani et al. (2018) predict $M(<1.8\,R_{\rm h})=3.2\pm 0.7\times 10^{7}~\rm M_{\odot}$, while our inclination averaged measurement is $M(<1.8\,R_{\rm h})=5\pm 2\times 10^{7}~\rm M_{\odot}$. – to subhalos in the Aquarius simulations (Springel et al., 2008), but this does not guarantee agreement in the extrapolated outer halo. Because their approach links the halo outskirts directly to the mass within $1.8\,R_{\rm h}$, cusp and core models can produce very different $M_{\rm total}$, $r_{\rm max}$, and $v_{\rm max}$. Our Jeans-based fits, in contrast, allow the outer profile to vary independently of correlations inherent to cosmological simulations, although our choice of a steeper outer slope333333$\beta=5$ in the generalized $\alpha\beta\gamma$ model. naturally yields lower total masses, even when matching at smaller radii. If Sculptor’s halo was once more massive and lost substantial DM to tidal interactions (see Section V.1), this steep slope remains physically justified (Peñarrubia et al., 2010).

Beyond inclination-sensitive parameters, some constraints remain tight despite their correlations with inclination. This is the case for our inclination-averaged distance estimate, $D=85.83^{+2.02}_{-1.91}$ kpc. This measurement is valuable in its own right, as dynamical distances provide an independent benchmark for stellar evolution models and parallax-based methods. The following sections assess the robustness of these results and their broader implications for cosmology and equilibrium-based mass modeling.

V Robustness of Modeling Results

V.1 The Effect of Tides

Perhaps the most fundamental assumption of our model is that the Sculptor dSph is in dynamical equilibrium. However, this galaxy orbits a significantly more massive system (i.e. the Milky Way) and is therefore subject to tidal interactions that could perturb its equilibrium state. This possibility has been raised in previous studies to account for the inflated velocity dispersion profiles observed in dSphs, which appear larger than expected from their baryonic content alone (Klessen & Kroupa, 1998; Hammer et al., 2018; Wang et al., 2024). Such arguments have even been used to question the necessity of DM in these systems, although the existence of satellites exhibiting high velocity dispersions while orbiting at relatively large pericentric distances (e.g. Li et al., 2021; Battaglia et al., 2022; Pace et al., 2022) remains a strong counterargument to this interpretation.

In our case, the total dark mass inferred for Sculptor is $1.58\times 10^{8}~\rm M_{\odot}$, nearly an order of magnitude lower than the $4.85\times 10^{9}~\rm M_{\odot}$ predicted by $\Lambda$CDM scaling relations for low-mass dwarfs (Read et al., 2017). This discrepancy could be somewhat accounted if a large fraction of Sculptor’s DM halo has been stripped by the Milky Way’s tidal field over its lifetime. Such a scenario is consistent with orbital reconstructions (e.g., Sohn et al., 2017), which indicate several close pericentric passages around the Milky Way, that would have favored substantial mass loss.

Sculptor’s dynamical stability under such conditions was examined by Iorio et al. (2019), who ran simulations to evaluate how tidal interactions might affect the system and the validity of equilibrium-based mass estimates. Even under extreme orbital scenarios permitted by current observed data, they found little evidence for significant tidal disturbance in the stellar distribution. While the DM halo may undergo stripping, the stars themselves remain largely unaffected, showing only gradual changes over time. In the last 2 Gyr, they found these variations to be minimal, well within a few percent. Such findings, alongside similar conclusions from more recent studies (e.g., Tchiorniy & Genina, 2025), support the use of equilibrium-based models for estimating Sculptor’s inner mass profile.

V.2 Tests with mock data

As a final robustness check, we apply our modeling strategies to mock datasets to evaluate both their performance and the influence of key assumptions. We focus on two primary tests: First, we simulate a system that adheres to the assumptions of our modeling framework, and examine whether elements such as the use of Gaussian likelihoods for projected velocities or our data-cleaning procedures – particularly the local PM corrections applied to mitigate systematics – introduce significant biases in the inferred parameters. Second, we assess the impact of multiple kinematic populations within Sculptor, despite our modeling being based on a single-component assumption.

To conduct these tests, we employ the distribution function-based code Agama (Vasiliev, 2019) to generate spherical343434At present, Agama does not support fully axisymmetric distributions with a given anisotropy and density profile, as would be ideal for comparison with the observed data. stellar systems in equilibrium, given a potential and mass distribution. Each mock setup is described below, while the general procedure used to align the mocks with the observed data is detailed in Appendix A.

V.2.1 Data and Likelihood Handling

For the first test, we construct a rotating mock system consisting of an isotropic Plummer stellar component with the same total mass and scale radius as Sculptor (see Appendix A for details). The DM halo is modeled using a cuspy generalized Plummer profile (i.e., $\gamma=-1$), with the total mass and scale radius set according to the expectations from $\Lambda$CDM cosmology, as outlined in Section III.1.1.

By comparing our fits to this simulated data to the true input values, we can assess the impact of our local PM corrections to eliminate HST systematic effects, as well as evaluate the performance of our fitting strategy. Figure 13 shows this comparison for a selection of the most relevant fit parameters: the slope of the DM profile over the region where 3D velocities are available, $\Gamma_{\rm dark}$, and the globally-averaged $\overline{\beta_{\rm B}}$ (cf. Eq. 5). The close agreement between the smaller black markers (true values) and the larger red markers (estimated values) highlights the proper handling of the data and demonstrates satisfactory convergence of our fit, assuming the parametrizations are correct. Additionally, the posterior distributions in blue exhibit similar widths to those obtained from fitting the true data, which further supports the reliability of our error bars presented in Table 2.

V.2.2 Distinct kinematic populations

One of the key outcomes of spectroscopic campaigns targeting dSphs, with Sculptor frequently serving as a benchmark system, has been the identification of at least two distinct chemo-dynamical populations (e.g. Tolstoy et al., 2004). These subpopulations have proven instrumental in mitigating degeneracies inherent to mass modeling based solely on LOS velocity measurements (e.g. Walker & Peñarrubia, 2011). More recently, expanded observational efforts (Tolstoy et al., 2023; Barbosa et al., 2025) have substantially advanced the characterization of Sculptor’s internal kinematics, yielding improved constraints on the spatial distributions, membership probabilities, and LOS velocity dispersion profiles of its distinct stellar populations (e.g. Arroyo-Polonio et al., 2024).

This naturally raises the question of whether modeling the system as a single dynamical component introduces bias – particularly in the estimation of key parameters such as $\Gamma_{\rm dark}$ and the mean $\beta_{\rm B}$ for the data. To assess the potential impact of multiple stellar populations, we simulate a composite system comprising metal-poor and metal-rich components. These components are designed to approximate the spatial and kinematic structure inferred by Arroyo-Polonio et al. (2025) for Sculptor, as illustrated in their figure 1, but without rotation, as not to mix different test diagnostics. To implement this, we extracted the relevant data from that figure using online plot digitization tools,353535Specifically, https://web.eecs.utk.edu/~dcostine/personal/PowerDeviceLib/DigiTest/index.html. and fitted each population: we used a Plummer profile for the surface density and constant velocity anisotropy model for their inferred $\beta_{\rm B}$.363636As discussed in Section III.1.1, separate tests showed that assuming a radially varying anisotropy did not significantly alter the results presented in Table 2, thereby justifying the simpler assumption of a constant value per population. Moreover, the large $\beta_{\rm B}$ uncertainties in figure 1 of Arroyo-Polonio et al. (2025) are statistically consistent with this simplified treatment. We then assigned tracer counts based on the normalization of each fitted density profile.

Figure 14 presents the results of our fits to this chemo-dynamically distinct mock dataset, focusing on the same key parameters evaluated in the previous test. Overall, we find that the estimates remain consistent within 1-$\sigma$ uncertainties, albeit with slightly reduced performance compared to the single-population case. Notably, the posterior distribution of the $\Gamma_{\rm dark}$ parameter exhibits a different shape and central value, suggesting that neglecting multiple populations may slightly inflate the formal uncertainties. This behavior implies a potential bias in precision and accuracy. In summary, while the simplified single-population modeling does not introduce substantial bias in the recovery of the DM density profile, it is not entirely neutral and may slightly favor shallower inner slopes relative to Sculptor’s true underlying profile. Meanwhile, the average velocity anisotropy remains relatively well constrained, benefiting from the full 3D velocity information.

VI Discussion

VI.1 Implications for mass modeling

Nearby dSphs provide a rare opportunity to test small-scale predictions of cosmological models using resolved stellar kinematics. This has driven extensive LOS velocity surveys with ground-based telescopes (e.g., Tolstoy et al., 2004; Walker et al., 2007; Gilmore et al., 2022), aimed at constraining the underlying DM distributions. Many mass-modeling studies, however, have relied on simplifying assumptions: most notably, spherical symmetry (e.g., Battaglia et al., 2008; Walker & Peñarrubia, 2011; Read et al., 2018; Arroyo-Polonio et al., 2025), despite clear observational evidence that dSphs are flattened systems (e.g., Muñoz et al., 2018). Even with growing access to more complex PM data (e.g., Massari et al., 2017, 2020), few works have incorporated intrinsic geometry into their models (see Hayashi et al. 2020 and Pace et al. 2020).

Our results, consistent with Paper I, demonstrate that accounting for dSph flattening significantly alters the inferred DM density profile and stellar orbital structure. While Paper I emphasized how the projected velocity anisotropy varies with position angle, potentially biasing DM slope estimates under spherical symmetry,373737In spherical models, observables are assumed to be identical at a given radius regardless of position angle. our current analysis highlights a strong dependence on the galaxy’s inclination. This introduces a key degeneracy overlooked by spherical, non-rotating models. As shown in Table 2, several modeled parameters shift with inclination, in some cases beyond formal uncertainties.

These findings support the case for more flexible modeling frameworks, even those possibly including non-spherical or triaxial DM halos and baryonic content. Yet, such approaches remain difficult in practice due to the lack of independent constraints on 3D shapes for individual dSphs. Axisymmetric models offer a pragmatic middle ground, requiring only projected flattening and inclination, with the former already well measured. Independent inclination constraints may soon be accessible through future Gaia data releases, which will deliver higher precision maps of mean PM, hopefully with higher signal than Gaia’s systematic uncertainties. Combining these with LOS velocity data via established formalisms (e.g., Evans & de Zeeuw, 1994; Bianchini et al., 2018) may finally enable robust, model-independent estimates of dSph inclinations.

Another important complexity – often not explicitly modeled, including in this study – is the presence of multiple chemically-distinct stellar populations within dSphs. Many systems are known to host at least two such populations, typically metal-poor and metal-rich, which exhibit different LOS velocity dispersions (e.g., Pace et al., 2020; Arroyo-Polonio et al., 2024). Whether these populations require separate anisotropy prescriptions in Jeans-based modeling remains an open question, with only a limited number of systems studied in detail (e.g., Arroyo-Polonio et al., 2025). In the specific case of Sculptor, our mock data tests suggest that treating its stars as a single dynamical component does not lead to significant bias in the recovered DM density profile, though it may introduce a modest loss of accuracy within the formal 1-$\sigma$ uncertainty range.383838Besides, the locations of the 16th and 84th percentiles may change, while their relative distance to each other seems somewhat the same as the one from the single population test.

This underscores the advantages of better associating PMs with distinct chemical populations to fully exploit 3D kinematics. Progress will depend on joint astrometric and spectroscopic observations, especially at faint magnitudes where most dSph stars lie and where facilities like HST, JWST, Euclid (e.g. Libralato et al., 2024) and the upcoming Nancy Grace Roman Space Telescope will be crucial.

VI.2 Implications for the nature of dark matter

As discussed in Sections IV.1 and IV.3, significant degeneracies arise due to Sculptor’s unknown inclination, preventing firm constraints on the nature of DM across the full range of tested configurations. These limitations stem from both the data, such as the lack of odd PM moments by design, and our modeling approach, which in the Jeans formalism does not account for velocity moments beyond second order. However, as shown in Section IV.1.1 and in Section 5.4 of Paper I, higher-order moments of the LOS velocities provide additional information on velocity anisotropy and the galaxy’s shape, helping to alleviate these degeneracies (e.g. Read & Steger, 2017; Read et al., 2018). While our analysis does not yield a precise estimate of Sculptor’s inclination, comparisons with scalefree models from de Bruijne et al. (1996) suggest that Sculptor likely favors a more flattened and less inclined configuration, deviating from an edge-on geometry.

At low inclination, we find that the DM slope within the PM-covered region is constrained toward shallower values, deviating from the cuspy profiles predicted by $\Lambda$CDM DM-only simulations (Navarro et al., 1997). This result supports previous findings that favor core-like DM profiles in dSphs, based on stellar kinematics (e.g. Battaglia et al., 2008; Walker & Peñarrubia, 2011; Breddels et al., 2013; Zhu et al., 2016; Amorisco & Evans, 2012; Brownsberger & Randall, 2021; Arroyo-Polonio et al., 2025) and HI rotation curves (e.g. Moore, 1994; Flores & Primack, 1994; Burkert, 1995). Within $\Lambda$CDM, the most common explanation involves the influence of baryonic processes, which can transform central cusps into cores by redistributing mass and energy via supernova feedback (Read & Gilmore, 2005; Pontzen & Governato, 2012; Brooks & Zolotov, 2014) or bursts of star formation (Read et al., 2018). Indeed, Fitts et al. (2017) identify a threshold around $M_{\star}\approx 2\times 10^{6}~\rm M_{\odot}$, above which these effects become efficient.

The comparison of total stellar masses between Draco ($4.7\times 10^{5}~\rm M_{\odot}$, Martin et al. 2008, Paper I), where we found a clear preference for a cusp, and Sculptor ($8.08\times 10^{6}~\rm M_{\odot}$, de Boer et al. 2012, this work), where our analyses seem to favor a core, qualitatively supports this picture and presents no direct challenge to $\Lambda$CDM. Nonetheless, underlying quantitative details remain uncertain. For instance, is there a predictable analytical or empirical link between the degree of DM slope flattening (i.e., the value of $\gamma_{\rm dark}$) and a galaxy’s baryonic content? How does the host halo mass, fixed at $M_{\rm halo}\approx 10^{10}~\rm M_{\odot}$ in Fitts et al. (2017), modulate this process? And how long would it take for a system like Sculptor to regenerate a cuspy density profile from the apparent present-day core? (e.g. Laporte & Penarrubia, 2015) These questions point to the need for more nuanced exploration of stellar feedback in cosmological simulations to reproduce the configurations favored by our fits. In parallel, our results serve as valuable observational constraints for testing both $\Lambda$CDM and alternative models in which core formation arises naturally, such as self-interacting DM (SIDM, e.g. Nadler, Yang, & Yu 2023) with larger cross sections or Fuzzy DM with specific particle mass ranges (Luu et al., 2025).

VII Conclusions

We present a new proper motion catalog for the Sculptor dwarf spheroidal galaxy, comprising 119 stars with proper motion uncertainties smaller than the galaxy’s intrinsic velocity dispersion (i.e. $\sim 10~\rm km\,s^{-1}$) – making it the most precise proper motion dataset compiled for this system to date. We also combined line-of-sight velocities from multiple literature sources, resulting in a sample of 1,760 stars with similarly low velocity uncertainties. Using this state-of-the-art public dataset,393939We will make the data public upon paper acceptance. we conducted axisymmetric Jeans-based mass and anisotropy modeling, complemented by robustness tests involving higher-order velocity moments and mock datasets. Below, we summarize the main conclusions of this work.

• •

We present updated estimates for Sculptor’s center, ellipticity, position angle, and scale radius, based on fits to stellar counts from Gaia EDR3 (Table 1).

• •

We confirm mild oblate rotation, with major-axis amplitudes reaching up to $\sim 2~\rm km\,s^{-1}$ beyond $20\farcm 0$ from the galaxy center, consistent with recent line-of-sight based findings (Arroyo-Polonio et al., 2024).

• •

We demonstrate that unresolved binaries have a negligible effect on the line-of-sight velocity dispersion and thus do not significantly impact the dynamical modeling.

• •

From proper motions measured in our HST fields along the projected minor axis, we find average observed velocity dispersion ratios of $\left<\sigma_{\rm POSt}\right>/\left<\sigma_{\rm POSr}\right>=1.19\pm 0.19$ and $\left<\sigma_{\rm LOS}\right>/\left<\sigma_{\rm POS}\right>=0.93\pm 0.08$, where $\sigma_{\rm POS}$ is the average dispersion across both proper motion directions. The first ratio is independent of distance, while the second scales inversely with it.404040We assumed here $D=85.83$ kpc.

• •

We perform tests with mock datasets constructed using Agama to confirm that our current treatment of HST proper motion systematics does not introduce significant biases in the fitted parameters of our Jeans modeling. Additionally, for the specific case of Sculptor, we show that modeling its stellar component as a single dynamical population does not significantly bias the recovered DM density profile, though it may introduce a modest reduction in precision and accuracy within the formal 1-$\sigma$ uncertainty range.

• •

Our non-spherical modeling reveals a significant degeneracy tied to the unknown inclination of the galaxy – an effect missed under spherical symmetry assumptions. This degeneracy allows acceptable fits across a wide range of dark matter profiles, from cuspy to cored and even centrally-decreasing density configurations. While we do not directly constrain the inclination, higher-order line-of-sight velocity moments provide useful limits, and comparisons with scalefree models from de Bruijne et al. (1996) favor highly flattened (i.e., less inclined) configurations.

• •

For a representative low inclination of $i=57\fdg 1$, near the median of the inferred inclination PDF (Figure 2), we find, alongside well-constrained radial velocity anisotropy, a dark matter density slope of $\Gamma_{\rm dark}=0.29^{+0.31}_{-0.41}$ within the radial extent of our 3D velocity data, ruling out $\Gamma_{\rm dark}\leq-1$ at 99.8% confidence.414141An asymptotic dark matter density slope of $-1$ represents the classic inference from dark matter only $\Lambda$CDM simulations (Navarro et al., 1997) This confidence level increases with lower inclinations and decreases for close to edge-on configurations. At the same inclination, we infer a maximum circular velocity of $\mathrm{max}(v_{\rm circ})=24.18^{+1.29}_{-1.19}~\rm km\,s^{-1}$.

• •

At lower inclinations, our results qualitatively agree with $\Lambda$CDM predictions in which feedback in higher-mass dwarf galaxies produces shallower dark matter profiles (Fitts et al., 2017). These models predict a transition from cusp to core at $M_{\star}\approx 2\times 10^{6}~\rm M_{\odot}$, while Sculptor’s total stellar mass is $8.08\times 10^{6}~\rm M_{\odot}$ (de Boer et al., 2012). Likewise, alternative DM models such as SIDM and Fuzzy DM also predict core formation for specific ranges of cross sections and particle masses, respectively, making our measurements a useful means of narrowing these parameter spaces.

• •

Averaged over all inclinations, we estimate an enclosed mass within the radial extent of our data of $M_{\rm dark}^{R_{\rm max}}=1.39^{+0.68}_{-0.58}\times 10^{8}~\rm M_{\odot}$, and a total dark matter mass of $M_{\rm dark}^{\rm total}=1.58^{+1.01}_{-0.62}\times 10^{8}~\rm M_{\odot}$. The modest increase beyond $R_{\rm max}$ reflects our choice of a steep outer halo slope ($\beta=5$), consistent with expectations for tidally affected satellites (Peñarrubia et al., 2010). We also constrain a dynamical distance of $D=85.83^{+2.02}_{-1.91}$ kpc.

Following the conclusions from the first paper in this series (Paper I), we reaffirm the growing potential for a deeper understanding of dwarf spheroidal galaxy internal kinematics. This progress relies on assembling long time-baseline proper motion datasets with the most precise astrometric facilities, such as HST, JWST, and Gaia, and on the development of future missions like Roman. Ongoing work includes observations of additional systems (e.g. Vitral et al., 2023b) and the refinement of HST-based measurements using JWST (e.g. van der Marel et al., 2023; Bennet et al., 2025a, b).424242These proposals can be retrieved in the following links:

https://ui.adsabs.harvard.edu/abs/2023hst..prop17434V/.

https://ui.adsabs.harvard.edu/abs/2023jwst.prop.4513V/.

https://ui.adsabs.harvard.edu/abs/2025hst..prop17926B/.

https://ui.adsabs.harvard.edu/abs/2025jwst.prop.9225B/.

With such data, significant advances in galactic dynamics and dark matter studies are within reach.

VIII Acknowledgments

We thank the anonymous referee for providing constructive comments that improved the quality of the paper. Eduardo Vitral acknowledges funding from the Royal Society, under the Newton International Fellowship programme (NIF\R1\241973). Support for this work was provided by NASA through grants for programs GO-12966 and GO-16737 from the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under NASA contract NAS5-26555, and through funding to the JWST Telescope Scientist Team (PI: M. Mountain) through grant 80NSSC20K0586. We thank Justin Read and Thomas de Boer, for sharing the algorithms used in Read et al. (2017) and de Boer et al. (2012), respectively.

The Digitized Sky Surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This project is part of the HSTPROMO (High-resolution Space Telescope PROper MOtion) Collaboration (https://www.stsci.edu/~marel/hstpromo.html), a set of projects aimed at improving our dynamical understanding of stars, clusters, and galaxies in the nearby Universe through measurement and interpretation of PMs from HST, Gaia, JWST, and other space observatories. We thank the collaboration members for the sharing of their ideas and software.

Data availability: Some of the data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via https://doi.org/10.17909/chzv-p290 (catalog DOI: 10.17909/chzv-p290). See Appendix B for the data products of this work.

Appendix A Mock data recipe

In this section, we outline the procedure used to construct realistic mock datasets from the raw output generated by the Agama distribution function-based software (Vasiliev, 2019). This output provides phase-space information for a given kinematical population, with stars distributed according to a specified input density profile within a given potential. For context, we adopted Agama’s QuasiSpherical distribution function,434343Details of this module are available in the software documentation: https://github.com/GalacticDynamics-Oxford/Agama. assuming constant velocity anisotropy. This was set to isotropic in mock test #1 (cf. Section V.2.1), and matched to the mean anisotropies inferred for Sculptor’s metal-poor and metal-rich populations from figure 1 of Arroyo-Polonio et al. (2025) in mock test #2 (cf. Section V.2.2). In mock test #1, we also introduced a mild rotational component by flipping the velocity vectors of 20% of the stars444444This choice reflects our findings from Section II.2.2, where a similarly low rotation amplitude was observed. with $L_{z}<0$. To ensure that this rotation would produce observable signatures in the plane of the sky, we adopted an intermediate inclination of $i=57\fdg 1$, deviating from an edge-on configuration where PM streaming would become undetectable. This setup allowed us to reproduce a scenario broadly consistent with the inclination constraints from Section IV.1.1, while still enabling a meaningful test of our artificial removal of streaming motions in the PM data.

The first step once the Agama outcome was ready involved converting the physical quantities of velocities and sky positions into observed, angular-based quantities, assuming a specific distance to the galaxy. We adopted the same distance value used in the priors for our Jeans modeling. Next, we applied local corrections to the PM components of the mock data, following the same procedure outlined in Section II.3.2 and in Paper I. This correction effectively removed any trace of transverse streaming motions from the mock data.

The next step was to select stars located within the same sky region as the real data for each subset (i.e., the LOS and PM subsets). This was accomplished by first identifying the stars within a convex hull defined by the spatial boundaries of each observed dataset. We then randomly selected the same number of stars from these regions as were present in the observations, aiming to assess the influence of spatial incompleteness on our estimates. Figure 15 illustrates the spatial distribution of the stars in each of the three components: the original mock data sampled by Agama (shown in green), the selected LOS subset (in blue), and the PM subset (in red).

The final step was to assign realistic uncertainties to both the $v_{\rm LOS}$ and PM quantities. For the LOS velocities, we sampled errors from the 1D distribution of observed $\epsilon_{v_{\rm LOS}}$. For the PMs, we drew from the observed 2D distribution of $\epsilon_{\mu_{\alpha,*}}$ vs. $\epsilon_{\mu_{\delta}}$, which reflects both nominal measurement uncertainties and the additional scatter introduced by local systematic corrections in the real dataset. Although this procedure does not exactly replicate the real-data reduction pipeline – where local corrections are applied after centering-related errors are already present – it provides an effective approximation. Under the assumption of uncorrelated, Gaussian-distributed errors (which we tested in Section II.3.2), the statistical properties of the final error distribution remain robust to the order in which uncertainties are introduced.454545This is a consequence of the fact that a convolution of two Gaussian distributions is still a Gaussian. Moreover, drawing directly from the real observed uncertainties ensures that the mock error budget reflects realistic observational effects such as CTE degradation, which are absent in the idealized Agama mocks. We then added zero-centered Gaussian noise with these sampled standard deviations to the respective $v_{\rm LOS}$ and PM values, completing the construction of the mock dataset.

Appendix B Online material

In this section, we briefly describe the kinematic datasets used in this study, which are provided as online material. These include a proper-motion table (Table 3) derived from the HST observations presented here, and a line-of-sight velocity table (Table 4) assembled by cross-matching the original catalogs of Tolstoy et al. (2023) and Walker et al. (2023). For further details on how the data was assembled, see Sections II.2 and II.3.