Stellar streams as gravitational experiments I. The case of Sagittarius

TL;DR

This paper uses N-body simulations to explore how the Sagittarius dwarf galaxy's tidal stream can test gravitational theories, specifically MOND, by matching observed structures and velocities.

Contribution

First N-body simulation of Sagittarius dwarf galaxy's tidal disruption in MOND, comparing results with observations to assess its viability as a gravitational model.

Findings

Reasonable match to observed stellar structures and velocities

Inclusion of a flattened hot gas corona improves velocity reproduction

First rigorous test of tidal dissolution in MOND at galactic scales

Abstract

Tidal streams of disrupting dwarf galaxies orbiting around their host galaxy offer a unique way to constrain the shape of galactic gravitational potentials. Such streams can be used as leaning tower gravitational experiments on galactic scales. The most well motivated modification of gravity proposed as an alternative to dark matter on galactic scales is Milgromian dynamics (MOND), and we present here the first ever N-body simulations of the dynamical evolution of the disrupting Sagittarius dwarf galaxy in this framework. Using a realistic baryonic mass model for the Milky Way, we attempt to reproduce the present-day spatial and kinematic structure of the Sagittarius dwarf and its immense tidal stream that wraps around the Milky Way. With very little freedom on the original structure of the progenitor, constrained by the total luminosity of the Sagittarius structure and by the observed stellar mass-size relation for isolated dwarf galaxies, we find reasonable agreement between our simulations and observations of this system. The observed stellar velocities in the leading arm can be reproduced if we include a massive hot gas corona around the Milky Way that is flattened in the direction of the principal plane of its satellites. This is the first time that tidal dissolution in MOND has been tested rigorously at these mass and acceleration scales.

Tables

Table 1. Table 1: Properties of the Sgr dSph remnant in terms of position, velocity, half-light radius along the minor axis r h subscript 𝑟 ℎ r_{h} , central velocity dispersion, and total luminosity. The sources are : 1 = 1 absent 1= Ibata et al. ( 1994 ) , 2 = 2 absent 2= Law & Majewski ( 2010 ) , 3 = 3 absent 3= Ibata et al. ( 1997 ) , 4 = 4 absent 4= Majewski et al. ( 2003 ) . The value of the distance and luminosity of the galaxy differs from the value of LM10 since we chose a Sun-GC distance of 8.5 kpc.

Parameter	Value	Source
RA	${18}^h{55}^m{19.5}^s$	1
Dec	$-30°{32}^{\prime }43.0"$	1
Distance	$28.5$ kpc	2
${\mu }_{\alpha }$	-2.45 mas.yr-1	2
${\mu }_{\delta }$	-1.30 mas.yr-1	2
Vrad	$+140\pm 2.0$ $\mathrm{km}{s}^{-1}$	3
$r_h$	0.6 kpc	3
${\sigma }_c$	11.4 $\mathrm{km}{s}^{-1}$	3
$L_V$	$2.4\times {10}^7$ L⊙	4

Table 2. Table 2: Initial positions and velocities of the Sgr dSph in our different simulations, to match the position and velocities of the remnant (as listed in Table  1 after 4 Gyr, X = 19.0 𝑋 19.0 X=19.0 , Y = 2.7 𝑌 2.7 Y=2.7 , Z = − 6.9 𝑍 6.9 Z=-6.9 kpc, V ​ x = 231.6 𝑉 𝑥 231.6 Vx=231.6 , V ​ y = − 40.3 𝑉 𝑦 40.3 Vy=-40.3 and V ​ z = 200.0 𝑉 𝑧 200.0 Vz=200.0 km.s -1 , where the coordinates ( X , Y , Z ) 𝑋 𝑌 𝑍 (X,Y,Z) are Galactocentric and in the right handed coordinate system.

Model	X (kpc)	Y (kpc)	Z (kpc)	Vx ($\mathrm{km}{s}^{-1}$)	Vy ($\mathrm{km}{s}^{-1}$)	Vz ($\mathrm{km}{s}^{-1}$)
N1	-53.93	-11.54	45.37	-25.74	25.30	-67.74
N2	-37.05	-24.18	51.34	-84.06	27.60	-64.41
M1	-30.71	1.00	-4.19	205.01	34.63	-138.96
M2	-16.60	11.41	-35.40	146.88	43.51	-139.14

Table 3. Table 3: Initial (t=0) and final (t=4 Gyr) stellar mass, half-light radius and central velocity dispersion of the Sgr dwarf in the two MOND simulations M1 and M2.

	Initial progenitor			Final remnant
Model	Mass (M⊙)	$r_h$ (kpc)	${\sigma }_c$ (km.s-1)	Mass (M⊙)	$r_h$ (kpc)	${\sigma }_c$ (km.s-1)
M1	$1.2\times {10}^8$	0.61	24.0	$5.1\times {10}^7$	0.64	11
M2	$1.4\times {10}^8$	0.62	25.0	$5.7\times {10}^7$	0.66	11

Equations

$${\nabla }^2\Phi =\nabla .\left[\nu \left(\frac{|\nabla {\Phi }_{\mathrm{N}}|}{a_0}\right)\nabla {\Phi }_{\mathrm{N}}\right],$$

$${\nabla }^2\Phi =4\pi G\left[{\rho }_{\mathrm{b}}+{\rho }_{\mathrm{ph}}\right],$$

$$\nu (x)=\frac{1+(1+4x^{-1}{)}^{1/2}}{2}.$$

$${\rho }_k\propto e^{\psi /\sigma }erf\left(\frac{\sqrt{\psi }}{\sigma }\right)-\sqrt{\frac{4\psi }{\pi {\sigma }^2}}\left(1+\frac{2\psi }{3{\sigma }^2}\right).$$

$${\rho }_{HG}(m)={\rho }_{0,HG}{\left(1+\frac{m}{r_{0,HG}}\right)}^{-3}\exp \left(\frac{m^2}{r_{t,HG}^2}\right),$$

Full text

11institutetext: Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de Strasbourg, 11 rue de l’Université, F-67000 Strasbourg, France

11email: guillaume.thomas@astro.unistra.fr 22institutetext: Helmholtz-Institut für Strahlen-und Kernphysik, Universität Bonn, Nussallee 14-16, D-53115 Bonn, Germany 33institutetext: Charles University in Prague, Faculty of Mathematics and Physics, Astronomical Institute, V Holešovičkách 2, CZ-180 00 Praha 8, Czech Republic

Stellar streams as gravitational experiments

I. The case of Sagittarius

Guillaume F. Thomas 11   

B. Famaey 11   

R. Ibata 11   

F. Lüghausen 22   

P. Kroupa 2233

Tidal streams of disrupting dwarf galaxies orbiting around their host galaxy offer a unique way to constrain the shape of galactic gravitational potentials. Such streams can be used as “leaning tower” gravitational experiments on galactic scales. The most well-motivated modification of gravity proposed as an alternative to dark matter on galactic scales is Milgromian dynamics (MOND), and we present here the first ever N-body simulations of the dynamical evolution of the disrupting Sagittarius dwarf galaxy in this framework. Using a realistic baryonic mass model for the Milky Way, we attempt to reproduce the present-day spatial and kinematic structure of the Sagittarius dwarf and its immense tidal stream that wraps around the Milky Way. With very little freedom on the original structure of the progenitor, constrained by the total luminosity of the Sagittarius structure and by the observed stellar mass-size relation for isolated dwarf galaxies, we find reasonable agreement between our simulations and observations of this system. The observed stellar velocities in the leading arm can be reproduced if we include a massive hot gas corona around the Milky Way that is flattened in the direction of the principal plane of its satellites. This is the first time that tidal dissolution in MOND has been tested rigorously at these mass and acceleration scales.

1 Introduction

The nature of the dark sector of the Universe represents one of the most pressing questions of modern physics. Over the years, we have built a large-scale picture in which the Universe is composed of only 5% of ordinary baryonic matter and 25% dark matter (DM), the rest being accounted for by a cosmological constant in the Einstein equation (Planck Collaboration, 2016). While arguably very succesful on large scales, this current $\Lambda$CDM picture is nevertheless plagued by a certain number of problems on small scales, especially on galaxy scales. Among those are the now famous ‘too-big-to-fail’ and ‘satellite planes’ problem (e.g., Kroupa et al., 2005, 2010; Boylan-Kolchin et al., 2011; Papastergis et al., 2015; Pawlowski et al., 2015), but also a general fine-tuning problem encapsulated in the relation between the surface density of baryons and the gravitational field in galaxies. This relation is now often referred to as the Radial Acceleration Relation (RAR McGaugh, 2016; Lelli et al., 2016b), and involves an acceleration constant $a_{0}\simeq 10^{-10}{\rm m}\,{\rm s}^{-2}$: this relation is also refelected in the tight baryonic Tully-Fisher relation (McGaugh et al., 2000; Lelli et al., 2016c; Papastergis & Shankar, 2016), in the relation between the stellar and dynamical surface densities in the central regions of galaxies (Lelli et al., 2016a; Milgrom, 2016), in the relation between the central rotation curve slope and the baryonic surface density (Lelli et al., 2013; Renaud et al., 2016a), or in the diversity of shapes of rotation curves at a given maximum velocity scale (Oman et al., 2015). While there have been attempts to explain the RAR in the classical DM picture (e.g., Navarro et al., 2016), these are far from convincing for a variety of reasons. For instance, Navarro et al. (2016) significantly overpredict the mass discrepancies, and assume, to get there, a tight stellar mass-size relation, whereas the latter is not actually observed. In fact, as shown by the comprehensive study of Desmond (2017), when using realistic correlations in the galaxy-halo connection, the RAR scatter is still overpredicted even when the abundance matching scatter is switched off (see also Wu & Kroupa, 2015, for comparisons of $\Lambda$CDM galaxy simulations with the RAR). One possible explanation for this conundrum would be that gravity is effectively different in the extremely weak field regime, and accounts for the effects usually attributed to particle DM in galaxies. This hypothesis is known as Modified Newtonian Dynamics (MOND), or Milgromian dynamics (Milgrom, 1983), which has predicted all the aforementioned galaxy scaling relations (and pushed observers to look for them), and especially the RAR, well before they were precisely assessed by observations (Famaey & McGaugh, 2012). However, such a description has, of course, also its own problems. Some faint dwarf spheroidals deviate from the predicted relation (McGaugh & Wolf, 2010) which would mean that they are out of equilibrium in this context, or that the paradigm must be extended. Some tensions also exist in globular clusters which behave in a seemingly Newtonian way when MOND would a priori predict a deviation from it (Ibata et al., 2011). Even more problematic is the need for dissipationless non-baryonic DM to reproduce the angular power spectrum of the cosmic microwave background (Planck Collaboration, 2016), as well as the failure of the MOND relation in galaxy clusters, needing either dissipationless DM or additional baryonic DM. Various hybrid frameworks have thus been proposed (e.g., Blanchet & Heisenberg, 2015; Berezhiani & Khoury, 2015) in which a new degree of freedom can in principle play the role of cosmological DM on large scales, while effectively mediating precisely a MOND force in galaxies. However, those have not yet been shown to be as succesful as $\Lambda$CDM on large scales, and they are of course less minimal. Other theories posit that a scalar field is responsible for dark energy, and it has been proposed that this scalar field may interact with dark matter. The scalar field may mediate additional long-range forces between dark matter particles of comparable strength to the canonical gravitational force, although the theory currently makes no predictions about the strength of the additional force. Of course, since the dark matter particles accelerate differently to baryons, this implies that the weak equivalence principle would have to be broken. Kesden & Kamionkowski (2006a, b) argued that disrupting satellite galaxies could allow one to investigate this effect.

Here we concentrate on MOND, a modification of gravity whose predictions are known to be quite succesful in galaxies, without dark matter (Famaey & McGaugh, 2012). However, these predictions have been mostly limited until recently to rather symmetric and static configurations (but see, e.g., Brada & Milgrom, 1999; Tiret & Combes, 2007, 2008a; Angus et al., 2014; Nipoti et al., 2007, for a few exceptions). This has changed with the advent of numerical codes, in particular the recent patches to the RAMSES code (Teyssier, 2002) developed by Lüghausen et al. (2014) and Candlish et al. (2015). The Phantom of Ramses patch (Lüghausen et al., 2014) has for instance recently led, following the seminal work of Tiret & Combes (2008b) on the topic, to the first MOND simulations of galaxy encounters with a detailed Eulerian hydrodynamical treatment of gas physics, including star formation and stellar feedback, finding that the star formation history is significantly more extended in time and space in MOND encounters than in the classical case (Renaud et al., 2016b).

Other very powerful dynamical probes of the gravitational potential, which have not yet been thoroughly investigated in MOND, are the tidal stellar streams of disrupting satellite galaxies. These are especially important probes of the three-dimensional shape of the potential outside of the Galactic plane, and at large distances. Such an analysis is not straightforward because streams do not delineate orbits and because the non-linear external field effect of MOND is likely to play a role, hence the availability of a MOND N-body code is crucial for correctly tackling the problem. This is precisely what we aim to achieve here with this series of papers on stellar streams in modified gravity.

The stream of the Sagittarius (Sgr) dwarf galaxy (Ibata et al., 1994) is the most prominent stellar structure around the Milky Way (MW) and the one for which we have data of the most exquisite precision. The orbit of the Sgr dwarf around the Galaxy is nearly polar, and the resulting tidal stream wraps a full 360∘ on the sky (Ibata et al., 2001). The detailed investigations of the stellar stream and its kinematics have led to a lot of confusion about the corresponding shape of the gravitational potential (Helmi, 2004; Johnston et al., 2005; Law et al., 2005). No model to date is satisfactory, and the benchmark model to compare with is still the one of Law & Majewski (2010, hereafter LM10) based on the spatial and kinematic structure of M-giant stars of the stream (Majewski et al., 2003, 2004), which relies on a triaxial almost oblate ellipsoid for the DM halo, but with its minor axis contained within the Galactic plane, which is not natural in $\Lambda$CDM (Debattista et al., 2013; Pearson et al., 2015). In the context of a scalar field mediating a long-range force on dark matter particles, Kesden & Kamionkowski also undertook a series of simulations of the Sagittarius dwarf galaxy in which they studied the influence of changing the strength of gravity for the dark matter. If the additional force is stronger than Newtonian gravity, the dark matter particles accelerate faster into the MW potential leaving the stars slightly behind. Because of this, any stars that do leave the system during the tidal disruption process are more likely to leave through the L2 Lagrange point, so that the resulting star stream appears asymmetric, with a less populated leading arm than the trailing arm. With observations of the stellar stream of the Sagittarius dwarf available at that time, Kesden & Kamionkowski (2006a, b) were able to rule out a 9% higher acceleration for the dark matter.

More than a decade ago, Read & Moore (2005) had shown that the orbit of the Sgr dwarf in a Milky Way MOND potential was barely distinguishable from that in a nearly-spherical to mildly oblate DM halo. However, as the stream does not follow the orbit, and as non-linearities in the MOND Poisson equation (notably the ‘external field effect’ breaking the Strong Equivalence Principle, see Section 6.3 of Famaey & McGaugh 2012) can lead to a priori unexpected effects, and it is thus urgent to revise this problem by using our modern simulation tools. Let us emphasize that a framework such as MOND has very little freedom to achieve the right final configuration of the stream, hence this exercise has a huge potential for falsification of the paradigm. Indeed, one must start from a progenitor dwarf sitting on the observed mass-size relation, and then hope that the global shape of the stream, the structure and kinematics of the remnant, and the total luminosity of the stream all fit with observations: all this with a gravitational potential of the MW fully determined by its baryonic distribution. This is not a trivial task.

We describe our method in Section 2 hereafter, then run in Section 3 two Newtonian simulations for comparison purposes, and finally run two MOND simulations of the Sgr stream in Section 4. We conclude in Section 5.

2 Method

In the following study, all N-body simulations were made with the RAMSES code (Teyssier, 2002). For simulations in the MOND framework, we used the Phantom-Of-Ramses (POR) patch developed by Lüghausen et al. (2014), who generalized the Poisson equation in the following way (Milgrom, 2010):

[TABLE]

where $\Phi$ and $\Phi_{\mathrm{N}}$ are the MOND and Newtonian potentials respectively, $\nu(x)=1$ for $x\gg 1$ (Newtonian regime) and $\nu(x)=x^{-1/2}$ for $x\ll 1$ (deep-MOND regime), and $a_{0}=1.2\times 10^{-10}$ m.s*-2* is the acceleration constant of the MOND paradigm. We can then introduce the “phantom dark matter” density $\rho_{\mathrm{ph}}$ such that

[TABLE]

which is fully defined through Eq. 1 once the baryonic distribution $\rho_{\mathrm{b}}$ (and its associated Newtonian potential $\Phi_{\mathrm{N}}$) is known, and can be seen as the MOND equivalent to the DM contribution in the classical case. This phantom DM density, that does not correspond to real particles, is computed at each time-step in POR and used to compute the MOND potential. The code is based on Adaptative-Mesh-Refinement (AMR) that increases the resolution of the grid in higher density regions without a dramatic growth in computation time. In our specific case the resolution will be highest along the tidal stream and in its progenitor. Since the energy along an orbit is conserved in a static external potential, it is important to keep the same resolution in the inner region of the satellite during the whole simulation time to avoid a numerical modification of the energy of the progenitor that could modify the orbit. In this work, we chose a minimum resolution of the AMR grid of 31 kpc and a maximum resolution of 15 pc.

In the following sections, we will proceed as follows: first we will devise benchmark Newtonian models for comparison with our MOND results, based on a spherical halo and the triaxial halo of LM10. These will be labelled simulaitions N1 and N2. These simulations are presented purely for comparisons with the following MOND simulations, and should thus not be over-interpreted in the CDM context. We will then move on to MOND simulations, with and without a massive hot gaseous corona around the Galaxy (simulations M1 and M2, respectively). Let us note here that the behavior of streams in MOND depends both on their own internal gravitational field and on the external field of the host galaxy. For these reasons, there is no general ‘simple’ test case to present here, as each situation will actually be different based on the internal properties of the progenitor and the properties of the host and of the orbit.

All our models are based on the following baryonic matter distribution for the MW (Dehnen & Binney, 1998): a double exponential stellar disk of $3.52\times 10^{10}M_{\odot}$ for the thin and thick disk components, with a scale length of 2 kpc and two scales heights of 0.3 and 1 kpc. The bulge and the interstellar medium components have respectively a mass of $0.518\times 10^{10}$ M⊙ and $1.69\times 10^{10}M_{\odot}$. This MW model is not live in the following simulations, and is represented by $5.6\times 10^{5}$ static particles of 105 M⊙ each, that generate the static potential in RAMSES. Hereafter, we will call this distribution of matter the disk model.

In all our simulations, we will follow the disruption of the Sgr dSph for 4 Gyr. The stream is indeed composed of relatively ‘young’ M-giants and is dynamically young (Majewski et al., 2003). Moreover, we will assume that the MW does not have an important modification of mass due to major mergers or heavy accretions, and that its morphology stayed the same during the last 4 Gyr. To determine the orbit of the progenitor and the initial position of the dSph in Newtonian dynamics, it is common to launch a point mass test-particle from its present location and make it run backwards, assuming a progenitor with a negligible mass compared to the mass of the host galaxy, i.e. neglecting dynamical friction with the DM halo. This is reasonable if the progenitor is less than $\sim 10^{9}M_{\odot}$, which is the case for the Sgr dSph galaxy during the last 4 Gyr in the case of a spherical DM halo (Peñarrubia et al., 2006), but which is in stark contradiction with abundance matching, which requires the mass of the dwarf galaxy to be at least $10^{10.5}M_{\odot}$. This tension could actually be problematic for CDM-based models (Kroupa 2015 but see Dierickx & Loeb, 2017). But if one ignores this tension and stays with a smaller mass, the dynamical friction is also negligible in the case of a DM halo with a similar triaxiality as that of LM10, as we will show hereafter. In MOND, this problem is of course trivially circumvented, as there is no dynamical friction outside of the MW disk (but note however, as a caveat, that within the disk, dynamical friction is actually more efficient in MOND than in Newtonian dynamics, which would thus be very important for, e.g., in-plane accretions Ciotti & Binney, 2004; Nipoti et al., 2008).

The present-day position and velocity of the Sgr dSph are listed in Table 1, where the distance and the radial velocity are heliocentric. We use the same proper motion as LM10, however we choose to adopt here the Solar peculiar motion of Schönrich et al. (2010), namely $(U_{\odot},V_{\odot},W_{\odot})=(11.1,12.24,7.25)\,{\rm km}\,{\rm s}^{-1}$ in Local Standard of Rest coordinates, and a Sun-Galactic Center (GC) distance of 8.5 kpc. Also the present-day apparent magnitude of the Sgr dSph is $m_{V}=3.63$ (Mateo et al., 1998; Majewski et al., 2003) corresponding to an absolute magnitude of $M_{V}=-13.64$, i.e. a V-band total luminosity of $2.4\times 10^{7}$ L⊙. The half-light radius along the minor axis is $r_{h}=0.6$ kpc and the central velocity dispersion is $\sigma_{c}=11.4$ ${\rm km}\,{\rm s}^{-1}$ (Majewski et al., 2003; Ibata et al., 1997), as summarized in Table 1. These quantities here refer to the remnant dSph only and do not include the associated stream. The initial conditions for our four simulations N1, N2, M1 and M2, as obtained by integrating a test-particle backwards in time for 4 Gyr, are given in Table 2.

3 Newtonian simulations with dark matter

In this section, we first run two simulations in Newtonian dynamics for testing the code and for comparison purposes, first in the case of a spherical DM halo (simulation N1), and then in the case of the LM10 triaxial DM halo (simulation N2).

In simulation N1, we add to the MW disk model a spherical DM halo based on the double-power-law model of Dehnen & Binney (1998) with $\rho_{h0}=2.46\times 10^{8}$ M⊙.kpc*-3*, $\alpha=-0.87$, $\beta=2.36$, a${}_{h}=2.66$ kpc and r${}_{h}=1000$ kpc that corresponds to M$(r<100{\rm kpc})=6.7\times 10^{11}$ M⊙ and a virial mass of M${}_{vir}=1.5\times 10^{12}$ M⊙. This halo is modelled as an external potential with $1.5\times 10^{6}$ static DM particles of individual mass 106 M⊙. For the progenitor dwarf, the initial positions and velocities are given in Table 2, and the internal structure is given by a light (DM+stars) King profile (King, 1966; Binney & Tremaine, 2008) with $M=4\times 10^{8}$ M⊙, a core radius of $r_{c}=0.65$ kpc and a ratio between the central velocity dispersion and the potential of $W=5$. After running the simulation for 4 Gyr, the resulting apparent positions on the sky, heliocentric radial velocities and heliocentric distances of the particles are shown on Fig. 1 together with the observed 2MASS M-giant stars of Majewski et al. (2004), taking their estimate of $\sim$ 20% error on the distance. In this simple N1 case, while the apparent positions and distances are reasonably well-reproduced, this is not the case of the radial velocities in the bright leading arm of the stream between RA$=140^{\circ}$ and $200^{\circ}$. The modelled velocities typically reach $-210\,{\rm km}\,{\rm s}^{-1}$ while the observations do not go below $-50\,{\rm km}\,{\rm s}^{-1}$. This means that the modelled particles fall back too fast towards the Galactic plane in this leading arm. This is actually a known problem since the work of Law et al. (2005), who showed that the problem remains for an oblate or a prolate DM halo. These observed radial velocities can for instance also be compared to the recent model of Dierickx & Loeb (2017): as can be seen on their figure 10, the problem is at least as severe in terms of radial velocities, whilst apparent positions on the sky are worse. One possibility is to discard these data as non-members of the stream. Another is to consider an alternative halo, which is precisely what LM10 did with their proposed triaxial halo. Here we reproduce their results with simulation N2 in Fig. 1. Following e.g., Peñarrubia et al. (2010), we modelled this halo with a NFW profile with virial mass M${}_{vir}=1.1\times 10^{12}$ M⊙, substituting the spherical radius with an elliptical radius $m$ where $m^{2}=x^{2}/a^{2}+y^{2}/b^{2}+z^{2}/c^{2}$ and $a=0.44,\,b=1.0,\,c=1.0$, similar to LM10. The progenitor dwarf here is modelled with a King profile with a mass $M_{init}=6.8\times 10^{8}$ M⊙, a core radius $r_{c}=0.65$ kpc and a ratio between the central velocity dispersion and the potential $W=4$. Taking into account the loss of mass over 4 Gyr, we also use Chandrasekhar’s formula (Chandrasekhar, 1943) to estimate the effect of dynamical friction on this model: the difference is shown in Fig. 2 where it is clear that for a light Sgr dwarf, the dynamical friction does not modify the orbit significantly. With such a triaxial halo, the kinematics of stars in the bright leading arm at $140^{\circ}<{\rm RA}<200^{\circ}$ is reproduced, but note that the projected dispersion on the sky is clearly overproduced and that the stream is less extended in distance (Fig. 1, bottom panel). Also, Debattista et al. (2013) have shown that such a halo is not able to host a stable disk, and recently, Pearson et al. (2015) have shown that due to this triaxiality, the stars of the Palomar 5 (Pal 5, see e.g. Thomas et al., 2016) stream would be on chaotic orbits and create a fanning shape at the end of the stream, which is not observed: they show that a spherical DM halo is better at reproducing the thin and coherent observed structure of the Pal 5 stream. We will return to the topic of modelling Pal 5 in modified gravity in a subsequent contribution.

4 MOND simulations

Contrary to the above Newtonian dynamics case, in which there is significant freedom to choose the potential (because it is dominated at the distance of the Sgr dSph by the DM halo), the baryonic mass distribution of the host galaxy has a crucial importance in MOND, as it fully fixes the gravitational potential. The same is actually true for the progenitor dwarf galaxy: it must sit on the observed stellar mass-size relation, and at the end of the simulation, the remnant should resemble what is observed: again, there is very little freedom here, as the gravitational potential of the progenitor is entirely determined by its baryonic content. Finally, there is an additional effect, which is unique to MOND, which we would like to investigate here: the ‘external field effect’. This effect is absolutely unique to theories like MOND which break the strong equivalence principle, and are distinct from the usual tidal effects. It means that the internal dynamics of a satellite system does not decouple from the external field produced by its host system, drastically reducing the amount of ‘phantom dark matter’ at pericenter compared to the apocenter or isolated case, and even producing pockets of negative phantom DM densities at places. It is this effect that led to the successful prediction of the small velocity dispersion of the dwarf galaxy Crater II in MOND (McGaugh, 2016; Caldwell et al., 2016). One of the things we would like to understand is whether this effect plays a role in shaping stellar streams in MOND, and whether this provides a distinct signature from Newtonian gravity.

For our study, we need to choose a transition function $\nu$ in Eq. 1, between the MONDian and the Newtonian regime around the $a_{0}$ acceleration scale. The shape of this transition will not be absolutely crucial in our study, since the MW potential is already in the deep-MOND regime at the distance of the Sgr dSph. Following Famaey & Binney (2005) and Zhao & Famaey (2006), for galaxies111See however Hees et al. (2016) for tight constraints in the Solar System for the strong gravitation regime, we choose an interpolation function of the form:

[TABLE]

To estimate the external field effect, we first computed analytically the density of phantom dark matter around a Plummer sphere with Plummer radius of $r_{s}=0.85$ and mass of $5.1\times 10^{7}$ M⊙, corresponding roughly to the current baryonic mass of the Sgr remnant assuming a mass-to-light ratio of $\sim 2.1$, on a circular orbit around a point mass of $5.6\times 10^{10}$ M⊙ at a distance of 20 kpc and 80 kpc, hence roughy corresponding to the expected pericenter and apocenter of the Sgr dwarf. As shown in Fig. 3 the distribution of phantom-dark-matter around the Plummer sphere is not spherical due to the external field effect (Wu et al., 2010). There are even pockets of negative density which compress the dwarf. At apocenter, the phantom-dark-matter mass is slightly larger, and so is its negative phantom density counterpart, but this happens beyond $\sim 2$ kpc from the center of the dwarf, and thus it does not affect the stream. The disruption of the satellite appears to be influenced mostly by the shape of the potential of the progenitor at pericenter, which is not lopsided and is very similar to the Newtonian case with dark matter, as can be seen on the top-right panel of Fig. 3.

After these analytical preliminaries, we will now model the progenitor stellar distribution with a King profile, which is closer to the shape of the observed stellar distribution of dSph galaxies, and produce full N-body simulations of the Sgr disruption in MOND.

4.1 Modellling the Sgr dSph in MOND

For the initial conditions on the internal kinematics of the dwarf, we constructed a fully self-consistent MONDian King model of a dwarf galaxy sitting on the observed stellar mass-size relation (Dabringhausen & Kroupa, 2013). It will be represented in our simulation by $10^{5}$ N-body particles. The King model (see Section 4.3 of Binney & Tremaine, 2008, for details) is defined by a distribution function depending on energy, which once integrated over velocity space, gives a density $\rho_{k}$ proportional to the relative binding potential $\psi$:

[TABLE]

To construct a MOND King model, all we need to do is replace the Newtonian relative potential $\psi_{N}$ of the classical formulation by the MONDian relative potential $\psi$ in that equation, and we integrate from inside out, assuming that the central region is in the Newtonian regime. We halt the inside-out integration once outside of the desired radial range, and check with POR that the King model is indeed in equilibrium in isolation in MOND.

4.2 The Sgr stream in the disk model

In this section, we present a benchmark simulation of the Sgr stream in MOND, denoted M1. In this simulation, the only source of gravity is the (non-live) MW model presented in Sect. 2 and the live fully baryonic self-consistent King model devised with the method outlined above.

Our first constraint will be to reproduce the total luminosity of the Sgr structure (the Sgr stream + the Sgr dSph) , which is of the order of (or a bit less than) $L\sim 10^{8}$ L⊙ (Niederste-Ostholt et al., 2010). After a few tries, we choose an initial stellar mass of $1.2\times 10^{8}M_{\odot}$ for the model M1. We then choose a core radius such that the half-light radius matches the observed stellar mass-size relation of other dwarf galaxies (Dabringhausen & Kroupa, 2013). For this, we choose a core radius of $r_{c}=0.6$ kpc and $W=5$ for the ratio between the central velocity dispersion and potential. This leads to a central velocity dispersion of $\sigma_{c}=24$ km.s*-1* and a half-light radius of $r_{h}=0.61$ kpc (see Table 3). As we know that dynamical friction of the dwarf galaxy with an inexistent DM halo does not take place in MOND, we can safely integrate the orbit backwards in time in the MW disk model MOND potential to get the initial positions and velocities listed in Table 2. We then run our M1 N-body simulation forwards in time for 4 Gyr.

Let us insist here on the very little wiggle room we had to choose these initial conditions. The MW potential is entirely determined by its baryonic distribution and so are the initial positions and velocities. The Sgr initial model is only made of stars and should fall on the stellar-mass size relation for similar isolated dwarfs. It is then truly remarkable that the morphology of the remnant after 4 Gyr of disruption is in perfect agreement with that observed by Majewski et al. (2003) with the 2MASS survey. In Fig. 4 we compare the density of stars along the minor axis in our M1 simulation (in blue) and in the observed M-giants (in red). The projected minor axis is least affected by tidal effects, and is thus the best direction to evaluate the morphology of the observed remnant. In the M1 simulation, our remnant has a final mass $M_{\rm final}=5.1\times 10^{7}$, a position angle of $104\degr$, a half-light radius along the minor axis of $r_{h}=0.64$ kpc and a central velocity dispersion of $\sigma_{c}=11$ ${\rm km}\,{\rm s}^{-1}$ (see Table 3), which is in excellent agreement with observations, (see Table 1) assuming a mass-to-light ratio $\gamma_{*}=2.1$ in the V-band. In Fig. 5 we also show the surface brightness of the remnant for the M1 simulation on the upper panel. This morphology can be directly compared to Figure 4 of Majewski et al. (2003).

In Fig. 6, we show a movie of the M1 simulation, and also provide a projection of the final shape of the stream at present time in the Sgr orbital plane defined by Belokurov et al. (2014) (note that the x-axis is points in the opposite direction to that in their paper). We note that the modelled stream extends to slightly larger distances than 100 kpc, but not much beyond that. This could be different in a simulation where the MW model itself would be live and responding to the gravitational pull of the dwarf. We also display the resulting projected positions on the sky, radial velocities and distances of our M1 simulation in Fig. 7, together with the Newtonian simulation of Law & Majewski (2010). The positions and distances of the M1 stream are in reasonable agreement with the observations. But again, the radial velocities do not match well in the leading arm, exactly as in simulation N1 and in all Newtonian models with spherical, oblate, or prolate halos (e.g. Law et al., 2005; Dierickx & Loeb, 2017). While this could be due to the influence of other satellites such as the LMC, we consider hereafter another possible solution based on the influence of the hot gas corona around the Galaxy. We note that the M1 stream is actually very similar to that produced in the N1 simulation (see Fig. 1, the small differences are of the same scale as the resolution of the grid of the RAMSES code), except that it is more self-consistent as the progenitor consists only of stars, and is more constrained in the sense that the progenitor had to obey observed scaling relations. This similitude with the N1 stream means that the external field effect does not have an important effect on the morphology of the stream in MOND, and is dominated by the effective gravitational potential of the MW at the distances probed by the orbit of Sgr.

4.3 The Sgr stream in the hot corona model

Here we present our second MOND simulation, denoted M2, in which a massive diffuse hot gas (HG) corona is included. The presence of such a hot diffuse gaseous corona around the MW at a temperature of $\sim 10^{6}$ K has been proposed for many years as a significant reservoir of baryonic matter that can be traced by the O VII and O VIII emission and absorption lines (Paerels & Kahn, 2003) in the soft X-ray band. Recent measurements obtained with the XMM-Newton and Suzaku X-ray telescopes have estimated a mass for this hot gaseous spherical component between $(0.5\,$-$\,1)\times 10^{11}$ M⊙, and could even be as massive as $1.5\times 10^{11}$ M⊙ (but see Gatto et al., 2013; Salem et al., 2015, for lower estimates), while it should extend up to at least 100 kpc (Gupta et al., 2012; Fang et al., 2013). Such a corona may actually be the remnant of the formation of the vast polar structure of satellite galaxies (VPOS) if it consists of dwarf galaxies and star clusters that formed within a large gas-rich tidal arm which was pulled out about 10 Gyr ago of either the young MW or the passing other galaxy which may have been Andromeda (Pawlowski et al., 2011; Zhao et al., 2013; Hammer et al., 2013). The star-formation efficiency is at most a few per cent on the scales of molecular clouds and within a tidal arm of extend of 100 kpc or more it is likely to have been significantly smaller, probably less than 0.1 per cent. The stellar mass in the VPOS comprises about $10^{9}\,M_{\odot}$ in total, such that this very rough estimate implies about $10^{11}\,M_{\odot}$ to have been the mass of the tidal arm(s) within which formed the present-day VPOS constituents. Today the gas is likely to be oriented in a thick oblate structure aligned with the VPOS, partially being derived also from ram-pressure stripping and gas blown out from the young tidal dwarf galaxies which were the precursors of the present-day satellite galaxies. Note that the LMC alone could be the progenitor of this corona, as small disk galaxies are very gas rich. The corona is likely to be hot as its low density implies a long cooling time and possibly heating agents such as the motions of the satellite galaxies and intergalactic radiation may continue to heat this oblate corona. The existence of such an ancient structure which would be the remnant of the tidal arm(s) is rather speculative though, and theoretical research is needed to provide constraints on the stability and existence of such a structure, if it can be created from a tidal arm at all. It is therefore rather interesting that observational evidence has appeared which suggests the existence of a hot corona, and our M2 model here shows that such a structure significantly improves the reproduction of the Sgr orbit and its tidal arms.

In our simulation we modelled this hot gaseous component with a triaxial cored halo profile:

[TABLE]

where $r_{0,HG}$ is the core radius, $r_{t,HG}$ the truncation radius, and the oblateness is defined throug the elliptical radius $m$ such that $m^{2}=x^{2}/a^{2}+y^{2}/b^{2}+z^{2}/c^{2}$ and $a=0.44,\,b=1.0,\,c=1.0$. We restrict the total mass of the flattened corona within 100 kpc to be M$(<\,100\,kpc)\approx 1.5\times 10^{11}M_{\odot}$, which parametrizes this profile with a total mass up to $2.6\times 10^{11}$ M⊙ within the truncation radius, but the mass beyond 100 kpc is not necessary, as cutting off the mass beyond 100 kpc does not influence the stream’s formation and kinematics. We also show in Fig. 8 how little this additional corona affects the MW rotation curve.

Since the total baryonic mass of the MW is higher in this M2 model than in M1, the tidal effects on the Sgr dwarf are stronger, and one needs to increase the mass of the progenitor to keep a realistic remnant after 4 Gyr. This can still be done while remaining consistent with the observed total luminosity of the Sgr structure, by taking an initial mass of $M_{init}=1.4\times 10^{8}{\rm\,M_{\odot}}$ , an initial core radius of $r_{c}=0.6$ kpc and $W=5$, corresponding to a central velocity dispersion of $\sigma_{c}=25$ km.s*-1* and a half-light radius $r_{h}=0.62$ kpc, which is still in agreement with the mass-size relation of other dwarf galaxies (Dabringhausen & Kroupa, 2013). We neglect the dynamical friction due to the gas particles of the HG corona, precisely because these are hot and not very reactive to the perturbation from the dwarf, hence not prone to creating dynamical friction, and also because our Newtonian study has shown that even a massive triaxial DM halo does not affect much a light progenitor in Newtonian dynamics, so we can expect the same for a much less massive corona in MOND.

After 4 Gyr of disruption, the remnant has a similar morphology to that in the case of the M1 disk model, as can be seen on Fig. 5 where we show the surface brightness of the remnant for the M1 disk model on the top panel and for the M2 HG model on the bottom panel. The morphology in the two cases is very close to the observed morphology of the Sgr dSph in the 2MASS survey (see e.g. Figure 4 of Majewski et al., 2003). The M2 remnant has a final mass of $M_{\rm final}=5.7\times 10^{7}{\rm\,M_{\odot}}$, a position angle of $104\degr$, a half-light radius along the minor axis of $r_{h}=0.66$ kpc and a central velocity dispersion of $\sigma=11$ ${\rm km}\,{\rm s}^{-1}$ (see Table 3), that correspond to a stellar mass-to-light ratio of $\gamma_{*}=2.4$ to reproduce the observed luminosity of the remnant.

The fact that the remnant is realistic in both the M1 and M2 simulations is interesting, as it results from very different orbital and mass loss histories. The difference between the orbit of the Sgr dSph in the two MONDian models M1 and M2 is shown in Fig 9 where the orbit for the M1 model is in red and the orbit for the M2 model in green. The pericenter in both cases is of about $\sim 15$ kpc but the apocenter is much closer in the case of the HG corona, $55$ kpc instead $80$ kpc for the M1 model. This is a consequence of the more massive baryonic model of the MW, and could be problematic to reproduce tentative pieces of the stream detected at distances of the order of 100 kpc. It is interesting to see that in the case of the M2 simulation, the pericenter after 2.4 Gyr of disruption is as large as $30$ kpc and thus limits the mass-loss of the progenitor for this orbit, as shown in Fig. 10. The bound fraction of particles can also increase with time, which is especially the case for the M1 model. In both the M1 and M2 models, all these effects conspire to give a realistic remnant. Note how the mass of the remnant is predicted to increase again in the future within the M1 simulation.

Finally, Fig. 11 represents the observed projected quantities in the same way as Fig. 1. In this M2 model, the simulated radial velocities reproduce very well the observations of the M-giant stars of Majewski et al. (2004) in the leading arm. Contrary to the N2 simulation, the projected dispersion along the stream is also more consistent with the observations, making this model quite superior to that resulting from the LM10 triaxial halo. The much smaller mass of the HG corona, especially in the inner parts, is not likely to destabilize the disk as much as the LM10 DM halo, and could perhaps also be consistent with the Pal 5 stream, which will be the subject of another contribution. However, note that the heliocentric distance extent of the stream in this M2 model is much less than in the M1 model, and it does not fit the BHB stars of Belokurov et al. (2014), but this is also a shortcoming of the current best simulation in Newtonian gravity to reproduce the leading arm velocities, made by Law & Majewski (2010).

In Milgromian dynamics, the VPOS would likely have formed as a tidal interaction early in the history of our Galaxy, and could have associated with it a much larger component in hot gas than the present-day stellar mass in the VPOS. If this hot gaseous corona shares similar flattening and orientation as the VPOS, our model M2 shows that it would affect the orbit of Sgr and its tidal arms in the right way to make the model consistent with the observed radial velocities of the M-giants near ${\rm RA}=180\deg$. However, we note that recent estimates of the mass of the corona are significantly lower than the mass we assumed here (Gatto et al., 2013; Salem et al., 2015). Therefore it is possible that a more realistic model has properties in between M1 and M2.

5 Discussion and conclusions

In this paper, we have for the first time rigorously tested the MOND paradigm in a regime where it had never been tested before, namely using stellar streams as a gravitational experiment. We started a series of papers on this topic with an analysis of the most prominent stream of the MW, the Sgr stream. Reproducing such a stream is not a trivial task for a theory such as MOND, as it requires the global shape of the gravitational potential, fully determined by the baryon distribution, to conspire with the tidal effects on the disrupting dwarf galaxy originally (i) sitting on the observed stellar mass-size relation and (ii) reproducing the total luminosity of the stream, to produce both (I) the correct shape and kinematics of the stream, and (II) the correct internal structure and kinematics of the remnant dwarf spheroidal. There was thus very little wiggle room for our simulations to reproduce the observations.

Starting with a King model of total mass between $1.2\times 10^{8}M_{\odot}$ and $1.4\times 10^{8}M_{\odot}$ and half-light radii between 610 pc and 620 pc, in accordance with the observed stellar mass-size relation, our two MOND simulations M1 and M2, with and without a flattened hot gas corona, both gave a very realistic remnant and a quite realistic stream morphology, despite very different orbital and mass-loss histories.

Our M1 model is our fiducial model without hot gas, for which we provided a movie in Fig. 6. This M1 simulation predicts stream debris reaching out to distances slightly larger than 100 kpc, produces a remnant that matches well the observations of the Sgr dSph at the present epoch, and is as well a very good match to the positions of the bright stream arms on the sky. This is quite an achievement for a model with so little freedom. However, the M1 model does not seem to reproduce well the observed radial velocities of M giants in the leading arm, between RA$=140^{\circ}$ and $200^{\circ}$, a well-known problem in Newtonian dynamics for all DM halo models with a spherical, oblate, or prolate shape (e.g. Law et al., 2005; Dierickx & Loeb, 2017).

The similarity of the M1 stream with our Newtonian spherical DM halo model N1 means that the external field effect does not have an important influence on the morphology of the Sgr stream in MOND, and is dominated by the quasi-spherical ‘phantom DM’ halo of the MW at the distances probed by the orbit of Sgr. This will not necessarily be the case for less massive progenitors (such as Palomar 5), which will be the topic of a further paper in this series.

If the radial velocities in the leading arm are not considered as misidentified stream members, two possible solutions in MOND would be (i) the influence of other satellites such as the LMC (e.g., Vera-Ciro & Helmi, 2013; Laporte et al., 2016) or (ii) the influence of a flattened hot gas corona aligned with the VPOS. While solution (i) will be the topic of further study, we examined here the plausibility of the solution (ii) in our M2 simulation. In this M2 model with a massive flattened hot gas corona, the simulated radial velocities reproduce very well the observations of the M-giant stars of Majewski et al. (2004) in the leading arm. As opposed to the N2 simulation based on the triaxial DM halo of LM10, the projected dispersion along the stream is also more consistent with the observations in this M2 simulation. Nevertheless the extent of the stream is predicted to be much smaller in this M2 model than in the fiducial M1 case, which could be problematic (checking whether BHB stars at large distances could be fitted in such a model would need to consider a much longer orbit), and the assumed mass of the corona is high compared to some recent estimates. In the future, it should be interesting to test if other configurations of the progenitor or of the flattened hot gaseous corona are also able to reproduce the bifurcation seen in the Northern sky in SDSS and recently extended in the Southern hemisphere by Navarrete et al. (2016). In the same way, modelling other streams, especially Pal 5, with this hot gaseous corona in MOND should be done in the future, to see if the Pal 5 stream stays coherent in this model, contrary to its structure in the MW potential derived by LM10 (Pearson et al., 2015). It will also be interesting to see if the external field effect plays a more important role for lower mass progenitors such as Pal 5, and could leave a distinctive MOND signature in the stream. All this will be the topic of further papers in this series on using stellar streams as gravitational laboratories.

Acknowledgments

The authors would like to thank Stacy McGaugh and Nicolas Martin for very helpful discussions and precious remarks.