N-body simulations of the Carina dSph in MOND

Abstract

The classical dwarf spheroidals (dSphs) provide a critical test for Modified Newtonian Dynamics (MOND) because they are observable satellite galactic systems with low internal accelerations and low, but periodically varying, external acceleration. This varying external gravitational field is not commonly found acting on systems with low internal acceleration. Using Jeans modelling, Carina in particular has been demonstrated to require a V-band mass-to-light ratio greater than 5, which is the nominal upper limit for an ancient stellar population. We run MOND N-body simulations of a Carina-like dSph orbiting the Milky Way to test if dSphs in MOND are stable to tidal forces over the Hubble time and if those same tidal forces artificially inflate their velocity dispersions and therefore their apparent mass-to-light ratio. We run many simulations with various initial total masses for Carina and Galactocentric orbits (consistent with proper motions), and compare the simulation line-of-sight velocity dispersions (losVDs) with the observed losVDs of Walker et al. We find that the dSphs are stable, but that the tidal forces are not conducive to artificially inflating the losVDs. Furthermore, the range of mass-to-light ratios that best reproduces the observed losVDs of Carina is 5.3 to 5.7 and circular orbits are preferred to plunging orbits. Therefore, some tension still exists between the required mass-to-light ratio for the Carina dSph in MOND and those expected from stellar population synthesis models. It remains to be seen whether a careful treatment of the binary population or triaxiality might reduce this tension.

methods: numerical, galaxies: dwarf, galaxies: kinematics and dynamics, Local Group, dark matter

1 INTRODUCTION

The classical dwarf spheroidal (dSph) galaxies of the Milky Way (MW) are eight low surface brightness galaxies that are currently at distances between 60 and 250 kpc. They have total luminosities in the V band ranging from L_V ∼ 4 × 10⁵ to 1.7 × 10⁷ L_⊙ (Mateo 1998) and sizes of order a kiloparsec. For comparison, the MW luminosity and size are L_V ∼ 6 × 10¹⁰ L_⊙ (McGaugh 2008) and ∼30 kpc. Clearly, such puny luminosities within relatively large volumes earn the dSphs their low surface brightness moniker and also put them in a very interesting category since low surface brightness galaxies typically have large dark matter (DM) components.

Being spheroidal systems, information about their dynamical mass can be obtained from Jeans modelling of their stellar velocity dispersions (see Mamon & Boué 2010 for more information). For this reason, Walker et al. (2007) obtained hundreds of spectra of probable member stars for each of the dSphs, sampled over their full projected areas. Photometrically and spectroscopically identified interloper stars (non-members, typically foreground stars) were rejected and each dSph's projected velocity dispersion, as a function of projected radius, was computed. They then performed Jeans modelling of each dSph, which employs the observed stellar surface brightness profile and fits for the unknown DM profile, by comparing modelled with observed projected velocity dispersions. This blatantly showed that the dSphs are some of the most DM dominated (in Newtonian dynamics) galaxies in the Universe.

Although the dynamics of the dSphs can be easily explained by the presence of DM, there are other peculiarities related to their phase-space distribution around the MW which makes one question this conclusion. The major open questions relating to dSphs are comprehensively reviewed in Walker & Loeb (2014), but we restate them here. First of all, from comparison with cold dark matter (CDM) only cosmological simulations (like those of Klypin et al. 1999; Moore et al. 1999) one would naively expect a greater number of these satellite galaxies within 250 kpc of the MW. Certain authors like Benson et al. (2002), Muñoz et al. (2009), Macciò et al. (2010) and Li, De Lucia & Helmi (2010) have suggested that this lack of satellites may be due to star formation inefficiencies due to re-ionization and supernova feedback in these lower mass CDM haloes which only enables a fraction of all haloes to form stars. However, this fails to address the problem noted by Boylan-Kolchin, Bullock & Kaplinghat (2012) that associating the dSphs with the most massive MW subhaloes, as we expect in these models, is incompatible with the relatively low masses and densities of the measured DM haloes.

The other more pressing concern is that the dSphs are not isotropically distributed around the MW. Rather, they are distributed as a great rotationally supported disc that is surprisingly thin, with an rms thickness of 10–30 kpc (see Metz, Kroupa & Libeskind 2008 and the detailed review of Kroupa et al. 2010), which is substantially smaller than typical rms thicknesses in nearby groups of galaxies. If it were an isolated incident, this would be less troubling, but Ibata et al. (2013) have recently shown a similar structure in the satellite galaxy distribution surrounding the M31 galaxy with an rms thickness of less than 14.1 kpc (with 99 per cent confidence) to which half the satellites belong. Furthermore, Chiboucas et al. (2013) have recently identified a flattened distribution of satellites around M81.

These satellite distributions have been shown to be highly unlikely to arise from CDM cosmological simulations, although once in place they could naturally be stable (Angus, Diaferio & Kroupa 2011; Deason, Belokurov & Evans 2011; Pawlowski et al. 2012; Bowden, Evans & Belokurov 2013).

On the other hand, following a merger or a flyby (e.g. Zhao et al. 2013) between two galaxies, with mass ratios between 1:1 and 1:4, the probability of forming such a polar disc of satellites could easily reach 50 per cent (Pawlowski et al. 2012).

Separately, there are observations of dwarf galaxies forming out of the tidal debris produced from a wet galactic merger (Bournaud et al. 2007), which may demonstrate evidence for Modified Newtonian Dynamics (MOND; Gentile et al. 2007; Milgrom 2007).

Returning to the dSphs, if they are in fact tidally formed they should not have large DM abundances. Furthermore, they have very little neutral hydrogen (Mateo 1998) and no significant emission from molecular gas. However, these eight classical dSphs do require large DM abundances when interpreted with Newtonian dynamics, and they have a peculiar orbital distribution that may be difficult to explain within the CDM framework. Therefore, it is worth investigating their dynamics in an alternative theory of gravity that can, in principal, be consistent with the merger scenario and the large velocity dispersions without galactic DM. One such alternative is MOND (Milgrom 1983 and see Famaey & McGaugh 2012 for a thorough review).

Brada & Milgrom (2000b) used a particle-mesh N-body solver to study the influence of the MW on the dSphs. Their work preceded the high-quality velocity dispersion data, but demonstrated that there are orbital regions where dSphs can orbit with adiabatic (reversible) changes to their velocity dispersion and density profiles. In addition, there are non-adiabatic regions where the rapid change of the external gravitational field of the MW disturbs the density profile at pericentre and this does not recover by the time the dSph returns to apocentre. Finally, there are tidal regions where mass will be stripped from the dSphs at pericentre.

Using the data of Walker et al. (2007): Angus (2008) and Serra, Angus & Diaferio (2010) performed Jeans modelling in MOND. There, the goal was to isolate the two free parameters: the mass-to-light ratio of the stellar population and the velocity anisotropy. Velocity anisotropy is the a priori unknown relationship between the probability of radial and tangential stellar orbits within the dSph. This can also be used as a free parameter in the context of DM halo fitting, but is somewhat redundant given the freedom of possible DM halo profiles. In MOND, it is an essential ingredient to alter the shape of the projected velocity dispersion profile, whereas all the mass-to-light ratio can do is raise or lower the amplitude of the velocity dispersions. Angus (2008) found that the four dSphs with the highest surface brightness (highest internal gravities) had reasonable mass-to-light ratios, but the other four required mass-to-light ratios that were larger than the expected range of 1 to 5 in the V band found from stellar population modelling (Maraston 2005).

Much simulation work has been done in this vein in the standard paradigm (see e.g. Kroupa 1997; Klessen, Grebel & Harbeck 2003; Read et al. 2006; Klimentowski et al. 2009; Peñarrubia et al. 2009). More specifically, the work of Muñoz, Majewski & Johnston (2008) focused on a very similar thesis as ours, which was whether tidally disturbed mass follows light models of a DM dominated Carina dSph are consistent with the observed projected surface density and projected velocity dispersion profile. Those authors found that there were indeed combinations of mass and orbital parameters that could faithfully reproduce the Carina dSph.

Sánchez-Salcedo & Hernandez (2007) investigated the likelihood of survival for the dSphs in MOND after successive orbits over a Hubble time. They found that only Sextans was likely to dissolve in less than a few Gyr, but that the deduced dynamical mass-to-light ratios of Ursa Minor and Draco (out of the eight classical dSphs) were too large to be consistent with only the stellar populations. They also showed, based on their current positions, that tidal stirring might be an important consideration for Sextans, Sculptor and Ursa Minor, but not Carina. Other relevant work was carried out by Sánchez-Salcedo & Lora (2010) and Lora et al. (2013) who looked at the importance of cold kinematic substructures that are found in the Sextans and Ursa Minor dSphs. It was shown that their longevity can be used to discriminate between modified gravity and CDM.

Given the separation in surface brightness between dSphs that satisfied MOND and those that did not, it was suggested in Angus (2008) that the latter four dSphs may be subject to tidal forces that produce tidally unbound interloper stars and inflate the velocities of the bound stars.

Our aim here is to test this hypothesis by running high-resolution MOND N-body simulations of satellite galaxies orbiting the MW and comparing the simulated projected velocity dispersions with the observed ones. Insodoing we also hope to elucidate the zones of possible orbits open to the satellites without being torn to shreds by the MW. This is an essential sanity check for when high accuracy proper motions become available.

We focus on the Carina dSph because out of the four least luminous classical dSphs it has a well-measured surface brightness profile, large numbers of stellar line-of-sight velocities for Jeans modelling and relatively accurately measured proper motions.

In Section 2 we present the Jeans analysis, in Section 3 we discuss how to incorporate the external field and the setup of our simulations. In Section 4 we compare simulated with observed projected velocity dispersions, in Section 5 we give our results and finally in Section 6 we draw our conclusions.

2 JEANS ANALYSIS

2.1 Likelihood analysis of the Carina dSph's observed projected velocity dispersion profile

Dynamical interloper stars, i.e. those not identified photometrically or spectroscopically, can be dealt with using various algorithms (none of which is fully accurate). At least two such methods have been applied to the dSphs (Klimentowski et al. 2007; Serra et al. 2010) and removed significant numbers of presumed interlopers. These interlopers are typically tidally stripped stars that are mostly removed during pericentre. In Section 4, we will be comparing our simulation velocity dispersions with the observed ones. This observational data is expected to contain dynamical interlopers, and our simulations – if valid representations of the real dSphs – should therefore produce comparable numbers of interlopers (assuming they originate from the dSph and not the MW).

For the above reason, we re-bin data from Walker et al. (2007) without any removal of dynamical interlopers beyond that performed by those authors. We separated the data into projected radius bins of 50 pc with unequal numbers of stars per bin. The projected velocity dispersion and its associated uncertainty, in each bin, is calculated by making a Markov chain Monte Carlo (MCMC) analysis using equation B1 of Hargreaves et al. (1994)

\begin{equation} P(v_i)={1 \over \sqrt{2\pi (\sigma _i^2+\sigma _v^2)}}{\rm exp}\left[{-v_i^2 \over 2(\sigma _i^2+\sigma _v^2)} \right]. \end{equation}

(1)

This equation uses the different measurement errors of each individual star to weight their contribution to the line-of-sight velocity dispersion in each radius bin. It gives the probability of a sample of stars with zero systemic velocity, with projected velocity, v_i, and velocity uncertainty, σ_i, being chosen from a distribution with velocity dispersion, σ_v.

The likelihood is formed by the product of these probabilities for the number of stars in that bin. The MCMC analysis fits for the velocity dispersion and uncertainty in each bin by sampling the likelihood and producing a cumulative likelihood distribution. This allows us to ascertain the maximum likelihood velocity dispersion and 1σ error for each radius bin.

2.2 Approximate Jeans equation for MOND

The Jeans equation for a spherical galaxy solves for the radial velocity dispersion, σ_r, requiring knowledge of the logarithmic density slope |$\alpha (r)={{\rm d}\ln \rho _* \over {\rm d}\ln r}$| of the tracers (stars in the dSph's case), velocity anisotropy |$\beta =1-{\sigma ^2_t \over 2\sigma ^2_r}$| – where σ_t is the tangential velocity dispersion. It also assumes knowledge of the gravity profile, g(r), which is usually based on fitting for the unknown parameters of the DM halo, or the mass-to-light ratio in non-DM models. This gives

\begin{equation} {{\rm d} \over {\rm d}r}\sigma _r^2(r) + {\alpha (r)+2\beta \over r}\sigma _r^2(r) = -g(r). \end{equation}

(2)

Using geometrical arguments, it is possible to convert the radial velocity dispersion into a line-of-sight (or projected) velocity dispersion, which is used to compare with the observed velocity dispersion.

2.3 Likelihood analysis of the Carina dSph's mass-to-light ratio and velocity anisotropy in MOND

In order to re-emphasize the likelihood of the two free parameters in the MOND Jeans analysis of the Carina dSph galaxy, we performed another MCMC analysis. The modelled, projected velocity dispersion is a function of the surface brightness profile, stellar mass-to-light ratio, stellar velocity anisotropy, Galactocentric distance and MW mass (see Section 2.2 and equations 2 and 6 for the Jeans equation relating these parameters). We consider all parameters fixed except the constant (with radius) velocity anisotropy and mass-to-light ratio. These two we allow to vary, and we produce likelihood plots after comparing the modelled projected velocity dispersion with the observed one.

We have fixed the baryonic mass of the MW to be M_MW = 6 × 10¹⁰ M_⊙ (McGaugh 2008) and use a MW–Carina distance of r_MW = 101 kpc (Mateo 1998). We use a 3D light distribution of the form ρ_*(r) = M/L_V × ρ_{*, 0}(1 + (r/r_c)²)^−α which we fitted to the observed surface brightness profile of Carina from Muñoz et al. (2006) in Fig. 1. We use ρ_{*, 0} = 3.8 × 10^{− 3} M_⊙ pc^{− 3}, r_c = 410 pc and α = 3.5 to give an M/L_V = 1 total luminosity of L_V = 4.4 × 10⁵ L_⊙ (Mateo 1998, for which the uncertainty is ∼27 per cent). We also use equation (6) to include the external field effect of MOND, the appropriateness of which we discuss in Section 3.

Figure 1.

Here we plot the projected number density of stars in Carina taken from Muñoz et al. (2006). The solid line is the re-normalized projection of our best-fitting 3D stellar density model with ρ_*(r) = M/L_V × 3.8 × 10^{− 3}[1 + (r/410 pc)²]^{− 3.5} M_⊙ pc^{− 3}.

Open in new tab Download slide

In the left-hand panel of Fig. 2 we plot the probability distribution functions for Carina's mass-to-light ratio and velocity anisotropy using MCMC Jeans modelling of the projected velocity dispersion profile. One can see there is a strong preference for a mass-to-light ratio larger than 3 and the maximum likelihood with 1σ uncertainty is 5.2 ± 1.2. The most probable velocity anisotropy is β = −0.8 (see right-hand panel of Fig. 2) and isotropic (as well as radial) orbits are disfavoured. If we add Carina's distance as a free parameter, with a prior set by its observational error, we find no significant change in the maximum likelihood solution for the mass-to-light ratio. However, the 1σ confidence limits increase by 15 per cent. According to Maraston (2005), the upper limit for a mass-to-light ratio in the V band for an old population of stars is 5. Therefore, the modelled mass-to-light ratio is at the high end of the expected range. Furthermore, a mass-to-light ratio of 5 requires a particular initial mass function and Carina is not formed purely by an ancient stellar population. In fact, Mateo (1998) shows there are several star-forming epochs, with one strong burst around 6 ± 1 Gyr ago. Additionally, Ursa Minor, Draco and Sextans appear to require higher than expected mass-to-light ratios.

Figure 2.

Here we plot the normalized probability distributions of the MOND stellar mass-to-light ratio (left-hand panel) and velocity anisotropy, β (right-hand panel) from a Jeans modelling MCMC analysis of the projected velocity dispersions measured by Walker et al. (2007).

Open in new tab Download slide

3 INCORPORATING THE EXTERNAL FIELD

MOND is an alternative theory of gravity which removes the need for galactic DM by appealing to stronger than Newtonian gravity in regions of weak acceleration. Its phenomenological basis lies in the constancy of the outer parts of rotation curves of spiral galaxies, the baryonic Tully–Fisher relation and the apparent one-to-one correspondence between baryonic density and the observed dynamics of the large majority of galactic systems. This last point means that knowledge of the baryonic matter distribution gives full knowledge of the DM distribution when Newtonian dynamics are used. In MOND, an isolated spherical galaxy with Newtonian internal gravity g_{n, i} = GM(r)r⁻² will actually produce a gravitational field according to |$g_{i}=\nu (|g_{n,i}|/a_o)g_{n_i}$|⁠, where ν is an interpolating function that allows the gravitational field to transition from the Newtonian dynamics we experience in the Solar system, or the bright nuclei of galaxies, to the necessarily amplified gravity at the edges of spiral galaxies.

In this work, we have chosen the interpolating function

\begin{equation} \nu (y)=0.5+0.5\sqrt{1+4/y} \end{equation}

(3)

(see Famaey & Binney 2005; Famaey & McGaugh 2012 for a discussion). This factor ν is therefore greater than or equal to unity and it will become larger the weaker the internal gravity of the system is. This means that the central regions of dense elliptical galaxies will have ν ≈ 1, MW like spiral galaxies will have ν of around 1.5–2 near the Sun's position and rising thereafter. Finally, galaxies with very low stellar densities will have large ν, where |$\nu \rightarrow \sqrt{a_o/|g_{n,i}|}$| which leads directly to flat rotation curves and the baryonic Tully–Fisher relation.

There is of course one final complication, which is the result of MOND breaking the strong equivalence principle (Milgrom 1986). In MOND, the internal gravity of a satellite galaxy of the MW is determined not only by the satellite's stellar distribution (and therefore mass distribution), but also by the local strength of the MW's gravitational field. This effect should not be confused with tidal forces, which also exist.

The tidal force imposes a force gradient across the satellite from the far to near side (relative to the MW). But in MOND, if the magnitude of the MW's gravitational field is constant across the satellite, the satellite's internal gravity will still be altered by that constant field. This phenomenon does not exist in Newtonian dynamics. Increasing the magnitude of the MW's gravitational field, by locating the satellite closer to it, will reduce the internal gravitational field of the satellite towards the Newtonian limit, i.e. ν → constant. Decreasing the magnitude of the MW's gravitational field, by moving the satellite further from it, will increase the internal gravitational field of the satellite towards the MONDian limit. The crucial qualitative corollaries are explored next and can be gleaned from a comparison with a satellite galaxy in Newtonian dynamics.

In Newtonian dynamics a satellite galaxy must have a DM halo that outweighs the stars by a factor between 10 and 100. It is this DM halo that provides the supplementary force to bind the stars, without which the satellite would dissolve in a few dynamical times. Furthermore, the MW has a DM halo which provides the enhanced gravity, over the stars, to bind the satellite galaxy to it.

In MOND, the satellite galaxy has no DM halo, but the boosted gravity of MOND provides the supplementary force to bind the stars. The MW also has no DM halo, but again MOND – sourced by the MW's baryons – provides the supplementary gravity to bind the satellite to the MW.

To illustrate the difference between the dSphs in MOND and Newtonian dynamics, let us say the satellite makes an orbit from apocentre to pericentre. In Newtonian dynamics, the satellite starts the orbit at apocentre with a particular internal gravitational field that drives the large velocity dispersions and at pericentre this gravitational field is virtually unchanged since the DM halo distribution has barely changed. The only difference at pericentre is that tidal forces from the MW start to influence it.

In MOND, the satellite starts the orbit at apocentre with strong internal gravity, but as it approaches pericentre the internal gravity becomes progressively weaker because of the increasing external gravity of the MW, which diminishes the internal gravity of the satellite. In addition, the tides get stronger the closer the dSph gets to pericentre. Therefore, precisely at the moment when the satellite requires additional gravity to protect it from the tidal forces, it loses it, making observations of tidal harassment of the Carina dSph by Battaglia et al. (2012) a natural expectation. This makes satellite galaxies in MOND far more susceptible to tidal destruction than those in Newtonian dynamics with DM. As such, the dearth of satellite galaxies near to the MW might be no surprise in MOND, but many more satellites are expected in CDM.

In the remainder of this section we discuss various methods for including the external field of a host galaxy acting upon a satellite galaxy in MOND. There are three main ideas we introduce, the first two being only approximations. The first one is the simple case where the external field only enters the ν function (equation 3). This approximation of the external field can be defined by equation (6) and we only use this in the Jeans analysis of Section 2.3 and to compare the gravitational field of a test satellite with the other methods – it is not used in any simulations.

The second method is to include the external field through the boundary conditions of the particle-mesh Poisson solver. This is described in Section 3.2 but, as per the first case, this method is only used to make comparisons of the gravitational field of a test satellite with the other two methods.

The final method is the one we use in all the simulations. We call this two-component simulations because we not only simulate particles representing the satellite (dSph–Carina), but also we simulate a coarse representation of the host (MW). These two-component simulations can account for both the external field and tidal forces and we describe it in more detail in Section 3.3.

We compare the gravitational fields produced by these three methods in Section 3.4.

3.1 Approximated MOND external field prescriptions for satellite galaxies

In this section we present the well-known equations that govern the incorporation of the external field in MOND. The field equation to solve for the MOND potential (Milgrom 2010), Φ is

\begin{equation} {\vec\nabla}\cdot ({\vec\nabla}\Phi ) ={\vec\nabla }\cdot \left[ \nu (|{\vec\nabla }\Phi _n| / a_o) {\vec\nabla }\Phi _n \right], \end{equation}

(4)

a_o is the MOND acceleration constant chosen here to be 3.6(km s^{− 1})²pc^{− 1}. Φ_n is the Newtonian potential which is solved from the matter density, ρ, using the Poisson equation |${\vec\nabla} ^2\Phi _n=4\pi G \rho$|⁠.

If we are considering a satellite in orbit of a host galaxy, we can separate |${\vec\nabla}\Phi$| into an internal, |${\vec\nabla}\Phi_i$|⁠, and an external, |${\vec\nabla}\Phi_e$|⁠, gravitational field. This gives, after removing the divergences and ignoring the curl field

\begin{eqnarray} &&{{\vec\nabla }\Phi _i = \nu \left({ |{\vec\nabla }\Phi _{n,i}+ {\vec\nabla }\Phi _{n,e}| \over a_o}\right) {\vec\nabla }\Phi _{n,i}}\nonumber\\ &&\,\,\,\,\, +\;\left(\nu \left({| {\vec\nabla }\Phi _{n,i}+ {\vec\nabla }\Phi _{n,e}| \over a_o}\right)-\nu \left({ |{\vec\nabla }\Phi _{n,e} |\over a_o}\right) \right){\vec\nabla }\Phi _{n,e}. \end{eqnarray}

(5)

If we now crudely consider only directions in the plane perpendicular to the external field we can ignore the second term in equation (5) because the external gravitational field's magnitude in that direction is zero. However, in the argument of the ν function for the first term of equation (5) we must include the modulus of all gravitational fields regardless of direction. In the direction perpendicular to the external field we can add the external field and the internal field in quadrature in the argument of the ν function, which gives

\begin{equation} {\vec\nabla }\Phi _i=\nu \left({ \sqrt{|({\vec\nabla }\Phi _{n,i})^2+ ({\vec\nabla }\Phi _{n,e})^2|} \over a_o}\right) {\vec\nabla }\Phi _{n,i}. \end{equation}

(6)

Thus, in the limit |${\vec\nabla}\Phi_{n,i} \gg {\vec\nabla}\Phi_{n,e}$|⁠, we can ignore the external gravitational field, |${\vec\nabla}\Phi_{n,e}$| and we reduce equation (5) (or equation 6) to the standard MOND formula in spherical symmetry, |${\vec\nabla }\Phi _i=\nu ({| {\vec\nabla }\Phi _{n,i} |\over a_o}) {\vec\nabla }\Phi _{n,i}$|⁠. If on the other hand |${\vec\nabla}\Phi_{n,e} \gg {\vec\nabla}\Phi_{n,i}$| everywhere, then the gravitational field of the satellite is simply a scaled up version of the Newtonian internal gravitational field, i.e. it is as if Newton's gravitational constant has been revised to a value |$G\times \nu ({ |{\vec\nabla }\Phi _{n,e}| \over a_o})$|⁠. As with most things, the interesting cases lie in the middle ground, so we typically use equation (6) if we perform a curl-free analysis including the external field effect.

3.2 Constant external field

On the topic of handling the external field numerically, we modified our openMP parallelized QUMOND galactic Poisson solver code, that was introduced in Angus et al. (2012) to fit the rotation curves of a sample of spiral galaxies, to account for the external field. The code uses a refinement strategy to go from the coarsest grid to the finest grid, each time halving the size of the box and thus doubling the spatial resolution. This strategy was able to handle the difficult boundary conditions of galactic MOND Poisson solvers (see also Brada & Milgrom 2000a,b; Nipoti, Londrillo & Ciotti 2007; Wu et al. 2007; Llinares, Knebe & Zhao 2008). We also introduced the ability to handle multiple populations with different particle masses as used in Angus et al. (2012) and also to update the positions and velocities, giving it the capability to handle evolving simulations. In the code there is a section that computes the QUMOND source density, right-hand side of equation (4), which we have repeated in the appendix.

In that section, the gradients of the Newtonian potential are taken in the x, y and z Cartesian directions to find the Newtonian gravitational field at one-half cell from the node (i, j, k) in all six directions. On top of this, the magnitude of the ν function must be evaluated at each of these six locations. All we do is change |$g_{z_1} \rightarrow g_{z_1}+g_{z,e}$| and |$g_{z_2} \rightarrow g_{z_2}+g_{z,e}$| meaning the external field is exclusively in the z-direction. Here g_{z, e} is the Newtonian external field. This has a knock-on effect for equations (A2)–(A6). One can see in Fig. 3 a comparison of the y-direction (perpendicular to the external field) internal gravitational field of a satellite galaxy when the strength of the external gravitational field increases, up to the limiting point which is the Newtonian field. Clearly, the external field has a huge potential impact on the internal satellite dynamics.

Figure 3.

Here we plot the internal gravity profile of Carina using various strengths for the external field. We use the parameters presented in Section 2.3. The external field is handled using the boundary conditions method described in Section 3.2 (not the two-component simulations). Starting with the top curve and going to the lowest, the strengths of the external field in units of a_o are 10⁻⁴, 10⁻³, 2.5 × 10⁻³, 5 × 10⁻³, 7.367 × 10⁻³ (red), 0.01, 0.02, 0.03, 0.1, 1, 10 and 100. The upper curve is the MOND limit and the bottom curve is the Newtonian limit. The red curve is the gravity profile for the external field Carina is currently experiencing.

Open in new tab Download slide

3.3 Two-component simulations setup

3.3.1 The code

Recall that our simulations are performed using a particle-mesh grid that is centrally refining. The number of particles is fixed at 128³, which is sensible since there are roughly that many stars in the Carina dSph galaxy, given that the luminosity is merely 4.4 × 10⁵ L_⊙ (Mateo 1998). Our coarsest grid has a box length of 4 Mpc in one dimension, with 65 cells per dimension. Our finest grid is 650 pc in length with still 65 cells. Thus, the spatial resolution improves greatly on the smaller grids.

In our case, we used a finest box of length 650 pc and a coarsest box of 4 Mpc. The smallest scale involved is roughly 10 pc – which is the size of the smallest cell on the finest grid, and of the order of the central mean particle/star separation. The largest scale involved is roughly 100 kpc, which is the current distance of Carina from the MW. This means several of the coarsest grids are only used to accurately find the boundary conditions for the finer grids. Particles are advanced in their orbits depending on where they are in real space. If they are within 250 pc of the centre they are advanced using the gravity calculated on the finest grid. Between 250 and 500 pc they are advanced using the second finest grid and between 500 pc and 1 kpc we use the third finest grid, etc.

The particles are separated into 64 equal batches for input/output reasons. This made it convenient to assign 63–64ths of the particles to represent the mass distribution of the Carina dSph and the final 1/64th to represent the MW. For the MW, the spatial information is not carefully set, but rather is just spherically symmetric. The MW's internal velocity information is not used because we use the dSph's frame of reference and so each time step, every particle representing the MW has its velocity and position updated by the same amount to mimic the true orbit the dSph would be executing around the MW. We use the dSph's frame of reference since we need the dSph to be centred on the most accurate part of the code.

3.3.2 The initial conditions

Since the gravity profile of a dSph embedded in an external field is non-axisymmetric, and its velocity distribution is anisotropic, it is not trivial to produce initial conditions that are stable. As an example, the procedure of Brada & Milgrom (2000b) was to generate a spherical King profile from the distribution function for an isolated, isotropic dSph in deep MOND. They then increased the magnitude of the external field gradually until it attained the value the field should have at apocentre. This includes the external field in a shrewd way, but unfortunately the observations of Walker et al. (2007) and Fig. 2 have subsequently shown that the dSphs require anisotropic velocity distributions. Also, King models produce poor fits to the surface brightness of Muñoz et al. (2006). We opted against following the procedure of Brada & Milgrom (2000b) for an anisotropic distribution function because we have no guarantee that the anisotropy will be unaffected by the increasing external field.

We found that we could evolve to a relaxed dSph satisfying the surface brightness distribution of Fig. 1 starting from a spherical distribution with α = 5.5, r_c = 550 pc and β = −4; the internal gravity was found using equation (6) and the radial velocity dispersion as a function of radius was found from solving the Jeans equation (see Section 2.2). As well as the spatial distribution quickly changing as the simulation evolves, the velocity distribution changes too: from β = −4 to ∼−0.8 in a time-scale of 300 Myr. We used the standard rejection-sampling technique to define the positions of the particles within the dSph and assumed a Gaussian distribution for the velocities. We always initially offset the MW from the dSph along the z-axis, typically by +100 kpc, and give the velocity of the MW relative to the dSph along the +y-axis. The x-axis points out of the orbital plane.

We ran the simulations for up to ∼6 Gyr. We felt that simulating the dSphs for two full Galactocentric orbits would be sufficient to demonstrate the impact of the tides. We were also constrained by available computing resources.

3.4 Comparison of external field parametrizations

There are clearly many parametrizations of the external field effect, as we noted at the beginning of this section. In order to demonstrate the differences in gravity profiles between them, we plot a series of curves together in Fig. 4 using the parameters presented in Section 2.3. The left-hand panel is in the z-direction (parallel to the external field) and the right-hand panel is in the y-direction (perpendicular to the external field). In the left-hand panel, the blue line shows the gravity profile for an isolated Carina with no external field. The positive distance from the centre is towards the MW. The thick red curve is the z-direction gravity (towards or away from the MW) using the two-component simulations with dSph and MW. The black line, which lies on top of the thick red line, is for the dSph only simulation with the external field included via the boundary conditions. The dashed green line does not include the curl field, i.e. only solves equation (5). This clearly is a very poor estimation of the gravity profile and should be avoided at all costs.

Figure 4.

Here we plot various computations of the internal gravity profile of Carina, using the parameters presented in Section 2.3. The left-hand panel is in the z-direction (parallel to the external field) and the right-hand panel is in the y-direction (perpendicular to the external field). In the left-hand panel, the blue line shows the gravity profile for an isolated Carina with no external field. The positive distance from the centre is towards the MW. The thick red curve is the z-direction gravity (towards or away from the MW) using the two-component simulations with dSph and MW. The black line, which lies on top of the thick red line, is for the dSph only simulation with the external field included via the boundary conditions. The dashed green line does not include the curl field, i.e. only solves equation (5). In the right-hand panel, the overlapping thick red and black lines are y-direction curves using the two-component simulation and one component with boundary conditions, respectively. The blue line is the gravity profile for an isolated Carina with no external field, which is larger in amplitude at all radii than the other curves. The dashed green line does not solve for the curl field (i.e. only solves equation 6).

Open in new tab Download slide

In the right-hand panel, the overlapping thick red and black lines are y-direction curves (perpendicular to the external field) using the two-component simulation and one component with boundary conditions, respectively. The blue line is the gravity profile for an isolated Carina with no external field, which is larger in amplitude at all radii than the other curves. The dashed green line does not solve for the curl field (i.e. only solves equation 6). Inspecting the left-hand panel's thick red (two components) and black (boundary conditions) lines, the centre of gravity in the z-direction is not at Carina's mass weighted centre.

In summary, using the constant external field, boundary conditions gravitational field will be quite accurate as a probe of the instantaneous dynamics of a dSph, but using the curl-free solution will introduce large errors, especially in the direction of the external field. Regardless, the two-component simulations must be used to account for tidal fields which become more important the closer the dSph approaches the MW.

4 COMPARISON OF SIMULATED WITH OBSERVED VELOCITY DISPERSIONS

4.1 dSph stability

To check the stability of a dSph in isolation, we performed a simulation with no second component (no MW nor external field) which lasted just over 1 Gyr with M/L_V = 5. 1 Gyr is ∼50 dynamical times, and so any severe changes should already be prominent. In Fig. 5 we plot the spherically averaged density profile of the dSph at regular intervals of time, up to just over 1 Gyr. Our first density curve (black line) is made after 300 Myr to ensure the dSph has had time to relax. The density changes only very slightly during this 1 Gyr period, suggesting the dSph is stable. The outer density fluctuates somewhat due to the slow ongoing leakage of particles. We found that time steps of 0.01 Myr were required to reach convergence. Time steps longer than this caused the dSph to dissolve on a 1 Gyr time-scale. Later we show many plots of the surface density as a function of time for the dSph on orbits around the MW.

Figure 5.

Here we plot the evolution in the 3D density profile of an isolated Carina with no external field or simulated MW over the period of 1 Gyr. The initial conditions are described in Section 3.3.2. The black line is the density after 300 Myr and the red lines show evenly spaced evaluations of the density over a 1 Gyr period. The blue line is the final density after 1 Gyr.

Open in new tab Download slide

4.2 Sampling the mock catalogue

Our comparison with the data of Walker et al. (2007) using our simulations is novel. In our re-binned projected velocity dispersion profile there are different numbers of stars per 50 pc bin. The number of stars per the central radius of the bin from 25 to 525 pc are 17, 43, 75, 79, 100, 83, 77, 67, 45, 33 and 29.

We create a mock catalogue of the line-of-sight velocities of a sample of stars in Carina by projecting the 3D velocities of the particles within the simulated dSph along a direction between the simulated dSph and the solar position in the simulated MW.

Note that if this was for external dSphs, like one of the Andromeda dSphs (McConnachie 2012), then this approach would not work since the satellite orbits are not MW centric. It is essential to translate the positions and velocities because the velocity dispersion in the z-direction of the external field is lower than in the orthogonal directions. However, if we just use the untranslated z-direction velocity dispersions the velocity dispersions will not just change with Galactocentric radius, but also angle.

These above procedures were applied to all particles corresponding to the dSph. In the following, we sample a small number of the particles corresponding to the numbers of stars observed by Walker et al. (2007).

Our next step is to randomly sample particles from the ensemble. With each sampled particle we calculate which projected radius bin it belongs to according to its newly translated projected radius, |$R=\sqrt{x^2+y^2}$|⁠. Using this radius bin we add the square of the line-of-sight velocity, |$v_z^2$|⁠, to the accumulated squared velocity dispersions in that bin. Our final condition is that the velocity of the particle with respect to the systemic velocity is less than 30 km s⁻¹ since there is no star in the observed data with a relative velocity larger than this. We continue to sample a random sequence of particles until we have the same number of particles in each bin as we have stars noted at the start of this section. We then divide the summed z-direction squared velocities by the number of stars in each particular bin and then take the square root. Next, we compare this simulated projected velocity dispersion in each bin with the observed projected velocity dispersion in each bin and calculate the |$\chi ^2_{\rm red}$|⁠. We repeat this process 10k times and each time we use a different random sequence of particles.

4.3 Fraction of good fits

In the left-hand panel of Fig. 6 we plot 25 (from the full 10k) random realizations of the simulated velocity dispersion profile of the static model (from Section 2.3, see also Fig. 1), without allowing for any evolution. This is simply to show the variation in the simulated projected velocity dispersions even before tidal forces have influenced the dSph. One can see that at each radius bin there is a spread of velocity dispersions around the mean profile. Therefore, the possibility exists to have significantly larger or smaller velocity dispersions than the average in each bin. However, the presence of 11 radius bins precludes a velocity dispersion profile with a low amplitude from randomly producing a good fit.

$In the left-hand panel, we plot our re-binning of the projected velocity dispersion data from Walker et al. (2007) with blue filled circles. Overplotted with green lines are 25 random realizations of our best-fitting static model of Carina as explained in Section 4. In the right-hand panel, we plot the $\chi ^2_{\rm red}$ distribution of 10k random realizations of the same model. $\chi ^2_{\rm red}$ levels are highlighted with different colours in order to clearly see the number of fits within the ranges $\chi ^2_{\rm red}<1$ (large black dots), $1<\chi ^2_{\rm red}<1.5$ (red dots), $1.5<\chi ^2_{\rm red}<2$ (blue dots) and $\chi ^2_{\rm red}>2$ (small black dots).$

Figure 6.

In the left-hand panel, we plot our re-binning of the projected velocity dispersion data from Walker et al. (2007) with blue filled circles. Overplotted with green lines are 25 random realizations of our best-fitting static model of Carina as explained in Section 4. In the right-hand panel, we plot the |$\chi ^2_{\rm red}$| distribution of 10k random realizations of the same model. |$\chi ^2_{\rm red}$| levels are highlighted with different colours in order to clearly see the number of fits within the ranges |$\chi ^2_{\rm red}<1$| (large black dots), |$1<\chi ^2_{\rm red}<1.5$| (red dots), |$1.5<\chi ^2_{\rm red}<2$| (blue dots) and |$\chi ^2_{\rm red}>2$| (small black dots).

Open in new tab Download slide

In the right-hand panel we plot |$\chi ^2_{\rm red}$| for 10k random realizations of the same static model. |$\chi ^2_{\rm red}=0.9$| is the best-fitting smooth curve to the re-binned data of Walker et al. (2007) (either with MOND or Newtonian gravity and a DM halo). There are, however, a number of realizations for which |$\chi ^2_{\rm red}<0.5$|⁠, demonstrating the potential to have a significantly better match to the data than the smooth fit allows.

The fraction of realizations with |$\chi ^2_{\rm red}<1$| is 0.06. For the rest of the paper we use this as our figure of merit because it gives a robust likelihood of Carina at its given orbital position having a good match to the data. |$\chi ^2_{\rm red}<1$| was simply a suitable number that would yield a statistically useful return after 10k realizations. Using a lower |$\chi ^2_{\rm red}$| threshold would lead to low number statistics and a much higher one, say |$\chi ^2_{\rm red}<1.5$| or 2 (the dashed red and blue lines, respectively, in Fig. 6), would contain information about poorer fits that we are less interested in. It remains that a realization which generally produces a higher number of |$\chi ^2_{\rm red}<0.5$| than another one, will produce a higher fraction of |$\chi ^2_{\rm red}<1$| than that other realization. We are aware that the correlation of the errors in each radius bin means we have overestimated those errors, however, since we only use the |$\chi ^2_{\rm red}$| to compare fits in a relative sense, we do not see this as a problem.

5 RESULTS

Our goal is to determine whether tidal effects from the MW have a meaningful influence on the Carina dSph and if they are conducive to lowering the inferred mass-to-light ratio relative to that garnered from Jeans modelling. It is also crucial to confirm that dSphs like Carina are stable for several orbits with realistic orbital parameters. To accomplish this, we ran a series of simulations of a Carina-like dSph orbiting the MW with initial conditions as described in Section 3.3. We use a variety of total stellar masses for the dSph and orbital parameters with respect to the MW.

The mass-to-light ratio that this stellar mass corresponds to changes with time since the dSph begins to lose a small fraction of particle mass due to tidal stripping as soon as the simulations start. This is because the initial conditions quickly transform to an equilibrium distribution to which not all particles are bound. The mass of our particles of course do not change, but we renormalize their corresponding luminosity while the dSph loses mass in order to fit its observed total luminosity.

The orbital parameters that we discuss refer to the initial radial distance from the MW and the initial velocity: both radially from and tangentially to the MW. Obviously, the closer the pericentric distance from the MW, the more mass will be stripped and the higher the mass the dSph has, the more resilient it will be to tides.

Since Carina currently appears to be approaching its apocentre, with current distance of ∼100 kpc (see Piatek et al. 2003; Metz et al. 2008), this is where we must always make our comparison with the other simulations – even if in the simulation the apocentre is larger than 100 kpc. The reason for this is that increasing the radial distance from the MW decreases the external gravitational field acting on Carina and thus increases the boost to the internal gravity due to MOND. This would allow smaller mass-to-light ratios to be consistent with the observed dynamics than possible at 100 kpc.

5.1 Measured proper motions

The Carina dSph has measured proper motions from the observations and analysis of Piatek et al. (2003). They use the proper motions, estimated from two separate stellar fields within Carina, along with the well-measured line-of-sight velocity to estimate the radial (V_r) and tangential (V_t), with respect to the MW centre, orbital velocities. The first field gives V_r = 18 ± 32 km s^{− 1} and the second field gives V_r = 22 ± 36 km s^{− 1} and are therefore consistent with each other and produce a weighted mean V_r = 20 ± 24 km s^{− 1}. The tangential velocity for the first field is V_t = 40 ± 53 km s^{− 1} and for the second field is V_t = 140 ± 59 km s^{− 1}, which is a considerable difference that barely allows an overlap within the errors. The weighted mean tangential velocity is therefore V_t = 85 ± 39 km s^{− 1}, but not much credence should be given to the formal error since only two, vastly differing, fields have been measured.

Following up on this measurement, Metz et al. (2008) corrected for the advanced charge transfer inefficiencies of the Space Telescope Imaging Spectrograph (Bristow, Piatek & Pryor 2005) and found the updated weighted means to be V_r = 22 ± 3 km s^{− 1} and V_t = 120 ± 50 km s^{− 1}.

We do not consider the radial velocity as being important in our simulations. Our reasoning is that it is small relative to the tangential velocity and the radial velocity only influences the apocentre, whereas the tangential velocity sets the pericentre. Given that only the pericentre sets the impact from tides, we can safely ignore the small radial velocity.

5.2 Fraction of good fits as function of orbit

To investigate the suitability of a given Carina stellar mass and orbital path to matching the projected velocity dispersion data of Walker et al. (2007), we ran three simulations with different orbital parameters for each of the total Carina stellar masses m = 2.2, 2.64 and 3.08 × 10⁶ M_⊙, which correspond to 5, 6 and 7 times the luminosity with mass-to-light ratio of unity. The three different orbital parameters are simply the initial tangential velocity relative to the initial offset along the z-axis of 100 kpc. These were V_y = 125, 150 and 175 km s⁻¹. V_y = 175 km s⁻¹ leads to an almost circular orbit with pericentre of ∼95 kpc and orbital period of ≈2.4 Gyr. V_y = 125 and 150 km s⁻¹ give pericentres of 51 and 71 kpc, respectively, and orbital periods of ≈1.9 and 2.1 Gyr. In Fig. 7 one can see six rows of plots for the m = 2.2 × 10⁶ M_⊙ simulations, where the left-hand, middle and right-hand columns refer to the initial tangential velocities V_y = 125, 150 and 175 km s⁻¹, respectively. The Figs 8 and 9 show the same plots but for the m = 2.64 and 3.08 × 10⁶ M_⊙ simulations.

$Here we plot six different quantities (each row) for a different orbit (each column) using the same initial conditions and the same initial dSph mass (here m = 2.2 × 106 M⊙). The three different orbits are defined by their initial tangential velocity (Vy) which is in the direction perpendicular to the initial separation between the dSph and the MW of 100 kpc: from left to right Vy = 125, 150 and 175 km s−1. (ii) The Galactocentric distance as a function of time is the second row of plots and (i) the top row shows the fraction of random realizations of the projected velocity dispersion which, when compared with the observed one, yield $\chi ^2_{\rm red}<1$ as a function of time. (iii) The third row gives the fraction of projected mass within radial shells of 0.6, 0.8, 1.0, 1.2, 1.5 and 1.8 kpc as a function of time. (iv) The fourth row give the 1D rms size in each of the three directions (x, y and z; where z is the line of sight) and (v) gives the 1D rms velocity. (vi) The final row shows surface density profiles (red lines) for evenly spaced snapshots in time. The normalization is the same for every snapshot. The blue line in panel (r) is the initial surface density at the start of the simulation. Also overplotted (filled circles with error bars) is the observed surface density of stars from Muñoz et al. (2006).$

Figure 7.

Here we plot six different quantities (each row) for a different orbit (each column) using the same initial conditions and the same initial dSph mass (here m = 2.2 × 10⁶ M_⊙). The three different orbits are defined by their initial tangential velocity (V_y) which is in the direction perpendicular to the initial separation between the dSph and the MW of 100 kpc: from left to right V_y = 125, 150 and 175 km s⁻¹. (ii) The Galactocentric distance as a function of time is the second row of plots and (i) the top row shows the fraction of random realizations of the projected velocity dispersion which, when compared with the observed one, yield |$\chi ^2_{\rm red}<1$| as a function of time. (iii) The third row gives the fraction of projected mass within radial shells of 0.6, 0.8, 1.0, 1.2, 1.5 and 1.8 kpc as a function of time. (iv) The fourth row give the 1D rms size in each of the three directions (x, y and z; where z is the line of sight) and (v) gives the 1D rms velocity. (vi) The final row shows surface density profiles (red lines) for evenly spaced snapshots in time. The normalization is the same for every snapshot. The blue line in panel (r) is the initial surface density at the start of the simulation. Also overplotted (filled circles with error bars) is the observed surface density of stars from Muñoz et al. (2006).

Open in new tab Download slide

Figure 8.

As Fig. 7, but for m = 2.64 × 10⁶ M_⊙.

Open in new tab Download slide

Figure 9.

As Fig. 7, but for m = 3.08 × 10⁶ M_⊙.

Open in new tab Download slide

For Figs 7–10, the top row (panels a–c) shows, as a function of time, the fraction of random realizations of the projected velocity dispersion which, when compared with the observed one, yield |$\chi ^2_{\rm red}<1$|⁠. This is the most important plot, which shows whether this particular combination of total stellar mass and tangential velocity will well reproduce the observed projected velocity dispersions (this procedure is discussed in more detail in Section 4 and Fig. 6). The second row (panels d–f) in Fig. 7 gives Galactocentric distance as a function of time for the simulations.

As Fig. 7, but for three odd models. From left to right m = 1.32 × 106 M⊙, Vy = 175 km s−1; m = 3.08 × 106 M⊙, Vy = 90 km s− 1, rmw = 160 kpc; m = 3.96 × 106 M⊙, Vy = 125 km s−1.

Figure 10.

As Fig. 7, but for three odd models. From left to right m = 1.32 × 10⁶ M_⊙, V_y = 175 km s⁻¹; m = 3.08 × 10⁶ M_⊙, V_y = 90 km s^{− 1}, r_mw = 160 kpc; m = 3.96 × 10⁶ M_⊙, V_y = 125 km s⁻¹.

Open in new tab Download slide

Looking specifically at Fig. 8 (m = 2.64 × 10⁶ M_⊙) for the V_y = 175 km s⁻¹ (right-hand column) simulations, one can see from panel (c) that there are not dramatic changes in the fraction of good fits with time because the orbital distance from the MW does not change significantly during the orbit and we know how crucial the Galactocentric radius is to the internal dynamics in MOND (e.g. Fig. 3). On the other hand, the simulations with V_y = 125 and 150 km s⁻¹ (left-hand and middle columns of Fig. 8, respectively) reach significantly smaller pericentres and this means they are exposed to stronger external gravities which cause the internal gravities to be reduced. In panel (c) of Figs 7 and 9 (i.e. for V_y = 175 km s⁻¹), the fraction of good fits appears to vary more with the orbit than for Fig. 8 panel (c). Actually, the magnitude of the change is similar for all three, it is simply that for Fig. 8 (m = 2.64 × 10⁶ M_⊙) the variation is relative to a larger number.

Another important thing is that on the smaller pericentre orbits, more mass is stripped and this is shown in the third row of plots (Figs 7–10 panels g–i) where the projected enclosed mass in several shells is given. The outermost shell is the projected mass within R = 1.8 kpc and should contain most of the bound mass. For the V_y = 125 km s⁻¹ simulation of Fig. 8 (panel g; m = 2.64 × 10⁶ M_⊙), the mass in the outer shell (R < 1.8 kpc) drops from more than 95 per cent to less than 80 per cent after three orbits, but the V_y = 175 km s⁻¹ (panel i) simulation only loses about 2 per cent of the mass in that shell.

The reduction of the internal gravity, due to the varying external field strength, can be clearly seen in the fourth and fifth rows of plots of Fig. 8 (panels j–o) which show the 1D rms sizes in the three orthogonal directions and the 1D rms velocities. For V_y = 175 km s⁻¹, the rms velocity (panel o) in each of the three directions is very constant as are the rms sizes (panel l). However, for V_y = 125 and 150 km s⁻¹, the sizes (panels j and k) and rms velocities (panels m and n) change according to the orbit. This is why the fraction of good fits is also a function of time.

The final row of plots for Figs 7–10 (panels p–r) show the surface density profiles for evenly spaced snapshots in time. The normalization is the same for every snapshot. The blue line in panel (r) is the initial surface density at the start of the simulation. There is very little change in shape for any of the orbits, but the V_y = 175 km s⁻¹ (panel r) simulation is particularly constant and this demonstrates the stability of the dSphs in MOND and our simulations. We show for comparison the surface density of stars found by Muñoz et al. (2006). One final point to take from these bottom three rows of plots in Figs 7–10 is how the z-direction size (panels j–l) and rms velocity (panels m–o) are different to the similar x and y direction rms velocities and sizes. This is obviously due to the direction of the external field which points along the z-direction at every time step because in the post-simulation analysis we rotate our frame of reference such that the z-direction always points towards the MW. This stretching along the external field direction is a well-known effect investigated by Milgrom (1986), Zhao & Tian (2006) and Wu et al. (2008), and will be the topic of future study.

The top rows of Figs 7–10 (panels a–c) show the suitability of the combination of mass and tangential velocity V_y. Here it is important that we only compare the fractions of good fits at apocentre (100 kpc) and clearly the highest fraction of good fits is found using a mass of m = 2.64 × 10⁶ M_⊙ and V_y = 175 km s⁻¹, which easily gives 0.07 – slightly larger than found from the sampling of the isolated model (0.06). Using m = 3.08 × 10⁶ M_⊙ and V_y = 175 km s⁻¹ gives roughly 0.03 and m = 2.2 × 10⁶ M_⊙ and V_y = 175 km s⁻¹ is roughly 0.02 after one orbit and much lower at 0.012 after two orbits. The reason m = 3.08 × 10⁶ M_⊙ does not work as well as m = 2.64 × 10⁶ M_⊙ is because it is too massive and leads to excessively high velocity dispersions ∼6.6 km s⁻¹ in the z-direction (Fig. 9, panel o). This can clearly be seen because the agreement improves at 94 kpc over 100 kpc where the larger external field reduces the velocity dispersion. For m = 2.2 × 10⁶ M_⊙, the dSph is not massive enough and generates too small velocity dispersions ∼5.8 km s⁻¹ in the z-direction (Fig. 7, panel o). Furthermore, even on the nearly circular orbit with V_y = 175 km s⁻¹, it is stripped gradually by tidal forces and this means that the number of good fits (with |$\chi ^2_{\rm red}<1$|⁠) to the observed velocity dispersions decreases during each orbit.

Using m = 2.64 × 10⁶ M_⊙ and V_y = 150 km s⁻¹ (Fig. 8 panel b) has a high fraction of good fits after one orbit (∼0.06), but this drops after each orbit because mass is stripped (panel h) reducing the gravitational field and velocity dispersions. With m = 2.64 × 10⁶ M_⊙ and V_y = 125 km s⁻¹ (Fig. 8 panel g), too much mass is stripped after one orbit and a very low number of good fits is left (panel a). The lower V_ys are also ruled out for m = 2.2 × 10⁶ M_⊙. For m = 3.08 × 10⁶ M_⊙ and V_y = 125 km s⁻¹ (Fig. 9), the correct amount of mass is stripped (panel g) by the end of the second orbit, but the fraction of good fits (panel a) is still quite mediocre (∼0.02). This appears to be a result of the tidal forces changing the velocity anisotropy. To clarify this, in Fig. 11 we plot the projected velocity dispersion using all particles (not just the number of observed stars) for three simulations after roughly 5 Gyr each. The first curve (red) is for the simulations with m = 3.08 × 10⁶ M_⊙ and V_y = 125 km s⁻¹ and is taken after three full orbits. The other two curves both use V_y = 175 km s⁻¹ after two full orbits, but have different masses: m = 2.2 and 2.64 × 10⁶ M_⊙ (black and blue curves, respectively). Clearly, the m = 2.2 and 2.64 × 10⁶ M_⊙ curves have the same shape, but the m = 2.64 × 10⁶ M_⊙ curve has a larger amplitude and significantly better |$\chi ^2_{\rm red}$|⁠; however, the m = 3.08 × 10⁶ M_⊙ curve has a different shape because it has a more isotropic velocity anisotropy. Therefore, the plunging orbit exposes the m = 3.08 × 10⁶ M_⊙ dSph to tides that distort the velocity distribution towards less tangentially biased orbits which makes it slightly less consistent with the observed projected velocity distributions, according to the fraction of good fits and half as likely to produce the observed velocity dispersions.

Here we plot the projected velocity dispersions using all particles of three simulations (m = 2.2 × 106 M⊙ and Vy = 175 km s−1; black), (m = 2.64 × 106 M⊙ and Vy = 175 km s−1; blue) and (m = 3.08 × 106 M⊙ and Vy = 125 km s−1; red) after two, two and three full orbits, respectively. The data points are our re-binned projected velocity dispersions using the data of Walker et al. (2007).

Figure 11.

Here we plot the projected velocity dispersions using all particles of three simulations (m = 2.2 × 10⁶ M_⊙ and V_y = 175 km s⁻¹; black), (m = 2.64 × 10⁶ M_⊙ and V_y = 175 km s⁻¹; blue) and (m = 3.08 × 10⁶ M_⊙ and V_y = 125 km s⁻¹; red) after two, two and three full orbits, respectively. The data points are our re-binned projected velocity dispersions using the data of Walker et al. (2007).

Open in new tab Download slide

Finally, we have three extra simulations for which we plot (in Fig. 10) the relevant quantities, as per Figs 7–9. The m = 1.32 × 10⁶ M_⊙ plot (Fig. 10 panel a) shows that if m = 1.32 × 10⁶ M_⊙ there is no likelihood of being consistent with the observed projected velocity dispersions. The right-hand column is for m = 3.96 × 10⁶ M_⊙ with V_y = 125 km s⁻¹ and also shows (panel a) that at apocentre (100 kpc), there is no likelihood of it being consistent. The central column is for a simulation with m = 3.08 × 10⁶ M_⊙ that does not start at the fiducial 100 kpc, but rather at 160 kpc and plunges to ∼50 kpc. Both on the inbound and outbound sections of the orbit at 100 kpc, the fraction of good fits (panel b) is mediocre at ∼0.015. This is simply because the mass is too high. We do not have a similar simulation for this orbit and m = 2.64 × 10⁶ M_⊙.

Therefore, the most appropriate model appears to be m = 2.64 × 10⁶ M_⊙ and the closer the orbit is to circular the better the match. Conversely, the more plunging the orbit, the more the velocity anisotropy is transformed towards less tangentially biased orbits creating a poorer match. An important point to bear in mind is that this initial mass m = 2.64 × 10⁶ M_⊙ actually means a final mass which is slightly lower due to the stripped stars/particles. The fraction of stars left bound to the dSph depends on the size of the shell we consider, but is somewhere between 90 and 95 per cent. This means a final M/L of between 5.4 and 5.7, which is close to the best-fitting value found using Jeans analysis (Fig. 2). On top of this is the uncertainty in the observed luminosity of Carina. The reason the preferred M/L is slightly larger than the maximum likelihood of 5.2 from Jeans modelling is that the external field causes elongation in the z-direction and likewise causes the velocity dispersion to be smaller than in the two orthogonal directions. This effect is not taken into account in the Jeans modelling.

Battaglia et al. (2012) have demonstrated that the projected stellar distribution of Carina has tidal tails which suggest ongoing harassment of the dSph from the MW. In Fig. 12 we show contours of the projected particle distribution for Carina on a near circular orbit after 5 Gyr. The tidal tails outside the circular isodensity contours of the bound particles lie in the plane of the orbit and occur naturally in MOND even on a near circular orbit at 100 kpc. Whether the same is true for such a distant orbit in CDM simulations remains to be seen.

Here we plot projected density contours for a simulation with m = 2.64 × 106 M⊙ and Vy = 175 km s−1 after 5 Gyr. Up and down is the x-axis (out of the orbital plane), left and right is the y-axis and the z-axis is along the line of sight, respectively.

Figure 12.

Here we plot projected density contours for a simulation with m = 2.64 × 10⁶ M_⊙ and V_y = 175 km s⁻¹ after 5 Gyr. Up and down is the x-axis (out of the orbital plane), left and right is the y-axis and the z-axis is along the line of sight, respectively.

Open in new tab Download slide

6 CONCLUSION

Here we have run a suite of MOND N-body simulations of a dSph like Carina with various total masses (m = 1.32, 2.2, 2.64, 3.08 and 3.96 × 10⁶ M_⊙) and orbital paths around the MW. We have shown that they are stable and long lived on nearly circular orbits at 100 kpc regardless of mass (≥m = 1.32 × 10⁶ M_⊙) and even on orbits that plunge to 50 kpc. However, the model most likely to give a good fit to the observed projected velocity dispersions is one with an initial m = 2.64 × 10⁶ M_⊙, which means an M/L in the range of 5.4 and 5.7 after two orbits (∼5 Gyr). The more circular the orbit, the less disturbed the internal velocity distribution is. This is important because the observations require substantially negative (tangentially biased) velocity anisotropies. After plunging orbits, the velocity anisotropy becomes slightly more radially biased, reducing agreement with the observations. Considering that an M/L in the range of 5.4 and 5.7 is potentially at odds with stellar populations synthesis models, we considered a model with m = 2.2 × 10⁶ M_⊙, which after a single orbit corresponds to an M/L between 4.5 and 4.7. This model has a likelihood of matching the observations that is roughly 3.5 times smaller than the model with M/L between 5.4 and 5.7. This range of mass-to-light ratio is slightly above that found from basic Jeans analysis because the isopotential contours are stretched (see e.g. Milgrom 1986; Zhao & Tian 2006; Wu et al. 2008) in the direction away from the MW (which coincides here with our line of sight) due to the external field effect. This leads to a stretching of the dSph along the line of sight, relative to the plane perpendicular, and a reduction of the velocity dispersions.

As for the compatibility of different orbits, it would appear that after two orbits with initial V_y = 125 km s⁻¹, the lower masses m = 2.2 and 2.64 × 10⁶ M_⊙ are not capable of generating a sizable fraction of good fits. m = 2.2 × 10⁶ M_⊙ would give less than 0.001, m = 2.64 × 10⁶ M_⊙ less than 0.01, but m = 3.08 × 10⁶ M_⊙ would produce roughly 0.03. This is because mass has been stripped leaving the true M/L after two orbits to be somewhere between 5.9 and 6.2. Using m = 2.64 × 10⁶ M_⊙ after only one orbit with V_y = 125 km s⁻¹ gives a fraction of good fits of only 0.015 with a true M/L between 5.1 and 5.4. So the best-fitting M/L for V_y = 125 km s⁻¹ is likely somewhere between these two limits. However, it will probably still be somewhat less likely than the more circular orbits since the tides adversely affect the velocity anisotropy. For the intermediate orbit with V_y = 150 km s⁻¹, m = 2.64 × 10⁶ M_⊙ leads to a fraction of 0.04 good fits after three full orbits with a true M/L of ∼5.3–5.4. Therefore, for V_y ≥ 125 km s⁻¹ the preferred M/L remains fairly constant (5.3–5.7), but obviously on the more plunging orbits mass is more rapidly stripped and thus it is required that the current M/L is in this range, not the initial one. A parallel observation is that the fraction of stripped mass during a period of almost half the age of the Universe is not more than half on any of the simulated orbits. Therefore, it must be the case that the dSph was formed with a mass very close to its current one and this is likely also true in the CDM paradigm.

Although the preferred M/L is between 5.3 and 5.7, there is still a reasonable probability that the M/L is lower than 5. From the various orbits it would seem that even on a near circular orbit, panel (c) of Fig. 7 shows (after one orbit) that M/L ∼ 4.8 is more than three times less likely than the best model. Panel (a) of Fig. 7 suggests that on an orbit with a 50 kpc pericentre, an M/L ∼ 4.5 has an insignificant probability of producing a good fit.

A larger sample of stellar line-of-sight velocities might subdue the errors here to distinguish between different mass-to-light ratios. Therefore, higher precision proper motions, larger samples of stars, ultraprecise photometry for the total luminosity, and more sophisticated and reliable stellar population synthesis models, as well as a full-fledged treatment of binaries for dSphs would be enormously useful for future studies. Another factor that should be built into future studies of Carina is the possibility for triaxiality in the 3D stellar distribution. This must be an important factor because all dSph surface brightnesses are observed to be moderately elliptical (Irwin & Hatzidimitriou 1995).

The authors are indebted to the referee for considerably improving the content and readability of the paper. GWA is a post-doctoral fellow of the FWO Vlaanderen (Belgium). Part of the research was carried out while GWA was a post-doctoral fellow supported by the Claude Leon Foundation. AD acknowledges partial support from the INFN grant Indark (PD51) and from the grant Progetti di Ateneo/CSP TO_Call2_2012_0011 ‘Marco Polo’ of the Universita’ di Torino.

REFERENCES

Angus

G. W.

. ,

MNRAS

2008

, vol.

387

pg.

1481

Month:	Total Views:
February 2017	3
March 2017	2
April 2017	1
May 2017	5
June 2017	1
August 2017	3
October 2017	1
December 2017	9
January 2018	6
February 2018	5
March 2018	6
April 2018	13
May 2018	5
June 2018	10
July 2018	9
August 2018	14
September 2018	3
October 2018	4
November 2018	6
December 2018	8
January 2019	5
February 2019	12
March 2019	12
April 2019	14
May 2019	8
June 2019	3
July 2019	6
August 2019	15
September 2019	2
October 2019	4
November 2019	11
December 2019	1
January 2020	7
February 2020	1
March 2020	3
April 2020	3
May 2020	1
June 2020	8
July 2020	5
August 2020	1

Article Contents

Abstract

1 INTRODUCTION

2 JEANS ANALYSIS

2.1 Likelihood analysis of the Carina dSph's observed projected velocity dispersion profile

2.2 Approximate Jeans equation for MOND

2.3 Likelihood analysis of the Carina dSph's mass-to-light ratio and velocity anisotropy in MOND

3 INCORPORATING THE EXTERNAL FIELD

3.1 Approximated MOND external field prescriptions for satellite galaxies

3.2 Constant external field

3.3 Two-component simulations setup

3.3.1 The code

3.3.2 The initial conditions

3.4 Comparison of external field parametrizations

4 COMPARISON OF SIMULATED WITH OBSERVED VELOCITY DISPERSIONS

4.1 dSph stability

4.2 Sampling the mock catalogue

4.3 Fraction of good fits

5 RESULTS

5.1 Measured proper motions

5.2 Fraction of good fits as function of orbit

6 CONCLUSION

REFERENCES

APPENDIX A: PARTICLE-MESH EXTERNAL FIELD INCLUSION

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only