Dynamical models for the Sculptor dwarf spheroidal in a \LambdaCDM universe

Dynamical models for the Sculptor dwarf spheroidal in a CDM universe

Dynamical models for the Sculptor dwarf spheroidal in a CDM universe

Louis E. Strigari Mitchell Institute for Fundamental Physics and Astronomy, Dep. of Physics and Astronomy, Texas A & M Univ., College Station, TX 77843, USA Carlos S. Frenk Institute for Computational Cosmology, Dep. of Physics, Univ. of Durham, South Road, Durham DH1 3LE, UK Simon D. M. White Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85740 Garching bei München, Germany


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 Navarro-Frenk-White (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 self-consistent 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 large-scale 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 non-linear (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 gas-poor 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 line-of-sight 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 metal-rich (MR) population and a more extended metal-poor (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 Michie-King 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, Plummer-profile populations of differing metallicity, together with a contaminating Galactic component. They then inserted the half-light 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 half-light radius and thus the slope of the density profile between the two half-light 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 re-examine 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 line-of-sight 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 non-circular on the sky, and Sculptor orbits within the potential of the Milky Way, so the effective potential seen by its stars is time-varying. 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,


with corresponding gravitational potential


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:


Our second model is the cored profile proposed by Burkert (1995),


where here . This profile also has two parameters, the central density and the scale radius , and we define . The corresponding gravitational potential is


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


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 phase-space 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,


with both and positive definite and given by simple parametric forms. It would be easy to build more general, non-separable 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


where . Note that with the definition we are using here, the total velocity dispersion at radius is


Eqns 89, and 10 can be combined to give the projected stellar density profile and stellar line-of-sight velocity dispersion profile at a fixed projected distance :


where .

A particularly interesting and simple case occurs when the angular momentum dependence is taken to be a power law,


where is a constant. For this assumption, the integrals over and separate in Eqns 89 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


is independent of radius and depends on alone, . For an isotropic velocity distribution, . For near-radial orbits is close to unity and approaches its lower limit of , while for near-circular 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,


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


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:


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 scale-free distribution functions explored by White (1981) shows that at sufficiently small radii (where , and ) Eqns 2,  16 and 18 imply a power-law stellar density profile, , where


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 cut-off 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 717 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,


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,




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 within-chain and between-chain variance (Gelman & Rubin 1992).

Figure 1.— Best-fit model of Eqns. 17 and 18 to the Battaglia et al. (2008) data for the metal-poor and metal-rich populations in Sculptor, assuming an NFW potential. The left and middle panels show the surface brightness profiles of the MP  and MR populations respectively; the right panel shows the velocity dispersion profile of the MP  population (blue line at the top) and the MR population (red line at the bottom).

4. Results

B08 obtained spectroscopy for 470 stars in Sculptor from which they measured line-of-sight 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 semi-major 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 20-parameter 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 near-maximal likelihood are shown in Figure 1. The count profiles of both stellar populations exhibit well-defined 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.

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
Table 1An example distribution function model that provides a good fit to the Sculptor two-population data. and are in units of and, for the MR population, is in units of ; is in km/s and in kpc.

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 N-body 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 two-dimensional 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 (Boylan-Kolchin et al. 2011).

Figure 2.— The 68% and 90% c.l. regions in the plane for fits of the model of Eqns. 17 and 18 to the two metallicity populations in Sculptor (thick lines). A cross indicates the model of Table 1. Thin lines delineate the region which contains 90% of the subhalos from the Aquarius CDM simulations of galactic halos (Springel et al. 2008).

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 best-fit 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 best-fitting NFW potential.

Figure 3.— 68% and 90% confidence contours in the plane for independent fits to the two metallicity populations in Sculptor assuming a Burkert profile. The red contours are obtained by fitting the 11 parameter model of Eqns. 17 and 18 to the B08 MR data, while the blue contours are obtained from fitting the same 11 parameter model to the B08 MP data. These contour sets can be compared with Figure 4 for the NFW profile. The black contours are the results of a joint fit to the MR and MP data as shown for the NFW profile in Figure 2.

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 metal-rich, Sersic for the metal-poor) 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 best-fit model has a smaller than the model of B08. Like B08, we find that radially-biased orbits are required at large radius for the metal-rich (although not for the metal-poor) 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 pseudo-isothermal 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 Plummer-law 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 metal-poor 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 two-dimensional 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 Plummer-law 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).

Figure 4.— The 68% and 90% c.l. contours in the plane for independent fits to the two metallicity populations in Sculptor. The red contours are obtained from fitting the 11 parameter model of Eqns. 17 and 18 to the B08 MR data, while the blue contours are obtained from fitting the 11 parameter model to the B08 MP data. Crosses in each case indicate the best-fit model.

5.4. Walker & Penarrubia (2011)

WP11 analyzed a sample of 1497 stars with spectroscopy from which they derive line-of-sight 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 3-component model of these data. The metal-rich and metal-poor 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 half-light radius and velocity dispersion estimated by this analysis for each population into the mass estimator proposed by Walker et al. (2009):


which gives the mass, , inside a sphere with radius equal to the projected half-light radius, , in terms of the measured velocity dispersion, , and 1. The derived increase in estimated mass between the two values of appears too large to be consistent with an NFW profile and is close to that expected for a constant density core. WP11 conclude that NFW is excluded at the 99% c.l.

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 metal-rich and metal-poor 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 half-light 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.

Figure 5.— Constraints from WP11 and from the models of this paper in the () plane. The two straight lines and the contours traced by open circles are taken directly from Figure 10 of WP11 and indicate with and the 50% c.l. regions given by their MCMC analysis for the parameters of the two underlying populations. For comparison, the solid contours show 68% and 90% c.l. regions from our own MCMC chains constrained by the B08 data and assuming both populations to be in equilibrium within a single NFW potential.

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 p-value then defines the MR and MP populations with the most significantly distinct spatial distributions.

Figure 6.— Reduced Mg index, , versus radius (left) and velocity (right) for stars considered to be probable Sculptor members by WP11. The histogram plotted vertically (center) shows the distribution of for the sample as a whole. The red dashed line defines the cut which maximally separates the radial distributions of metal-rich and metal-poor populations.

Following this procedure, we find a minimum p-value 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 p-value 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 half-light radii for the MR and MP populations to be 0.18 and 0.22 kpc, respectively. The MR half-light radius is in good agreement with that derived by WP11, while the MP half-light 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.

Figure 7.— The 68% and 90% c.l. contours in the plane for joint fits to the two metallicity populations in Sculptor using the WP11 data split according to the cut in Fig.6. Solid contours are for an NFW and dashed are for a Burkert potential.

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.

Figure 8.— Results of joint fits to the  MP and MR  photometry and kinematics, using the WP binned data. In each panel, the solid line assumes an NFW, and the dashed line a Burkert potential.
NFW Burkert
photometry dispersion photometry dispersion
MR 5.5 3.9 6.0 2.8
MP 5.0 1.2 5.5 3.2
Table 2Total values for the best-fit to the photometry and kinematic profiles using the WP11 data.

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 self-consistent 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 Plummer-law models which are inconsistent in detail with the Sculptor data. However, since these are used only to estimate the half-light 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 phase-space coverage through measurement of internal proper motions.


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 PHY-1066293. LES acknowledges support from National Science Foundation grant PHY-1522717. This material is based upon work supported by the National Science Foundation under Grant No. CNS-0723054. 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 E-infrastructure 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 E-infrastructure.


  1. 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)


  1. Agnello A., Evans N. W., 2012, ApJ, 754, L39
  2. Amorisco N. C., Evans N. W., 2012, MNRAS, 419, 184
  3. Battaglia G., Helmi A., Tolstoy E., Irwin M., Hill V., Jablonka P., 2008, ApJ, 681, L13
  4. Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, Mon.Not.Roy.Astron.Soc., 415, L40
  5. Breddels M. A., Helmi A., van den Bosch R. C. E., van de Ven G., Battaglia G., 2013, MNRAS, 433, 3173
  6. Burkert, A., 1995, ApJ, 447, L25
  7. Frenk C. S., White S. D. M., 2012, Annalen der Physik, 524, 507
  8. Gelman S., Rubin D., 1992, Statistical Science, 7, 457?472
  9. Gilmore G., Wilkinson M. I., Wyse R. F. G., Kleyna J. T., Koch A., Evans N. W., Grebel E. K., 2007, ApJ, 663, 948
  10. Irwin M., Hatzidimitriou D., 1995, Mon.Not.Roy.Astron.Soc., 277, 1354
  11. Jardel J. R., Gebhardt K., 2013, Astrophys.J., 775, L30
  12. King I. R., 1966, AJ, 71, 64
  13. Laporte C. F. P., Walker M. G., Peñarrubia J., 2013, MNRAS, 433, L54
  14. Lewis A., Bridle S., 2002, Phys.Rev., D66, 103511
  15. Lianou S., Cole A. A., 2013, A&A, 549, A47
  16. Łokas E. L., 2009, MNRAS, 394, L102
  17. 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
  18. Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS, 283, L72
  19. Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  20. Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  21. Newman A. B., Treu T., Ellis R. S., Sand D. J., 2013, ApJ, 765, 25
  22. Peñarrubia J., Pontzen A., Walker M. G., Koposov S. E., 2012, ApJ, 759, L42
  23. Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS, 430, 105
  24. Pontzen A., Governato F., 2014, Nature, 506, 171
  25. Richardson T., Fairbairn M., 2014, MNRAS, 441, 1584
  26. Shao S., Gao L., Theuns T., Frenk C. S., 2013, MNRAS, 430, 2346
  27. Simon J. D., Geha M., 2007, Astrophys.J., 670, 313
  28. Springel V., Wang J., Vogelsberger M., Ludlow A., Jenkins A., et al., 2008, Mon.Not.Roy.Astron.Soc., 391, 1685
  29. Strigari L. E., 2013, Phys.Rept., 531, 1
  30. Strigari L. E., Frenk C. S., White S. D., 2010, Mon.Not.Roy.Astron.Soc., 408, 2364
  31. Tollerud E. J., Beaton R. L., Geha M. C., Bullock J. S., Guhathakurta P., et al., 2012, Astrophys.J., 752, 45
  32. 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
  33. Walker M., 2013, Planets, Stars and Stellar Systems. Volume 5: Galactic Structure and Stellar Populations, Oswalt, T. D. and Gilmore, G., 1039, p. 1039
  34. Walker M. G., Mateo M., Olszewski E. W., 2009, AJ, 137, 3100
  35. Walker M. G., Mateo M., Olszewski E. W., Penarrubia J., Evans N. W., et al., 2009, Astrophys.J., 704, 1274
  36. Walker M. G., Penarrubia J., 2011, Astrophys.J., 742, 20
  37. White S. D. M., 1981, MNRAS, 195, 1037
  38. Wolf J., Martinez G. D., Bullock J. S., Kaplinghat M., Geha M., et al., 2010, Mon.Not.Roy.Astron.Soc., 406, 1220
  39. Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description