© ESO, 2010

1. Introduction

Impressive observations by Walker et al. (2009a) and Mateo et al. (2008) have accurately catalogued thousands of stars in the vicinity of the eight classical dwarf spheroidal galaxies (dSphs) of the Milky Way. By binning these line-of-sight (los) velocities as a function of projected radius and by computing the root mean squared (rms) velocity within each bin, it is possible, through Jeans analysis, to reconstruct the mass profile with considerably less uncertainty than with only global los velocity dispersions.

By fitting mass models to these dSphs and to the ultra-faints (Belokurov et al. 2007; Koposov et al. 2007; Simon & Geha 2007; Willman et al. 2006) that were discovered with the aid of the SDSS survey, Strigari et al. (2008) have proposed that the mass within 300 pc of every dSph galaxy, over four orders of magnitude in luminosity, lies in a small range near 10⁷   M_⊙. This result is at odds with the naive expectation that more luminous galaxies reside in more massive halos; one therefore needs to conclude that the Strigari et al. (2008) result is of some fundamental significance to the dark matter particle which creates the halos at high redshift, in which the dSph galaxies form, or to feedback processes which control star formation (Koposov et al. 2009). Using high resolution simulations combined with semianalytical models, Macciò et al. (2009) compute both luminosity and the mass within 300 pc M₃₀₀ of galactic satellites in the ΛCDM model. They find that M₃₀₀ ~ 10⁷   M_⊙ can be explained by the narrow range of circular velocities covered by visible satellites at the time of accretion.

Several other analyses have been made that rely on the los velocity dispersions, for instance the calculation of the dSph mass-to-light ratios M/L (the only truly free parameter) in Modified Newtonian Dynamics (MOND) by Angus (2008) and the mass modelling performed by Walker et al. (2009a).

An issue may exist with these studies because they depend crucially on the dSph (and all the stars within the projected tidal radius) being in virial equilibrium, except for some field stars which can and have been rejected from metallicity cuts by Walker et al. (2009a). In general, there is no concrete way to know if a star is bound to the dSph or not. This is especially true of the cold dark matter (CDM) paradigm, since in principle we can have any density profile of dark matter halo which can simply be increased to accommodate a fast star. Nevertheless, a re-analysis of the SDSS discovered dSphs Segue I (Niederste-Ostholt et al. 2009) and Hercules (Adén et al. 2009), that have been performed with alternative techniques to the one presented here, has demonstrated that neither of these two ultra-faint dSphs are as massive as previously claimed, and unbound stars are surprisingly common.

The point is that if the bound stars have a certain rms velocity and just a few unbound stars have rms velocities significantly higher, the unbound stars can dominate the combined rms velocity. Since in Newtonian dynamics galaxy mass is proportional to the second power of the velocity dispersion (and to the fourth power in modified Newtonian dynamics) it is absolutely critical to make sure the stars used in the Jeans analysis are bound.

This is contrary to the mass modelling of spiral galaxies and clusters of galaxies. For instance, the masses of spiral galaxies are usually measured from the velocities of the neutral hydrogen gas which is forced to be on circular orbits and tidal features would leave very obvious signatures. For one thing, escaped gas along the los would have an entirely different density to the gas inside the galaxy. Secondly, warps are relatively easy to identify in gas, but not from the discrete distribution of stars (McConnachie et al. 2009).

The same is true for clusters of galaxies when the mass is inferred from the measured temperature and density of the X-ray emitting plasma. Alternatively, cluster masses can be measured with gravitational lensing (e.g., Clowe et al. 2006) or the caustic method (Diaferio & Geller 1997; Diaferio 1999; Serra et al. 2010), both of which are independent of the state of equilibrium.

Another important consideration is that the central regions of clusters are not tidally influenced by their surroundings, both because of their mass dominance and the time scales involved. Similarly, the majority of spirals evolve in the field, where galaxy-galaxy interaction are relatively rare. At the opposite end of the mass spectrum, dSphs are in close proximity to their hosts (~100 kpc) and the tidal field is rapidly changing as their orbits are shorter (~100 Myr). Furthermore, the dSphs have extremely low surface brightnesses and are very extended.

Since the wrong identification of bound stars can distort the real velocity dispersion of a system, an accurate detection of interlopers is essential. A traditional method for removing interlopers consists of eliminating those candidates whose los velocity is greater than 3σ from the mean velocity. The process is iterated until the number of remaining stars is stable (Yahil & Vidal 1977). A less conservative method is the jackknife that removes, in each step, the stars whose elimination causes the largest change in the virial mass estimator. But these methods are restricted to the information given by the los velocity, and they do not explore the combined information coming from the velocity and the position of the stars.

There are several methods available for performing a membership selection in velocity space. den Hartog & Katgert (1996) developed another iterative method taking into account that the maximum radial velocity of a member can be estimated if the mass distribution of the system is known. The mass profile is calculated by the application of the virial theorem and the stars whose velocities are greater than the maximum velocity allowed at its radius are excluded from the calculation of the mass profile in the next iteration step. This technique has been successfully tested by Klimentowski et al. (2007) in simulated dwarf spheroidal galaxies, and its application to real dwarf galaxies (the same ones as studied here) by Łokas (2009) resulted in the removal of significant numbers of potential interlopers and enabled them to find consistent mass-follows-light models that fit the data.

An alternative to these methods is the caustic technique which utilises the three known phase-space coordinates to delineate caustics (escape speed curves) in the los velocity-projected radius diagram which enclose the bound stars. Although the caustic technique was originally designed to study the mass profiles of galaxy clusters beyond the virial radius (for a review see Diaferio 2009), it can also function as an interloper identifier.

This is of particular significance for the ultra-faint dSphs because lack of stars, greater susceptibility to tides and more severe impact of interlopers aggravate the problem. However, in contrast to the ultra-faint dSphs, the classical dSphs have a highly statistically significant number of stars, which the caustic technique requires to operate effectively. So we focus on them and await more data on the ultra-faint dSphs.

After we have the confirmed stars as members, we recompute the los velocity dispersion as a function of projected radius and perform our Jeans analysis in both MOND (to check whether the M/Ls are consistent with the stellar populations) and in Newtonian dynamics to investigate the claim that the mass within a radius of 300 pc is similar for every dSph.

We make use of the recently published data collected by Walker et al. (2009a) for Sculptor, Fornax, Carina and Sextans, as well as Leo I from Mateo et al. (2008).

2. The caustic method

The caustic technique (Diaferio 1999) was originally used to identify galaxy members and compute mass profiles in clusters of galaxies, from the central region to beyond the virial radius. This technique has proven to be very effective in estimating the mass profiles of galaxy clusters from N-body simulations, independently of the dynamical state of the cluster. Considering that the technique can identify members in potentially any system, we have modified it to be used for interloper removal in dSphs. The principle is the same, i.e., from the velocity diagram (los velocity vs. projected radius from the centre), the technique locates the caustic curves which represent the escape velocity of the system. Therefore, the stars within these caustics are members of the system and they can be used to calculate the velocity dispersion profile.

In the velocity diagram of an astrophysical system, the components (stars or galaxies) populate a region whose amplitude decreases with increasing projected distance r from the centre. This amplitude  was identified by Diaferio & Geller (1997) to be equivalent to the escape velocity along the line of sight, i.e. . To locate the caustics, it is necessary to begin by determining a set of candidate member stars of the dSph. These stars will provide a centre, a size and a value for the los velocity dispersion, once they have been arranged in a binary tree (Fig. 1), according to their pairwise “projected” binding energy (1)In this expression, R_p is the projected separation on the sky between the stars i and j, m_j,i are their masses, Δv is their los velocity difference and G is the gravitational constant.

The binary tree is built as follows: (i) each star can be thought as a group g_α; then (ii) at each group pair g_α,g_β is associated the binding energy E_αβ = min { E_ij } , where E_ij is the binding energy between the star i ∈ g_α and the star j ∈ g_β; (iii) the two groups with the smallest binding energy E_αβ are replaced with a single group g_γ and the total number of groups is decreased by one; (iv) the procedure is repeated from step (ii) until only one group is left.

Figure 1 shows the binary tree of a random sample of stars in Leo I. The branches of the tree connect the root (at the top of the figure) to the leaves, namely the stars (at the bottom of the figure). The vertical positions of the nodes where the tree branches depart, between the root and the leaves, depend on the binding energy (1) of the stars that hang from their branches: the more bound the stars are, the lower the node’s location is.

Fig. 1 Binary tree of a random sample of 200 stars in the field of Leo I. Stars are the leaves of the tree (at the bottom). The thick path highlights the main branch. The horizontal lines show the upper and lower thresholds used to cut the tree. Some nodes are labelled on the left side, with the number of stars hanging from those nodes between parentheses.

The tree automatically arranges the stars in potentially distinct groups with only a single parameter, the star mass, which enters Eq. (1). In the following, we assign the same mass to every star. To get effectively distinct groups and to specifically define the set of candidates, we need to cut the tree at some level. This level sets the node from which the candidate members hang. The level is defined by the “σ plateau”, as described below.

2.1. The σ plateau

We identify the main branch of the binary tree as the branch that emerges from the root and contains the nodes from which, at each level, the largest number of leaves hang (see the bold path in Fig. 1). The leaves hanging from each node x of the main branch provide a velocity dispersion . When walking along the main branch from the root to the leaves,  rapidly decreases due to the progressive loss of stars most likely not associated with the dSph (Fig. 2);  reaches a plateau when most of the stars hanging from the main branch are members of the dwarf galaxy. In fact, the system is nearly isothermal and progressively removing the less bound stars does not affect the value of . Then we arrive at the most bound stars, whose binding energy is very small (as well as their velocity dispersion), so  drops again. The extreme of the σ plateau closest to the root is thus the appropriate level for the identification of the system.

The procedure to identify the two nodes x₁ and x₂ corresponding to the σ-plateau extrema is as follows (see Fig. 3). Consider the density distribution of the population of the ’s provided by the different nodes x. The  that provides the maximum of this density distribution defines the velocity dispersion σ_pl, corresponding to the σ plateau. We then identify the N nodes whose  satisfies the relation . The parameters N and δ clearly depend on each other and determine how well we sample the σ plateau. To identify the final nodes x₁ and x₂ from these N nodes, we first consider the five nodes closest to the root and the five closest to the leaves; the node x₁ and x₂ are the nodes, in the first and second sets respectively, with the minimum deviation from σ_pl. With this procedure, we find that N = 15 provides the nodes x₁ and x₂ with  sufficiently similar to σ_pl and close enough to the extrema of the σ plateau.

The σ plateau only appears when the stellar mass m in Eq. (1) lies within a proper range  [m₁,m₂] , which depends on the dSph considered. Outside this range the plateau is not evident enough and the procedure described in the preceding paragraph does not yield sensible results. In fact, each star is a tracer of the velocity and density fields, namely of the gravitational potential of the system; therefore, the parameter m entering Eq. (1) represents a fraction of the total mass of the system and is not necessarily close to the real mass of a star. The arrangement of the binary tree strongly depends on the parameter m because the binding energy (1) is a balance between the potential energy term (proportional to m²) and the kinetic energy term (proportional to m). The absence of a σ plateau indicates that m assumes an inappropriate value, because the binding energy is dominated either by the potential energy or by the kinetic energy. When m is in the correct range  [m₁,m₂] , the σ plateau appears and we can recognise the bound stellar system. By varying m within this range, we can estimate the systematic errors on the velocity dispersion derived from the caustic technique. Below, we apply the caustic technique for three values of m: m = m₁, m = m₂, and an intermediate value (Table 3).

Fig. 2 Line-of-sight velocity dispersion of the “leaves” (stars) of each node along the main branch of the binary tree shown in Fig. 1. Filled squares indicate the nodes x1 (left) and x2 (right).

Fig. 3 Line-of-sight velocity dispersion of the leaves of each node along the main branch (black circles). The filled red circles are the N nodes. The circles with a yellow center indicate the first and the last 5 nodes of the set N and the blue squares denote the final thresholds.

2.2. Final members

The stars hanging from the node x₁ represent the candidate members of the dSph. They determine the centre of the system. The systemic dSph velocity is the median of their velocities, while the celestial coordinates of the centre are the coordinates of the peak of the two-dimensional density distribution of the stars. To find the peak we compute the 2D density distribution on the sky f_q(α,δ) with the adaptive kernel method described in Diaferio (1999). The centres we determine coincide, within the uncertainties, with those provided in the literature.

In the plane (r,v) (Fig. 4), where r is the projected distance and v the los velocity of each star from the dSph centre, the caustics are the curves satisfying the equation f_q(r,v) = κ. Here f_q(r,v) is the star distribution in the plane, and κ is the root of the equation (2)The function is the mean caustic amplitude within R, ,  ⟨ v² ⟩ ^1/2 is the velocity dispersion of the candidate members and R is their mean projected separation from the centre. The stars within the upper and lower caustics are the dSph members, used to compute σ_los(r).

We can use this set of members as the input of an iterative procedure: we locate the caustics in the velocity diagram, remove the stars outside the caustic, and locate new caustics with only the stars left. We proceed until no star is rejected. For each dwarf, we apply this procedure to the velocity diagram using the outcomes provided by the binary tree corresponding to each stellar mass: the procedure converges after a number of iterations in the range  [3,13] , and the total number of rejected stars is in the range  [3,45] , while the percentage of the rejected stars is  [0.2,5.7] % of the stars within the initial caustics. We conclude that the iterative procedure adds little to the identification of members: both the initial and final caustics and the initial and final velocity dispersion profiles are basically indistinguishable. Nevertheless, we keep the set of stars provided by the iterative procedure as the final members.

Fig. 4 Velocity diagram for Leo I (data from Mateo et al. 2008). The black circles are all the stars in the catalogue, while the circles with a green star are the candidates selected by the binary tree. The red curves are the caustics identifying the members. The blue vertical line represents the mean separation R from the centre.

Fig. 5 Identified members and interloper fraction resulting from the application of the caustic technique to a simulated halo.

Fig. 6 Velocity dispersion along the main branch and redshift diagram of a mock dSph in a dense environment. The different colours represent the groups identified by the binary tree (the green circles are the stars of the main group). The thick solid lines delimit the 3σ range.

We performed extensive tests of the caustic technique as an identifier of interlopers of self-gravitating systems. We used dark matter halos extracted from N-body simulations of a cubic volume of 192 h^-1 Mpc on a side of a standard ΛCDM model (Borgani et al. 2004). We will present our results elsewhere (Serra et al., in prep.). Here, we provide in advance an example of our tests, to show that the caustic technique is indeed rather effective at identifying interlopers in simulated dark matter halos. For the sake of simplicity, we define the members as the particles within the sphere of radius 3 r₂₀₀, where r₂₀₀ is the radius within which the average mass density is 200 times the critical density . Figure 5 shows the cumulative member fraction at increasing radii r of a dark matter halo from our simulation: r is in units of the half-mass radius r_1/2, defined by the relation M(<r_1/2) = M(<r₂₀₀)/2. For this halo, r₂₀₀ = 2.25r_1/2. The caustic technique identifies more than 98% of the real members at all radii (left panel). Moreover, the fraction of interlopers, i.e. the particles misidentified as members, remains below 10% out to 5r_1/2 and increases to 20% only at radii as large as 7r_1/2 (right panel).

We emphasise the robustness of the caustic technique by comparing it to the 3σ-clipping method. Figure 6 shows a dSph embedded in a particularly dense environment. We have built this catalogue by taking the celestial coordinates of the stars in the field of Carina and assigning a los velocity to each star in such a way to create a uniform velocity distribution across the range  [−30,30]  km s^-1 and a uniform background in the velocity space along  [−150,150]  km s^-1. The membership determined by the 3σ-clipping procedure yields a velocity dispersion of  ~80 km s^-1, a factor of four larger than the real velocity dispersion of the main distribution  ~17 km s^-1, while the caustic technique gives  ~20 km s^-1. The worse performance of the 3σ-clipping procedure is a consequence of its assumption of an underlying Gaussian distribution of the line-of-sight velocities. This assumption is unnecessary in the caustic technique, which therefore appears to be less prone to failure when the foreground and background contamination is severe. We have chosen to perform this test, even if this system is unrealistic, to show that the technique yields reliable results in a dense velocity diagram.

The caustic technique was already used to remove interlopers both in real clusters (e.g. Rines et al. 2004, 2005), to investigate the dependence of galaxy properties on environment, and the stellar halo of the Milky Way (Brown et al. 2010), to estimate the velocity dispersion of the stars in the halo.

Fig. 7 Velocity diagrams for our 5 dSphs with the 3 sets of caustics overplotted (see Sect. 3.3). For Leo I: the dark pink diamonds indicate stars that are members of all three caustics, the black circles are stars that are members of no caustics and the blue stars are members of at least one, but not all the caustics. For the remaining dwarfs: the three set of caustics are the dotted, dashed and solid lines, the colours of the symbols indicate the probability p of being members. : red triangles, : yellow circles, : blue squares, : green diamonds according to Walker et al. (2007).

3. Results

For the five dSphs, we separate the member stars identified by the caustic technique into projected radius bins of 50 pc. We then evaluate the standard deviation of the velocities around the systemic velocity, which gives us the los velocity dispersion (losVD) as a function of the projected radius. Figure 7 shows the velocity diagram of the five dSphs and Fig. 8 shows the cleaned losVDs for the star members identified with the intermediate of the three sets of caustics that we locate for the three adopted masses (Table 3). Given that we use the entire sample of stars, in Fig. 7 we plot the probability of membership for each star, according to Walker et al. (2009a). In all cases there are a number of low-probability stars within the caustics; however, if we exclude these stars, the velocity dispersion profiles remain unaffected. In Fig. 8, the losVD as calculated by Walker et al. (2007) are also overplotted for reference.

Fig. 8 Line-of-sight velocity dispersion profiles for the Walker et al. (2007) uncleaned samples (grey circles) and our cleaned samples (red squares). In the left-hand panels we overplot the best fit MOND models (parameters listed in Table 2) and in the right-hand panels we plot the best fit cored dark matter halo models (parameters in Table 1). The two lines either side of the best fit are the 1-σ errors on the ρDM,0 and rDM parameters.

The losVD computed through the caustic technique are considerably smoother and in some cases lower than those calculated by Walker et al. (2007). Moreover, the number of members identified with the caustic technique in each dwarf is always smaller than the number of members identified by the algorithm of Walker et al. (2009b) based on the joint distribution of velocity and the magnesium strength. Our procedure of interloper removal also appears to be more strict than the technique of Klimentowski et al. (2009). According to Fig. 1 of Łokas (2009), who use Klimentowski et al. (2009)’s procedure, the number of members she finds is 5-10% greater than our results for Carina, Fornax and Sculptor; for Sextans, her number of members is 25% larger.

The caustic technique has not significantly altered the magnitude of the losVD for the most distant and luminous dSphs (Leo I and Fornax): for these larger dSphs, it has merely pruned some outlying stars. However, for Sextans and Sculptor, the losVDs have changed considerably.

The typical error bars on each losVD bin in Walker et al. (2007) are  ± 1.5   kms^-1, but we have used fewer bins, thus increasing the number of stars per bin and decreasing the uncertainty on each losVD value. Additionally, the three sets of losVDs identified by the caustic technique, found using different star masses, do not systematically differ by more than  ± 0.5   kms^-1 until the last data point for Carina, Sextans and Sculptor (see Fig. 13). The three caustics of Fornax and Leo I vary even less. Therefore, we take  ± 0.75 kms^-1 as a reasonable error bar for each of our radius bins.

For each of the five dSphs, we performed a Jeans analysis to infer their DM halo properties or to infer their M/L ratios in the case of MOND. Crucially, the Jeans equation is (3)which solves for the radial velocity dispersion, , making use of . From this we use the standard method to cast the radial velocity dispersions into los velocity dispersions. The vital point for our analysis is that g(r) is the only part of the equation to change between the MOND and dark matter scenarios. Using dark matter, we know g(r) = g _∗ (r) + g_DM(r). We discuss the case of MOND in Sect. 3.2.

All the dSph parameters, such as luminosity and distance, are taken from the table of Angus (2008) – where the main observational references are Mateo (1998) and Irwin & Hatzidimitriou (1995).

3.1. Dark halo properties

To demonstrate how quantitatively different the DM halos required to fit the cleaned losVDs are, we used the cored dark halo density profile ρ(r) = ρ_DM,0 [1 + (r/r_DM)²] ^-1.5 (cusped halos can achieve just as good fits as shown in Angus & Diaferio 2009; and discussed in Evans et al. 2009) applied in Angus (2010) (where we refer interested readers to find a detailed description of the models) along with constant, velocity anisotropies (β(r) = β_o). We tried models with radially varying velocity anisotropy, but found no improvement in the fits.

We began by evaluating the χ² on a finely spaced grid of the free parameters: ρ_DM,0 and r_DM in the dark halo scenario, and ρ _∗ ,0 (or M/L) in the case of MOND, as well as velocity anisotropy in both cases. The unreduced χ² is defined by (4)where r_i is the radius of the ith data point, σ_los is the fit and σ_obs is the observed velocity dispersion along the line of sight.

After finding the χ² minimum for each dSph, we re-evaluated the surrounding parameter space with greater resolution, and to determine the 1-σ error on each parameter. This is given by the difference from the best fit value a parameter must be altered in order to increase the χ² by 3.5 in the case of the dark halos with three free parameters. Figure 11 demonstrate how in the case of MOND, the χ² depends on the two free parameters M/L and velocity anisotropy and from this the 1-σ errors are easily deduced.

Fig. 9 The enclosed dark mass within 300 pc as a function of dSph luminosity. The purple dots are from Strigari et al. (2008) and the red squares are our cleaned samples presented in Fig. 8. The empty circles correspond to the fits of Angus (2010) to the uncleaned Walker et al. (2007) data, using the same cored dark halos as we implement here. The points referring to the same dSph are slightly shifted along the L axis for clarity.

Apart from Fornax, all the other dSphs are consistent at 1-σ with any β between −3 and −0.5, with a suitably chosen ρ_DM,0 and r_DM. For higher β (i.e. greater than −0.5) there is always a point where for any combination of ρ_DM,0 and r_DM the Δχ² above the χ² minimum reaches 3.5 (or 1-σ for 3 free parameters). On the other hand, for very negative β, it is not always possible to reach a Δχ² of 3.5. However, only in one case (Carina) does the best fit prefer lower values of β than −3, and even then the difference in χ² between β = −1.3 and −3 is negligible (see Fig. 10). For this reason, in Table 1 the lower (upper) error bar on ρ_DM,0 (r_DM) is not strictly the 1-σ error, but that value at β = −3 (and the noted χ²). Here we are applying a prior that the velocity anisotropy cannot be highly tangentially biased.

The parameters of all the models are given in Table 1 with their 1-σ error. Clearly there is a large range of central DM densities (ρ_DM,0) and DM scale radii (r_DM) which can achieve a good fit to the data and this range is only truncated for large ρ_DM,0 and low r_DM because we force the condition β > −3.

Although the range in ρ_DM,0 and r_DM are quite large, the defined range in enclosed mass is not nearly as vast. In Fig. 9 we plot the enclosed mass within 300 pc for each of our 5 dSphs and compare with the analysis of Strigari et al. (2008). One can see that our sample of five appears consistent with a constant value of 10⁷   M_⊙ within 300 pc although Sextans falls a factor of two short. We also overplot the enclosed mass from fits made by Angus (2010; where the original losVDs of Walker et al. 2007 were fitted with the same DM density used here) and the trend for the enclosed mass to decrease after interloper removal is obvious, except for Leo I and Sextans. The lower M₃₀₀ from Angus 2010 for Leo I is a subtle effect of fitting the outer 4 data points (from Walker et al. 2007), for which a fairly radial β along with a large core, but relatively low central density allow a superior fit. In comparison, our fit to the cleaned sample is more consistent with a smaller core and larger central density. This makes our M₃₀₀ larger, but at 400 pc our enclosed mass drops below his model as expected.

We emphasise that although the dSph may have different velocity anisotropies than presented here, these are the simplest models that can achieve good fits to the data. Interestingly, Wolf et al. (2010) used a technique which enables the measurement of the enclosed mass at a single radius (coincidentally the half-light radius) which is independent of velocity anisotropy. The half-light radius of Carina is close to 300 pc and their value is similar to ours. Finally, Adén et al. (2009) have shown the mass within 300 pc of the ultra low luminosity dSph Hercules is reduced from   ×  10⁷   M_⊙ to   ×  10⁷   M_⊙ using another interloper removal technique, making it incompatible with the common mass scale. For our sample, the M₃₀₀ of Carina is actually reduced to be consistent with 10⁷   M_⊙, but the low mass of Sextans is exacerbated. Excluding Sextans, the other 4 dSphs are consistent at 1-σ with M₃₀₀ = 1.2  ×  10⁷   M_⊙

Fig. 10 The lowest possible χ2 for Carina with a fitted dark halo as a function of velocity anisotropy. Notice the lowest χ2 decreases only very slightly from β = −1.3 to −3.

Fig. 11 The lowest possible χ2 for Fornax in the case of MOND as a function of velocity anisotropy (left panel) and central density (right panel).

3.2. MOND

In Angus (2008) – see also Kosowsky (2010) – the losVDs as measured by Walker et al. (2007) for all 8 classical dSphs were taken without further investigation of interlopers. A Jeans analysis was then performed for each dSph, in the context of Modified Newtonian Dynamics or MOND (Milgrom 1983, 1995). The only difference in the Jeans equation for MOND is that the gravitational acceleration must be modified, the stellar density and velocity anisotropy behave in the standard fashion. In the standard paradigm, the gravitational acceleration of a dSph is the combination of the weak stellar component and the dominant dark matter component. On the other hand, for MOND we must account for both the internal dynamics of the stellar component and the external gravitational field of the Milky Way which determine the MONDian gravity as a function of radius in the dSph.

The formulation we use to define g(r) for Eq. (3) is μ(x)g(r) = g _∗ (r) where . Here, a_o = 1.2  ×  10^-10   ms^-2 is the acceleration constant of MOND, is the MOND interpolating function we choose and g_MW(R_MW) = (170   kms^-1)²/R_MW is the external gravity of the Milky Way at the dSph’s Galactocentric radius assuming the terminal value of the Milky Way’s rotation speed to be 170   kms^-1.

Since MOND does not benefit from the additional freedom of a dark matter halo and has an extra constraint due to the Milky Way’s potential, the velocity anisotropy becomes a fitting function that enables one to uniquely determine the MONDian dynamical mass. Therefore, by fitting velocity anisotropy profiles to the losVDs, as well as accounting for the so called external field effect, Angus (2008) was able to deduce the necessary M/L of each dSph by comparing the dynamical mass with the luminosity.

From here, the required M/L was compared with the expectation from models of stellar evolution given the star formation history of the dSphs in question. It was found that while a sensible range would be a V-band M/L of 1−3, that only half (four out of eight) of the dSphs resided in this range (interestingly they are the most distant and most luminous dSphs), whereas a further two (Carina and Ursa Minor) had M/Ls that were high, but consistent at 1-σ and Draco and Sextans had exorbitantly high M/Ls (>10). It was suggested that these high M/Ls could be the result of tidal interactions with the Milky Way producing interloper stars in the foreground and background of the problematic dSphs.

Following exactly the same method, we have re-analysed the dSphs after removing the interloper stars and have found that it leaves all three more luminous dSphs (Fornax, Leo I and Sculptor) with appropriate M/Ls (Fig. 12), Sextans has a borderline M/L (2.7  ±  0.5), but the least luminous dSph Carina still has a suspiciously high M/L (5.9).

The removal of interlopers by the caustic technique does not fully relieve MOND of its problem with the high M/L of the low luminosity dSphs, regardless of the velocity anisotropy used (best fits given in Table 2). Even more severe caustics cannot prune enough stars to give Carina a significantly lower M/L. As can be seen in Table 2, the MOND M/L for Carina from the 3 other caustics we tried give 4.7 and 5.9. For the low luminosity dSphs to be consistent with MOND, they must either form from unusual initial mass functions or be tidally heated by their proximity to the Milky Way. A crucial point to grasp here is that the idea of dSph tidal heating is unusual in terms of the DM paradigm because tidal forces are too feeble at their current locations. Although there are numerical studies of tidal harassment of dSphs by e.g. Peñarrubia et al. (2008), Klimentowski et al. (2009) and Łokas et al. (2010), they all rely on the dSphs plunging to pericentres that are a factor of three or four nearer the Milky Way (≈ 20 kpc) than their current positions and observed proper motions allow (Piatek et al. 2006, 2007, 2008). Given the observed disk of satellites (Kroupa et al. 2005) and that the Magellanic Clouds are not bound to the Milky Way (Kallivayalil et al. 2006b,a) these orbits are improbable, so other mechanisms should be investigated.

Fig. 12 Inferred M/Ls in MOND using the old data (green asterisks, Angus 2008) and our cleaned samples (blue squares). M/Ls consistent with population synthesis models exist between the red dashed lines. The points referring to the same dSph are slightly shifted along the L axis for clarity.

Furthermore, in the DM paradigm, it is not obvious how stars can be easily stripped out of the dSph (although detailed studies are being made, e.g. Wetzel & White 2010), because the DM halo protects the dSph from tidal disruption, especially at the current distances of the dSphs. On the other hand, in MOND the dSphs are more loosely bound – than their Newtonian counterparts with dark matter – and the internal dynamics vary with distance from the Milky Way (due to the so-called external field effect, e.g. Milgrom 1983, 1995; Brada & Milgrom 2000; Angus 2008). When this external gravity (due to the Milky Way) reaches a similar magnitude to the internal gravity of the dSph, the internal dSph gravity decreases towards the simple Newtonian gravity of the stars and the dwarf and its velocities contract. For this reason, dSphs are more susceptible to tidal disruption in MOND. This idea is currently being studied with high-resolution N-body simulations in MOND (Angus et al., in prep.). Another idea to bear in mind is whether bottom heavy initial mass functions (see Mieske & Kroupa 2008) can give rise to the high M/L ratios. One would have to explain why only the currently nearby, low luminosity dSphs are affected by this.

Fig. 13 A comparison between the los velocity dispersions found using different star masses in the caustic technique as discussed in Sect. 3.3. The red triangles corresponding to the intermediate caustics are supplemented by blue squares and green circles in each radius bin that signify the two other, extremal sets of caustics. The best MOND fit is given to each series of data points.

3.3. Star masses

The star mass mentioned in Sect. 2 is the only free parameter included in the caustic technique. To estimate the response of the caustics to changes in the star mass, we produced three sets of caustics for each dSph by setting three different values of the star mass (as shown in Fig. 13). The highest and lowest values contain the range in which we find a recognisable σ plateau. The explicit values of the stellar mass used for each dwarf are given in Table 3.

The differences in the three sets of caustics for each dSph are not significant (Fig. 7, Table 4), and the number of stars lying within the widest caustics and outside the narrowest one do not cause a large variation in the losVD profiles. This emphasises that the method is robust. Nevertheless, to identify the spread in the context of the mass models, we fitted models to the losVD from each caustic set for all 5 dSphs specifically in the case of MOND because differences in the M/L are more obvious since the range is relatively tight. The three models used for Sculptor, Fornax and Leo are barely distinguishable (see M/L column in Table 2) and yield MOND M/Ls that differ from the main caustics at a much lower level than the errors we find at 1-σ. For Carina and Sextans the differences are moderate, but still lower than the 1-σ errors.

Fig. 14 Fit to the line of sight velocity dispersion of Sculptor in MOND with different velocity anisotropy models, as discussed in Sect. 3.4. The dotted blue line is for the best fit constant velocity anisotropy, β = −0.03 (χ2 = 5), the dashed red line is for a conservative (because an upper limit is imposed upon β) radially varying (χ2 = 3.75), and the solid green line is for an unconstrained (χ2 = 2.6).

3.4. Velocity anisotropy

The velocity anisotropies we initially tested were generalised Osipkov-Merritt models that allow for non-isotropic orbits at small radii (5)where β₀, β_∞ and r_β are free parameters (although −∞ < β < 1).

We found no preference for these models over constant velocity anisotropy, except in the case of Sculptor with MOND. In both cases (MOND and DM halos with our proposed density law), the tangential orbits (more negative β) at larger radii are necessary to fit the shape (in particular the flatness) of the losVD profile. The weightings of the tangential orbits are mild in the DM case (usually anywhere between β = −3 and −0.5) because of the freedom given by the DM halo. In MOND Fornax, Carina and Leo I prefer β < −0.5 (but only Leo I prefers β < −1.0). On the other hand, Sextans prefers almost isotropic orbits.

Sculptor definitely prefers a velocity anisotropy that not only varies with radius, but actually increases from very mildly tangential orbits at the centre (β_o = −0.13) to quite strongly radial orbits by the edge (β_∞ = 0.4). In fact, the χ² actually decreases slightly for even more radial orbits, but we imposed a maximum β_∞ of 0.4. The reason for this (demonstrated in Fig. 14) is to keep the velocity anisotropies from being abnormally radial, especially since the other dSphs prefer negative ones. In comparison, constant β (with β = −0.03) has a χ² of 5, has a χ² = 3.75 and has a χ² = 2.6.

Perhaps the reason for Sculptor’s odd losVD, and hence different velocity anisotropy, may be related to the reason the low luminosity dSphs have high M/Ls. Although Sculptor is relatively luminous, it is also much nearer the Milky Way (~87 kpc) than Fornax and Leo I (138 and 257 kpc respectively) making it potentially more susceptible to tidal effects. However, this can only be tested with numerical simulations.

Łokas et al. (2010), who used Klimentowski et al. (2009)’s technique to remove interloper stars from the catalogues of Walker et al. (2009a) that we also use here, found that reasonable mass-follows-light models (no dark halo is assumed, simply the stars are given arbitrary M/Ls) with constant velocity anisotropies could be fitted to the los VDs of the cleaned sample and their kurtoses. Interestingly, the losVDs emerging from the interloper removal technique used by Łokas et al. (2010) have a general trend to decline with radius, which allows for the good fits with constant velocity anisotropies and mass-follows-light models. On the other hand, four of our five (excluding Sculptor) stay very constant until the last measured point, meaning choice of interloper removal technique is key to constraining which models can be fitted to the data. Nevertheless, our ranges of β are still very similar to Łokas et al. (2010)’s, with the exception of Carina. Our larger ranges in the allowed βs stem from the fact that we use dark halos instead of mass-follows-light models.

4. Conclusion

We have applied the caustic technique to remove star interlopers in the 5 classical dSphs observed by Walker et al. (2007, 2009a) and Mateo et al. (2008). We have applied the technique to each dSph, taking into account the systematics introduced by the free parameter in the binary tree construction, namely the star mass. We have found that there are no major discrepancies between the caustics calculated from each star mass (as shown by Fig. 7 and Tables 1 and 2) in the mass range where a σ plateau appears. Consequently, the velocity dispersion profiles are not strongly dependent on this parameter of the caustic technique.

Using this technique we have shown that unbound/non-equilibrium stars are reasonably common in the low luminosity dSphs of Carina, Sextans and the relatively nearby (but more luminous) Sculptor. Removing these interlopers still leaves enclosed dark halo masses within 300 pc that appear to be constant and independent of the dSph luminosity (Strigari et al. 2008; Strigari 2010), except for Sextans.

We also investigated the implications for MOND, for which original analysis by Angus (2008) showed the M/Ls of the high luminosity dSphs were in perfectly reasonable agreement with stellar population synthesis models, but the low luminosity dSphs had abnormally high M/Ls. In particular, Sextans and Draco cast serious doubt on MOND’s validity. After removing the interloper stars identified by the caustic technique, Sextans’ M/L has dropped to 2.7  ±  0.5, but Carina still has a rather high M/L of 5.9 in solar units. On the other hand, Fornax, Leo I and Sculptor were found to have very similar M/Ls − 1.05, 1.28, 1.85, respectively. Further analysis of the final data release for Draco and Ursa Minor by Walker et al. (2007), and MONDian N-body simulations will help to evaluate this incongruity.

The agreement between our results and previous identifications of the dSph interlopers based on kinematic and/or spectroscopic information show that the caustic technique can be applied with success to both galaxies and clusters of galaxies. The disadvantage of this technique is the large number of stars or galaxies required to obtain robust estimates of the caustic location. However, the independence of the method on the dynamical state of the system and the lack of any assumption on the underlying distribution of the member velocities provide substantial advantages over other widely used techniques for removing interlopers.

Acknowledgments

G.W.A.’s research is supported by Università di Torino and Regione Piemonte. We also acknowledge partial support from the INFN grant PD51.

References

All Tables

All Figures

	Fig. 1 Binary tree of a random sample of 200 stars in the field of Leo I. Stars are the leaves of the tree (at the bottom). The thick path highlights the main branch. The horizontal lines show the upper and lower thresholds used to cut the tree. Some nodes are labelled on the left side, with the number of stars hanging from those nodes between parentheses.
In the text

	Fig. 2 Line-of-sight velocity dispersion of the “leaves” (stars) of each node along the main branch of the binary tree shown in Fig. 1. Filled squares indicate the nodes x1 (left) and x2 (right).
In the text

	Fig. 3 Line-of-sight velocity dispersion of the leaves of each node along the main branch (black circles). The filled red circles are the N nodes. The circles with a yellow center indicate the first and the last 5 nodes of the set N and the blue squares denote the final thresholds.
In the text

	Fig. 4 Velocity diagram for Leo I (data from Mateo et al. 2008). The black circles are all the stars in the catalogue, while the circles with a green star are the candidates selected by the binary tree. The red curves are the caustics identifying the members. The blue vertical line represents the mean separation R from the centre.
In the text

	Fig. 5 Identified members and interloper fraction resulting from the application of the caustic technique to a simulated halo.
In the text

	Fig. 6 Velocity dispersion along the main branch and redshift diagram of a mock dSph in a dense environment. The different colours represent the groups identified by the binary tree (the green circles are the stars of the main group). The thick solid lines delimit the 3σ range.
In the text

Fig. 7 Velocity diagrams for our 5 dSphs with the 3 sets of caustics overplotted (see Sect. 3.3). For Leo I: the dark pink diamonds indicate stars that are members of all three caustics, the black circles are stars that are members of no caustics and the blue stars are members of at least one, but not all the caustics. For the remaining dwarfs: the three set of caustics are the dotted, dashed and solid lines, the colours of the symbols indicate the probability p of being members. : red triangles, : yellow circles, : blue squares, : green diamonds according to Walker et al. (2007).

In the text

Fig. 8 Line-of-sight velocity dispersion profiles for the Walker et al. (2007) uncleaned samples (grey circles) and our cleaned samples (red squares). In the left-hand panels we overplot the best fit MOND models (parameters listed in Table 2) and in the right-hand panels we plot the best fit cored dark matter halo models (parameters in Table 1). The two lines either side of the best fit are the 1-σ errors on the ρDM,0 and rDM parameters.

In the text

Fig. 9 The enclosed dark mass within 300 pc as a function of dSph luminosity. The purple dots are from Strigari et al. (2008) and the red squares are our cleaned samples presented in Fig. 8. The empty circles correspond to the fits of Angus (2010) to the uncleaned Walker et al. (2007) data, using the same cored dark halos as we implement here. The points referring to the same dSph are slightly shifted along the L axis for clarity.

In the text

	Fig. 10 The lowest possible χ2 for Carina with a fitted dark halo as a function of velocity anisotropy. Notice the lowest χ2 decreases only very slightly from β = −1.3 to −3.
In the text

	Fig. 11 The lowest possible χ2 for Fornax in the case of MOND as a function of velocity anisotropy (left panel) and central density (right panel).
In the text

	Fig. 12 Inferred M/Ls in MOND using the old data (green asterisks, Angus 2008) and our cleaned samples (blue squares). M/Ls consistent with population synthesis models exist between the red dashed lines. The points referring to the same dSph are slightly shifted along the L axis for clarity.
In the text

	Fig. 13 A comparison between the los velocity dispersions found using different star masses in the caustic technique as discussed in Sect. 3.3. The red triangles corresponding to the intermediate caustics are supplemented by blue squares and green circles in each radius bin that signify the two other, extremal sets of caustics. The best MOND fit is given to each series of data points.
In the text

Fig. 14 Fit to the line of sight velocity dispersion of Sculptor in MOND with different velocity anisotropy models, as discussed in Sect. 3.4. The dotted blue line is for the best fit constant velocity anisotropy, β = −0.03 (χ2 = 5), the dashed red line is for a conservative (because an upper limit is imposed upon β) radially varying (χ2 = 3.75), and the solid green line is for an unconstrained (χ2 = 2.6).

In the text