Dynamical models for the Sculptor dwarf spheroidal in a CDM universe
Dynamical models for the Sculptor dwarf spheroidal in a CDM universe
Abstract
The Sculptor dwarf spheroidal galaxy appears to contain two distinct stellar populations of differing metallicity. Several authors have argued that in order for these two populations to reside in the same gravitational potential, the dark matter halo must have a core similar to that observed in the stellar count profile. This would exclude cuspy NavarroFrenkWhite (NFW) density profiles of the kind predicted for halos and subhalos by dark matter only simulations of the CDM cosmological model. We present a new theoretical framework to analyze stellar count and velocity observations in a selfconsistent manner based on separable models, , for the distribution function of an equilibrium spherical system. We use this machinery to analyze available photometric and kinematic data for the two stellar populations in Sculptor. We find, contrary to some previous claims, that the data are consistent with populations in equilibrium within an NFW dark matter potential with structural parameters in the range expected in CDM; we find no statistical preference for a potential with a core. Our models allow a maximum circular velocity for Sculptor between 20 and 35 km/s. We discuss why some previous authors came to a different conclusion.
1. Introduction
CDM has emerged as the standard model of cosmic structure largely because of its successful predictions for the temperature anisotropies of the cosmic microwave background radiation, the power spectrum of the largescale distribution of galaxies, the structure of the main body of galaxy and galaxy cluster halos, as inferred from weak gravitational lensing, and the broad features of the galaxy formation process (see Frenk & White 2012, for a review). However, the adequacy of the model remains controversial in the well observed inner regions of galaxies where the distribution of dark matter is strongly nonlinear (e.g. Walker & Penarrubia 2011; Newman et al. 2013). Disagreements on these scale are particularly interesting because they could provide clues to the nature of the dark matter (e.g. Lovell et al. 2012; Peter et al. 2013; Shao et al. 2013; Zavala et al. 2013).
A robust prediction of CDM is that, in the absence of baryonic effects, the spherically averaged radial density profiles of dark matter halos of all masses should approximately follow a universal form, the NFW profile (Navarro et al. 1996, 1997), which diverges as towards the centre. In galactic halos, baryonic effects associated with the formation of the galaxy could, in principle, flatten this central cusp through explosive events produced by supernovae, as proposed by Navarro et al. (1996) and seen more recently in a number of galaxy formation simulations (see Pontzen & Governato 2014, for a review). Energetic arguments suggest that such processes  if they do indeed occur in nature  should be ineffective in dwarf galaxies of stellar mass below or M which would then retain their NFW dark matter cusps (Peñarrubia et al. 2012). Although difficult to study because of their intrinsic faintness, the lowest mass dwarf galaxies are thus promising sites for testing CDM in the strongly nonlinear regime and learning about the identity of the dark matter and the effects of baryonic processes.
The most direct way to study the central density structure of a gaspoor galaxy is by fitting an equilibrium stellar dynamical model to a large sample of stars that have high resolution spectroscopy and good photometry. In recent years, data of the required quality have been obtained for a number of nearby dwarf spheroidal galaxies (dSphs) around the Milky Way (Simon & Geha 2007; Walker et al. 2009) and M31 (Tollerud et al. 2012). Simple dynamical analyses based on spherical symmetry and the Jeans equations suffer from degeneracies which preclude an unambiguous determination of the dark matter potential (Walker 2013; Strigari 2013). Thus, data for several dSphs have been shown to be equally consistent with flat central profiles (cores) (Gilmore et al. 2007) or with NFW cusps (Strigari et al. 2010; Jardel & Gebhardt 2013). In some cases, one can hope to break degeneracies by considering higher moments of the lineofsight velocity distribution (Richardson & Fairbairn 2014).
Sculptor, a dSph of stellar mass M located kpc from the Galactic Centre (Lianou & Cole 2013), is a particularly interesting case. Modelling using a variety of techniques but treating the available stellar data as sampled from a single stellar population and assuming spherical symmetry has shown that the kinematic data are consistent with an NFW halo potential, but also allow a core (Strigari et al. 2010; Breddels et al. 2013; Richardson & Fairbairn 2014). These studies suggest a dark matter halo mass of M for Sculptor though with rather large uncertainties (Łokas 2009; Strigari et al. 2010; Breddels et al. 2013).
The data for Sculptor are of sufficient quality that two distinct stellar populations of differing metallicity can be identified: a centrally concentrated metalrich (MR) population and a more extended metalpoor (MP) population (Battaglia et al. 2008, B08). The presence of two populations makes it possible to carry out more refined dynamical analyses. Thus, applying the Jeans equations to each population separately, B08 showed that their data could be fit by a model in which the orbital distribution of each population is isotropic near the centre and becomes radially biased in the outer regions. They found a best fit for a model potential with a core, but also found the data to be consistent with an NFW potential. Using MichieKing models for the stellar distribution function, Amorisco & Evans (2012) also found that while an NFW model provides an acceptable fit to the data, models with a core seem to be preferred. On the other hand, applying the projected virial theorem, Agnello & Evans (2012) concluded that it is not possible to fit both the MR and the MP populations with a single NFW model.
An independent dynamical analysis of Sculptor using a larger sample of stars was carried out by Walker & Penarrubia (2011, WP11). Rather than simply separating the observed stars into two populations according to their estimated metallicity, they devised a statistical method which fits the full dataset simultaneously with two constant velocity dispersion, Plummerprofile populations of differing metallicity, together with a contaminating Galactic component. They then inserted the halflight radius and velocity dispersion estimated for each population into the mass estimator proposed by Walker et al. (2009). This allowed them to infer the mass contained within each halflight radius and thus the slope of the density profile between the two halflight radii. They concluded that the slope is flatter than predicted for an NFW profile at the 99% c.l.
In this study we carry out a new analysis of Sculptor in an attempt to clarify the conflicting claims in the literature. The specific statistical question we ask is whether the kinematic and photometric data for this galaxy exclude potentials of the type predicted by CDM. We reexamine both the B08 and WP11 datasets from a different theoretical perspective and discuss how they compare. There are both similarities and differences between our analysis and those that have been undertaken previously. Like B08 and Amorisco & Evans (2012), but unlike Agnello & Evans (2012) and WP11, we exploit the full information contained within the lineofsight velocity dispersion and photometry profiles. Like Amorisco & Evans (2012), but unlike B08, we build models based on distribution functions. Our analysis differs from that of Amorisco & Evans (2012) primarily in that we use a more flexible form for the distribution function which allows a wider range of energy distributions and velocity anisotropies for the stars.
We conclude, in agreement with B08 and Amorisco & Evans (2012), that NFW potentials are not excluded by the B08 data. For our more general models the constraints used by Agnello & Evans (2012) are also no longer sufficient to exclude NFW potentials. Indeed, the implied peak circular velocity of the Sculptor dark matter halo and its concentration are consistent with the values predicted from CDM simulations. While an NFW potential gives an acceptable fit to the Sculptor data, our analysis cannot exclude potentials with a core; indeed, we find below that a potential of the form proposed by Burkert (1995) provides a slightly (though not significantly) better fit to the B08 data than an NFW profile. We find similar conclusions when we apply our analysis directly to the WP11 data. Thus, our discrepant conclusions reflect differences in analysis methods rather than in observational datasets.
This paper is organized as follows. In Section 2 we introduce our model for the stellar distribution function. In Section 3 we briefly discuss our methodology for fitting the theoretical model to the data. In Section 4 we present our results and, in Section 5, we compare them to previous studies, highlighting discrepancies where they exist.
2. Models
In this section we introduce the dynamical models we use to interpret the observed stellar populations in Sculptor. We assume each population to be spherically symmetric and to be in dynamical equilibrium within a static and spherically symmetric potential well. These are strong assumptions which should be treated as approximations. The observed stellar distribution is clearly noncircular on the sky, and Sculptor orbits within the potential of the Milky Way, so the effective potential seen by its stars is timevarying. Some aspects of the effects of flattened potentials on the dynamical analysis of dSph data are considered by Laporte et al. (2013).
2.1. Dark matter
For the total mass density profile of the system we adopt two standard models. First, an NFW model,
(1) 
with corresponding gravitational potential
(2) 
where we define and . This simple model is determined by just two scale parameters, the characteristic radius, , and the characteristic density, . Note that we define the potential to be zero at the centre of the system and to be at infinity. NFW models are often parametrized in terms of the maximum circular velocity, , and the radius, , at which this is attained. These are related to and through:
(3) 
Our second model is the cored profile proposed by Burkert (1995),
(4) 
where here . This profile also has two parameters, the central density and the scale radius , and we define . The corresponding gravitational potential is
(5) 
As for the NFW case, we define the potential to be zero at the centre of the system. With these definitions, as we have . For the Burkert model, the maximum of the circular velocity curve and the radius at which it is attained are related to the other parameters through
(6) 
2.2. Stellar distribution function
We define the specific energy and specific angular momentum of a star as and , respectively, where is the modulus of the velocity vector and is the angle between this vector and the star’s position vector relative to system centre. Given a static and spherically symmetric gravitational potential well, any positive definite function corresponds to the phasespace distribution function of some stable, dynamically mixed and spherically symmetric equilibrium for a stellar population. In this paper we will consider only models in which the dependence on and is separable,
(7) 
with both and positive definite and given by simple parametric forms. It would be easy to build more general, nonseparable models as a superposition of several individually separable components, but we will not pursue this further here.
The stellar density profile and the radial and tangential stellar velocity dispersion profiles of such models are given by
(8)  
(9)  
(10) 
where . Note that with the definition we are using here, the total velocity dispersion at radius is
(11) 
Eqns 8, 9, and 10 can be combined to give the projected stellar density profile and stellar lineofsight velocity dispersion profile at a fixed projected distance :
(12)  
(13) 
where .
A particularly interesting and simple case occurs when the angular momentum dependence is taken to be a power law,
(14) 
where is a constant. For this assumption, the integrals over and separate in Eqns 8, 9 and 10, and the ratio of the two velocity dispersions is independent both of and of . The lower limit on is required for the integrals to converge for small . For this choice of the orbital anisotropy of the stellar population model, usually parametrized as
(15) 
is independent of radius and depends on alone, . For an isotropic velocity distribution, . For nearradial orbits is close to unity and approaches its lower limit of , while for nearcircular orbits is very large and positive while is very large and negative.
In this paper we will investigate models where the orbital anisotropy varies with radius and we therefore need a more general form for . We consider the function,
(16) 
which interpolates between a power law of index at and a power law of index at . The parameter controls the rapidity of the transition between the two regimes at the characteristic scale, , which corresponds to a radius of order , where for NFW and for the Burkert case. In addition, is required to be positive for and to be negative in the opposite case.
For simplicity when comparing with the Sculptor data, we in this paper prefer to use a function with fewer free parameters and to assume that the velocity distribution is isotropic near the centre, as seems plausible on general theoretical grounds. We therefore set and , resulting in the simpler expression
(17) 
The upper and lower cases here correspond to radially and tangentially biased orbits for large angular momenta, respectively. Both produce isotropy for small angular momenta and so also at small radii. This simplified model retains only two parameters, which sets the extent of the inner isotropic region and which determines the velocity anisotropy for large angular momenta.
For the energy distribution, , we have found the following form to be sufficiently general for our purposes:
(18) 
where the restriction is required because orbits with are unbound. The normalisation, , in this expression sets the amplitude of the stellar density profile, while the exponent determines the behaviour at small energies, hence as . Comparison with the simple scalefree distribution functions explored by White (1981) shows that at sufficiently small radii (where , and ) Eqns 2, 16 and 18 imply a powerlaw stellar density profile, , where
(19) 
When , the density in the innermost regions is dominated by stars on orbits which are confined to those regions, while in the contrary case it is dominated by stars on orbits which extend well beyond them. Our model for thus produces a central cusp in the stellar density profile when Eqn. 19 gives .
At somewhat larger energies, (hence at radii larger than , where ) the density profile steepens to a new slope, , which is given by Eqn 19 with replaced by where we assume . The rapidity of the transition around is controlled by the parameter . The final factor in Eqn 18 allows for truncation of the stellar density at a radius, , defined by , which is directly analogous to the “tidal radius” in the classic King models for globular clusters (King 1966). The shape of this cutoff in the profile can be adjusted using the final parameter, . A special case arises when . Then and the density profile at large radii becomes a power law of slope (or for the simpler case of Eqn 17).
As a final remark, we note that when , the constraint forces stellar velocities and hence stellar angular momenta to be small as . The anisotropy in this region is thus determined by the form of for small rather than large , with the result that rather than (i.e. the distribution becomes isotropic again as if the simpler parametrisation of Eqn 17 is used). This complication does not arise when , in which case is indeed equal to at large radii ( for the simpler model of Eqn 17).
3. Data analysis
In this section we briefly detail the Bayesian analysis methods that we will apply to the datasets described in the following sections.
For a single stellar population, the model of Eqns 7, 17 and 18 is a function of nine parameters, . Including the two parameters that describe the NFW potential, Eqn 2, a fit to an individual stellar population has 11 free parameters, whereas a joint fit to both populations (each of which has independent distribution function parameters), has a total of 20 free parameters.
We use our theoretical model to fit to the binned velocity dispersion data. We define the quantities,
(20)  
(21) 
where the subscript, , denotes a specific population, either or . Additionally is the number of data points in the photometric profile of a population, and the number of data points in the velocity dispersion profile of a population, both of which are measured at projected distance, . The associated measurement uncertainties are and . Using these quantities we define a full likelihood function of the form,
(22) 
where
(23) 
The constant in Eqn 22 depends on the photometric and velocity dispersion measurement uncertainties but does not depend on the distribution function model.
We employ an MCMC algorithm which is an adapted version of the CosmoMC code (Lewis & Bridle 2002). From this algorithm we can extract two important quantities: (i) the maximum value of the likelihood, , which corresponds to a minimum value of and to the “best fit” set of parameters from a given chain, and (ii) the posterior probability distribution for each model parameter. For scans of a large and complex theoretical parameter space, MCMC algorithms are not necessarily effective at finding the true value of , so it is important to determine if the estimated value of corresponds to a set of theoretical parameters that provide a statistically good fit to the data. We thus run several chains from different starting points in the theoretical parameter space to ensure that the chains find similar values of and are thus not burning in at local maxima in the likelihood.
The fact that we are marginalizing over up to twenty parameters also means that we must test that the posterior probability distributions for the model parameters have appropriately converged. We test for convergence of the posterior probability distributions in a standard manner by estimating the variance of a parameter as a weighted sum of the withinchain and betweenchain variance (Gelman & Rubin 1992).



4. Results
B08 obtained spectroscopy for 470 stars in Sculptor from which they measured lineofsight velocities and metallicites derived from the Calcium triplet lines. They identified two distinct populations with different metallicity, spatial distribution and kinematics: a metal rich (MR) population, defined to have , and a metal poor (MP) population, defined to have . This clean separation and the large radial coverage of the two populations make this an attractive sample to analyze using our model distribution function described in Section 2. We fit our distribution function model to the photometric and velocity dispersion profiles reported by B08 for each population by performing a likelihood analysis using the MCMC technique, as described in Section 3. Note that the photometric profiles are obtained from a different and much larger population of stars than the velocity dispersion profiles.
Since the stellar distribution of Sculptor is elongated on the sky, B08 give the surface brightness profile of each population as a function of an “elliptical radius” which corresponds to the projected semimajor axis determined from the photometry (Tolstoy et al. 2004). To account for this in the context of our assumption of spherical symmetry, we take as the radial coordinate the geometric mean of the major and minor axes, which we expect to correspond best to the count profile for circular annuli. We also perform the same scaling on the B08 kinematic data. The ellipticity of Sculptor is (Irwin & Hatzidimitriou 1995).
We fit the full stellar density and velocity dispersion profiles of the two metallicity populations to the 20parameter model defined in Section 2, which here assumes an NFW potential, and in which the velocity anisotropy can vary with radius but is assumed isotropic at the centre. From the MCMC chains we obtain both the maximum likelihood value and the posterior probability distribution for each model parameter. The surface density and velocity dispersion profiles for a model that has nearmaximal likelihood are shown in Figure 1. The count profiles of both stellar populations exhibit welldefined cores. The MP population in this particular model is isotropic everywhere, while the MR population is isotropic in the centre, has a sharp transition at a scale radius of kpc to over the range kpc, and then transitions smoothly back to at larger radii in order to satisfy the boundary conditions at . The parameters of this model are listed in Table 1.
Figure 1 shows that the data for the two metallicty populations in Sculptor are very well fit by our model, and this impression is confirmed by the values of for the fits: for the MP photometry (left; 23 data points), for the MR photometry (middle; 15 data points), for the MP kinematics (upper right; 8 data points), and for the MR kinematics (lower right; 5 data points). Our analysis therefore demonstrates that the data are consistent with both populations residing in a single NFW potential.
Population  

MR  2.0  5.3  2.5  0.16  0.45  1.5  9.0  6.9  21  1.5  
MP  2.4  7.9  1.1  0.17  0.60  3.0  0  8.2  – 
Consistency with a potential of the NFW form does not, however, guarantee that the data are consistent with the predictions of the CDM model. For this to be the case, the parameters of the NFW density profile must lie within the theoretically predicted range. The comparison is most easily carried out in space, where the maximum circular velocity of the dark matter halo, , and the radius, , at which it is attained are defined in Eqn. 3 and are readily measured for subhalos in high resolution CDM Nbody simulations of galactic halos.
The region in the plane in which 90% of the subhalos in the “Aquarius” CDM simulations of Springel et al. (2008) lie is shown by the thin lines in Figure 2. The thick lines show the 68% and 90% contours of twodimensional joint posterior probability distributions of and derived from our fits to the B08 data. These contours overlap well with the theoretically predicted region, demonstrating that the kinematics of the populations in Sculptor are fully consistent with expectations in a CDM universe. The range of allowed by our fits, km/s, is significantly wider than the range estimated in some previous analyses (BoylanKolchin et al. 2011).

To compare the quality of our NFW fits to that obtained for a cored potential, Figure 3 presents results from a joint analysis of the two populations using a Burkert profile. In this case, we also have 20 parameters, but we replace by . The bestfit values in this figure are almost identical to those found in the NFW case, but the constraints are considerably tighter, reflecting the much more sharply defined characteristic scale of the cored potential. The total for the best fit is in the Burkert case which is slightly but not significantly smaller than the value of 40.2 which we found for the bestfitting NFW potential.

Given the similar quality of the fits for the two profiles, we consider the more general question of whether we can expect to distinguish them with kinematic and photometric datasets similar to those of B08. We start from distribution function parameters resembling those of our best Burkert fit to B08, and generate mock photometric and kinematic datasets with similar size and uncertainties to B08. We consider several values for the Burkert scale radius and scale density, chosen to give dispersion and count profiles similar to those observed. The largest value for the scale radius we consider is kpc. We fit these mock data to Burkert and NFW models, comparing the best fit values in the two cases. Even for kpc, we find NFW models with very similar to the best fitting Burkert model. The latter always has parameters very close to the input values, showing our procedure to be approximately unbiased. This test implies that the B08 data samples are not large enough to be able to distinguish between Burkert and NFW potentials on the basis of count and velocity dispersion profiles. We will return to the issue of distinguishing between these two models below when discussing the WP11 data.
5. Comparison to previous results
Several previous studies have constrained the potential in Sculptor by splitting its stars into high and low metallicity populations and requiring each separately to be in equilibrium. In this section we compare our results with these earlier analyses.
5.1. Battaglia et al. (2008)
In the first dynamical analysis to separate the two populations in Sculptor, B08 applied the Jeans equation to each population individually, fitting the star count profiles to standard forms (Plummer for the metalrich, Sersic for the metalpoor) and then predicting the velocity dispersion profiles for various assumed potentials and anisotropy profiles. They found that their data were consistent with both NFW and core potentials, but required radially biased orbits for both populations.
Our results largely corroborate the conclusions reached by B08. Our use of a flexible stellar distribution function removes the need to assume standard forms for the count and anisotropy profiles and ensures that the resulting model is physically realisable. As a result, our bestfit model has a smaller than the model of B08. Like B08, we find that radiallybiased orbits are required at large radius for the metalrich (although not for the metalpoor) population.
5.2. Amorisco & Evans (2012)
Of the previous studies, that of Amorisco & Evans (2012) is most similar to our own. They also assumed separable distribution functions for the two populations and fit predicted counts and velocity dispersion profiles to the data of B08. They considered both NFW potentials and pseudoisothermal potentials with a core. Although the form they assumed for their distribution functions is considerably less flexible than our own, the resulting best fit for an NFW potential is similar to ours, shown in Fig. 1, and has a value which is clearly insufficient to exclude the model. The best fit for the core case also looks similar (compare their figures 9 and 10). Nevertheless, its value is sufficiently smaller that a likelihood ratio test clearly prefers it over an NFW potential.
As Amorisco & Evans (2012) note, the preference for a core potential over a cuspy one is driven in their analysis by its lower prediction for the innermost points of the count profiles and, to a lesser extent, by a somewhat larger predicted difference in velocity dispersion between the two populations. With our more flexible distribution function model, the count discrepancy at small radius disappears for the NFW potential and the difference in velocity dispersions between the MP and MR populations is slightly enhanced (see Fig. 1), leading to a fit of very similar quality to that found by Amorisco & Evans (2012) for their core potential. A final difference with Amorisco & Evans (2012) is that their NFW fit required a halo concentration which is lower than expected in CDM. With our distribution function model, this problem has disappeared.
5.3. Agnello & Evans (2012)
Agnello & Evans (2012) applied the projected virial theorem separately to the two populations identified by B08 assuming that they reside in an NFW potential and have Plummerlaw surface brightness profiles. With these assumptions, the observational data for each population imply a relation between and for its halo. They then show that the regions of the (, ) plane allowed at 2 by the MR and MP data do not overlap. They therefore conclude that no single NFW potential can accommodate both populations.
B08 noted that the metalpoor population in Sculptor is poorly fit by a Plummer model. By applying our procedures to the MR and MP data separately, we can perform an analysis analogous to that of Agnello & Evans (2012). In Figure 4 we show the twodimensional joint posterior probability distributions of , . We indeed reproduce an offset similar to that seen by Agnello & Evans (2012). However, the greater freedom afforded by our relaxation of the Plummerlaw assumption, allowing instead any profile consistent with the observed counts, widens the confidence regions so that they are no longer exclusive. The two distributions are marginally consistent with each other for km/s and kpc which, not suprisingly, is the region also picked out by our single potential model. As before, these parameters are consistent with the CDM predictions of Springel et al. (2008).

5.4. Walker & Penarrubia (2011)
WP11 analyzed a sample of 1497 stars with spectroscopy from which they derive lineofsight velocities and an Mg index which they take as a proxy for metallicity. In practice, we are able to work with the subsample of 1307 stars for which the public online dataset lists a membership probability value. WP11 use MCMC techniques to map the parameter distributions for a 3component model of these data. The metalrich and metalpoor populations are each represented by a circularly symmetric Plummer profile, with Gaussian velocity and metallicity distributions independent of radius (note that these assumptions are not consistent with the B08 data), while the contaminating Galactic foreground is taken to be spatially uniform with broader distributions of velocity and metallicity. They insert the halflight radius and velocity dispersion estimated by this analysis for each population into the mass estimator proposed by Walker et al. (2009):
(24) 
which gives the mass, , inside a sphere with radius equal to
the projected halflight radius, , in terms of the measured
velocity dispersion, , and
This conclusion is incompatible with our conclusion derived above, based on the B08 data. In Fig. 5 we show the results of WP11 in the (, ) plane, together with lines corresponding to , with and . Clearly, these results agree much better with the dotted line representing a core than with the dashed line representing an NFW cusp. Our distribution function based MCMC analysis allows us to reconstruct and for all models consistent with the B08 data and residing in an NFW potential. Solid red and blue contours in Fig. 5 give the 68% and 90% confidence regions for the metalrich and metalpoor populations respectively. As expected, the centre points of these contours define a slightly shallower slope than the dashed line since only in the innermost regions of an NFW profile. The halflight radii found for the MR and MP populations in the two analyses agree well but there is an offset in the preferred values, although the contours do overlap. Hence, fitting our models to the B08 data has resulted in lower velocity dispersions for the MP and higher velocity dispersions for the MR population than estimated by WP11 from their own data.

To determine whether the WP11 data favor a cored or cusped halo when analyzed using our procedures, we construct binned photometric and velocity dispersion profiles directly from the WP11 data. Figure 6 shows the metallicity as a function of both radius and velocity for the 1160 stars in the tables published by WP11 that have high quality data and are assigned a membership probability %. The W distribution of the member stars is unimodel, with no obvious indication of two distinct populations. To separate the stars by metallicity into two populations with maximally distinct spatial distributions, we split them at a given value of in the left panel of Fig. 6, and then perform a KS test to determine the significance of the difference in radial distribution between the two populations. The value of which minimizes the KS pvalue then defines the MR and MP populations with the most significantly distinct spatial distributions.


Following this procedure, we find a minimum pvalue of at . This cut gives MP members with , and MR members with . The left panel of Fig. 6 shows that, with this cut, the fraction of MR stars at large radius is noticeably smaller than for the MP stars. Including all 1307 WP11 stars with a membership probability, regardless of its value, increases the pvalue by nearly two orders of magnitude, but still gives an optimal separation at about the same W’ value, and with a similar ratio MR to MP stars. For comparison, in their statistical separation, WP11 found that % of the Sculptor member stars belong to the underlying ”true” MR population, and the remainder to the MP population.
The velocity dispersion and the photometry profiles of the two populations are shown in Figure 8. In comparison to the B08 data in Figure 1, the WP11 data cover a narrower range of radii but have higher signal to noise for the velocity dispersion measurements. For all four profiles, we have corrected the raw numbers using the radial selection function as suggested in WP11. We find the halflight radii for the MR and MP populations to be 0.18 and 0.22 kpc, respectively. The MR halflight radius is in good agreement with that derived by WP11, while the MP halflight radius is significantly smaller than the value of 0.30 kpc derived by WP11. The velocity dispersion profile of the MR component derived from the WP11 data does not show a steep decline at large radii of the kind seen in the B08 data.

An MCMC analysis of the WP11 results in allowed regions in the parameter space shown in Figure 7 for both Burkert and NFW profiles. For the NFW profile, the derived values of are in good agreement with the values derived from the B08 data, while for the Burkert profile the central is larger, peaking at km/s. The corresponding best fitting photometry and velocity dispersion profiles are shown in Figure 8. The for the best fitting profiles for both joint and individual fits are given in Table 2. Note that in both cases we are fitting 34 data points to a 20 parameter model. so the values we find indicate fully acceptable fits and show a slight preference for NFW over Burkert.



NFW  Burkert  

photometry  dispersion  photometry  dispersion  
MR  5.5  3.9  6.0  2.8 
MP  5.0  1.2  5.5  3.2 
We thus conclude that when analyzed with the methods of this paper the WP11 data do not favour a core over a cusp. This result is in disagreement with that of WP11, who find NFW profiles to be ruled out with % confidence. At this point we are unable to determine why our conclusions differ from those of WP11. Our dynamical analysis is based on selfconsistent distribution function models with considerable flexibility which can be exploited to obtain good simultaneous fits to the photometric and kinematic data. In contrast, WP11 use simple, constant velocity dispersion Plummerlaw models which are inconsistent in detail with the Sculptor data. However, since these are used only to estimate the halflight radii and total velocity dispersions of the two populations, it is unclear whether this inconsistency can significantly bias their results. A second difference is that our likelihood analysis explicitly separates Sculptor stars into two populations based on the directly measured quantity W’ and rejects ab initio the small number of stars which are not high probability members, whereas WP11 treat W’ as an indicator of the relative probability of belonging to one or the other of two assumed underlying populations, and also treat background rejection probabilistically. However, at the present time it is unclear whether these differences can account for the difference in our conclusions. Fig. 8 does appear to demonstrate that our model can reproduce the WP11 data very well in an NFW potential.
6. Conclusions
In this paper we have presented a new distribution function based framework to study stellar populations in equilibrium within a spherical dark matter potential well. We use it to study the two metallicity populations in the Sculptor dwarf spheroidal galaxy, in particular, to explore the controversial question of whether their properties exclude a cuspy profile of the kind expected in the CDM cosmology.
The family of distribution functions we consider gives substantially more freedom than the models assumed in previous studies and, as a result, leads to a weakening of the constraints implied by the observations. Although in the absence of any prior on the shape of the inner potential, we concur with previous studies that the Sculptor data prefer a shallower profile than NFW, we find this preference to be far too weak to exclude the cosmological prediction. Indeed, in a sense, we are able to find equilibrium models which are a good fit to the datasets of both B08 and WP11 within an NFW potential with parameters that are fully consistent with CDM.
Since the inner structure of dwarf galaxies appears at present as one of the few significant challenges to the standard cosmological paradigm it is unsurprising that considerable attention has been focused on measuring this structure precisely. Unfortunately, the problem is underconstrained by currently available data, given the considerable freedom inherent in the equations of stellar dynamics. The analysis in this paper, while comparatively general, still makes at least two major assumptions which are known to be incorrect: dwarf spheroidal galaxies are clearly not spherically symmetric and their orbits within the Milky Way’s potential ensure that most cannot be static systems in equilibrium. Further theoretical progress will require these shortcomings to be addressed (See Zhu et al (2016) for a recent study of Sculptor which relaxes the assumption that the stelllar distribution is spherical, while continuing to assume a spherical potential.) Further observational progress may be achieved by reducing the statistical and measurement uncertainties and, in the more distant future, by increasing the phasespace coverage through measurement of internal proper motions.
Acknowledgements
We would like to thank Chervin Laporte for his critical comments on an earlier version of this work which prompted us to develop the distribution function methods used in this paper. We also thank Matt Walker for critical and constructive comments on the paper. Much of the theoretical modelling framework for this paper was developed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY1066293. LES acknowledges support from National Science Foundation grant PHY1522717. This material is based upon work supported by the National Science Foundation under Grant No. CNS0723054. CSF acknowledges ERC Advanced Investigator Grant COSMIWAY and SDMW ERC Advanced Investigator Grant GALFORMOD. This work was supported in part by an STFC rolling grant to the ICC. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National Einfrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National Einfrastructure.
Footnotes
 This estimator is constructed to be only weakly sensitive to the details of the density and velocity anisotropy profiles (see also Wolf et al. 2010)
References
 Agnello A., Evans N. W., 2012, ApJ, 754, L39
 Amorisco N. C., Evans N. W., 2012, MNRAS, 419, 184
 Battaglia G., Helmi A., Tolstoy E., Irwin M., Hill V., Jablonka P., 2008, ApJ, 681, L13
 BoylanKolchin M., Bullock J. S., Kaplinghat M., 2011, Mon.Not.Roy.Astron.Soc., 415, L40
 Breddels M. A., Helmi A., van den Bosch R. C. E., van de Ven G., Battaglia G., 2013, MNRAS, 433, 3173
 Burkert, A., 1995, ApJ, 447, L25
 Frenk C. S., White S. D. M., 2012, Annalen der Physik, 524, 507
 Gelman S., Rubin D., 1992, Statistical Science, 7, 457?472
 Gilmore G., Wilkinson M. I., Wyse R. F. G., Kleyna J. T., Koch A., Evans N. W., Grebel E. K., 2007, ApJ, 663, 948
 Irwin M., Hatzidimitriou D., 1995, Mon.Not.Roy.Astron.Soc., 277, 1354
 Jardel J. R., Gebhardt K., 2013, Astrophys.J., 775, L30
 King I. R., 1966, AJ, 71, 64
 Laporte C. F. P., Walker M. G., Peñarrubia J., 2013, MNRAS, 433, L54
 Lewis A., Bridle S., 2002, Phys.Rev., D66, 103511
 Lianou S., Cole A. A., 2013, A&A, 549, A47
 Łokas E. L., 2009, MNRAS, 394, L102
 Lovell M. R., Eke V., Frenk C. S., Gao L., Jenkins A., Theuns T., Wang J., White S. D. M., Boyarsky A., Ruchayskiy O., 2012, MNRAS, 420, 2318
 Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS, 283, L72
 Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
 Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
 Newman A. B., Treu T., Ellis R. S., Sand D. J., 2013, ApJ, 765, 25
 Peñarrubia J., Pontzen A., Walker M. G., Koposov S. E., 2012, ApJ, 759, L42
 Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS, 430, 105
 Pontzen A., Governato F., 2014, Nature, 506, 171
 Richardson T., Fairbairn M., 2014, MNRAS, 441, 1584
 Shao S., Gao L., Theuns T., Frenk C. S., 2013, MNRAS, 430, 2346
 Simon J. D., Geha M., 2007, Astrophys.J., 670, 313
 Springel V., Wang J., Vogelsberger M., Ludlow A., Jenkins A., et al., 2008, Mon.Not.Roy.Astron.Soc., 391, 1685
 Strigari L. E., 2013, Phys.Rept., 531, 1
 Strigari L. E., Frenk C. S., White S. D., 2010, Mon.Not.Roy.Astron.Soc., 408, 2364
 Tollerud E. J., Beaton R. L., Geha M. C., Bullock J. S., Guhathakurta P., et al., 2012, Astrophys.J., 752, 45
 Tolstoy E., Irwin M. J., Helmi A., Battaglia G., Jablonka P., Hill V., Venn K. A., Shetrone M. D., Letarte B., Cole A. A., Primas F., Francois P., Arimoto N., Sadakane K., Kaufer A., Szeifert T., Abel T., 2004, ApJ, 617, L119
 Walker M., 2013, Planets, Stars and Stellar Systems. Volume 5: Galactic Structure and Stellar Populations, Oswalt, T. D. and Gilmore, G., 1039, p. 1039
 Walker M. G., Mateo M., Olszewski E. W., 2009, AJ, 137, 3100
 Walker M. G., Mateo M., Olszewski E. W., Penarrubia J., Evans N. W., et al., 2009, Astrophys.J., 704, 1274
 Walker M. G., Penarrubia J., 2011, Astrophys.J., 742, 20
 White S. D. M., 1981, MNRAS, 195, 1037
 Wolf J., Martinez G. D., Bullock J. S., Kaplinghat M., Geha M., et al., 2010, Mon.Not.Roy.Astron.Soc., 406, 1220
 Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20