TY - JOUR
AU - Zhao,, Hongsheng
AB - Abstract We attempt to understand the planes of satellite galaxies orbiting the Milky Way (MW) and M31 in the context of Modified Newtonian Dynamics, which implies a close MW–M31 flyby occurred ≈8 Gyr ago. Using the timing argument, we obtain MW–M31 trajectories consistent with cosmological initial conditions and present observations. We adjust the present M31 proper motion within its uncertainty in order to simulate a range of orbital geometries and closest approach distances. Treating the MW and M31 as point masses, we follow the trajectories of surrounding test particle discs, thereby mapping out the tidal debris distribution. Around each galaxy, the resulting tidal debris tends to cluster around a particular orbital pole. We find some models in which these preferred spin vectors align fairly well with those of the corresponding observed satellite planes. The radial distributions of material in the simulated satellite planes are similar to what we observe. Around the MW, our best-fitting model yields a significant fraction (0.22) of counter-rotating material, perhaps explaining why Sculptor counter-rotates within the MW satellite plane. In contrast, our model yields no counter-rotating material around M31. This is testable with proper motions of M31 satellites. In our best model, the MW disc is thickened by the flyby 7.65 Gyr ago to a root mean square height of 0.75 kpc. This is similar to the observed age and thickness of the Galactic thick disc. Thus, the MW thick disc may have formed together with the MW and M31 satellite planes during a past MW–M31 flyby. methods: data analysis, methods: numerical, Galaxy: disc, Galaxy: kinematics and dynamics, galaxies: groups: individual: Local Group, dark matter 1 INTRODUCTION The standard cosmological paradigm [Λ cold dark matter (ΛCDM); Ostriker & Steinhardt 1995] faces several challenges in the relatively well-observed Local Group (e.g. Kroupa 2015 and references therein). In particular, the satellite systems of its two major galaxies – the Milky Way (MW) and Andromeda (M31) – are both highly flattened. For the MW, this has been suspected for several decades (Lynden-Bell 1976, 1982). In light of subsequently discovered satellites, this led Kroupa, Theis & Boily (2005) to argue that its satellite system is inconsistent with ΛCDM as it is too anisotropic. Later proper motion measurements showed that most of its satellites corotate within a well-defined plane (Metz, Kroupa & Libeskind 2008; Pawlowski & Kroupa 2013). Recently discovered ultra-faint satellites, globular clusters and tidal streams independently prefer a similarly oriented plane (Pawlowski, Pflamm-Altenburg & Kroupa 2012; Pawlowski & Kroupa 2014). Although some flattening is expected in ΛCDM (e.g. Butsky et al. 2016), it remains difficult to explain the very small thickness of the MW satellite system and its coherent rotation (Ahmed, Brooks & Christensen 2017). An analogous situation was suspected around M31 (Metz, Kroupa & Jerjen 2007, 2009) and recently confirmed by Ibata et al. (2013) using the homogeneous sky coverage of the Pan-Andromeda Archaeological Survey (McConnachie et al. 2009). Despite its greater distance, the detection of the M31 satellite plane is rather secure because we presently view it almost edge-on. It probably rotates coherently given the observed redshift gradient across it (Ibata et al. 2013), though proper motions are needed to confirm this. As with the MW, it is difficult to understand the observed properties of the M31 satellite system in ΛCDM (Ibata et al. 2014). The solution to these puzzling observations should also consider the two planes of non-satellite LG galaxies found by Pawlowski, Kroupa & Jerjen (2013). Just outside the LG, a satellite plane has recently been identified around Centaurus (Cen) A (Müller et al. 2018). This structure is likely corotating because its members exhibit a radial velocity trend similar to M31. The satellites of M81 also appear to lie in a flattened distribution, though their kinematics are not currently known (Chiboucas et al. 2013). These discoveries may be related to a vast plane of dwarf galaxies recently found near M101 (Müller et al. 2017). This is not a satellite galaxy plane as it extends over 3 Mpc. Such a structure appears difficult to find in cosmological simulations of the ΛCDM paradigm (González & Padilla 2010, fig. 8). Although ΛCDM is widely considered to work well on large scales and at early times (e.g. Planck Collaboration XIII 2016), it faces some tension regarding the assumptions of homogeneity and isotropy (Javanmardi et al. 2015; Kroupa 2015; Javanmardi & Kroupa 2017; Bengaly et al. 2018). The model also has difficulty in explaining the strength of the 21 cm absorption feature corresponding to the hyperfine transition of cold neutral hydrogen at redshift z ≈ 17 (Bowman et al. 2018). While those authors attributed their unexpectedly strong detection to baryon–dark matter interactions, only a narrow range of parameters remains viable once laboratory and other constraints are considered (Berlin et al. 2018).1 Ewall-Wice et al. (2018) considered a somewhat more conventional explanation based on a radio background from accreting black holes that formed out of stars or via direct collapse of gas clouds. They found that this scenario would have to be highly contrived to match other constraints, in particular the optical depth to reionization (Heinrich & Hu 2018). For example, the black holes required by their model would need to stop growing at z ≈ 16. Given these difficulties and the many a priori predictions of Modified Newtonian Dynamics (MOND) which were subsequently confirmed (e.g. Famaey & McGaugh 2012), it is worth considering whether the matter content of the Universe might be purely baryonic (McGaugh 1999). Indeed, this seems to be one of the few viable scenarios for explaining the strong observed 21 cm absorption feature without significant fine tuning (McGaugh 2018). Such a purely baryonic universe could have an expansion history very similar to ΛCDM (Zhao 2007). Although high-redshift observations like this are interesting, the non-linear dynamics of gravitational systems are much easier to observe in the LG. Here, ΛCDM faces another recently discovered problem concerning the very high radial velocities of non-satellite LG dwarf galaxies at distances of 1−3 Mpc (Pawlowski & McGaugh 2014). We investigated this problem based on a timing argument analysis of the LG (Kahn & Woltjer 1959; Einasto & Lynden-Bell 1982) extended to include test particles representing LG dwarfs. Following on from previous spherically symmetric dynamical models (Sandage 1986; Peñarrubia et al. 2014), we constructed an axisymmetric model of the LG consistent with the almost radial MW–M31 orbit (van der Marel et al. 2012) and the close alignment of Cen A with this line (Ma et al. 1998). Treating LG dwarfs as test particles in the gravitational field of these three massive moving objects, we investigated a wide range of model parameters using a full-grid search (Banik & Zhao 2016). None of the models provided a good fit, even when we made reasonable allowance for inaccuracies in our model as a representation of ΛCDM based on the scatter about the Hubble flow in detailed N-body simulations of it (Aragon-Calvo, Silk & Szalay 2011). This is because several LG dwarfs have Galactocentric Radial Velocities (GRVs) much higher than expected in our best-fitting model, though the opposite is rarely the case (Banik & Zhao 2016, fig. 9). We found that this high-velocity galaxy (HVG) problem should remain even when certain factors beyond our model are included, in particular the Large Magellanic Cloud (LMC) and the Great Attractor (GA). To test whether the issue might be resolved using a three-dimensional (3D) model of the LG, we borrowed an algorithm described in Phelps, Nusser & Desjacques (2013). The typical mismatch between observed and predicted GRVs in the best-fitting model is actually slightly higher than in the two-dimensional (2D) case, with a clear tendency persisting for faster outward motion than expected (Banik & Zhao 2017, figs 7 and 9). These results are similar to those obtained by Peebles (2017) using a similar algorithm. Despite using a very different method to Banik & Zhao (2016), the conclusions remain broadly similar. We recently revisited this analysis by performing a much more thorough search for the best-fitting model (Banik & Zhao 2018b, section 4.1). As our best-fitting model still yielded several HVGs that it could not account for, we considered possible ΛCDM-based explanations for them at some length (see their section 7.2). Not coming up with any plausible solutions, we looked for clues in the spatial distribution of the HVGs. This revealed that they lie very close to a plane extending over  ∼3 Mpc, with an isotropic distribution ruled out at 3σ (see their table 7). ΛCDM has difficulty explaining the kinematics of even one HVG, let alone the 6 we identified. We therefore considered how these observations might be explained in the alternative theory called MOND (Milgrom 1983). In MOND, galaxies are not surrounded by haloes of dark matter. Instead, their dynamical effect is provided by an acceleration-dependent modification to gravity. The gravitational field strength g at distance r from an isolated point mass M transitions from the usual |$\frac{GM}{r^2}$| law at short range to \begin{eqnarray*} g = \frac{\sqrt{GMa_{0}}}{r} \ \ \ \text{for } \ r \gg \sqrt{\frac{GM}{a_{0}}}. \end{eqnarray*} (1) MOND introduces Milgrom’s constant a0 as a fundamental acceleration scale of nature below which the deviation from Newtonian dynamics becomes significant. Empirically, a0 ≈ 1.2 × 10−10 m s–2 to match galaxy rotation curves (McGaugh 2011). Remarkably, this is similar to the acceleration at which the classical energy density in a gravitational field (Peters 1981, equation 9) becomes comparable to the dark energy density uΛ = ρΛc2 implied by the accelerating expansion of the Universe (Riess et al. 1998). Thus, \begin{eqnarray*} \frac{g^2}{8\rm {\pi }G}\,\,>\,\, \ u_{\Lambda } \ \ \Leftrightarrow \ \ g \ \lesssim \ 2\rm {\pi }a_{0}. \end{eqnarray*} (2) This suggests that MOND may arise from quantum gravity effects (e.g. Milgrom 1999; Pazy 2013; Verlinde 2016; Smolin 2017). Regardless of its underlying microphysical explanation, it can accurately match the rotation curves of a wide variety of both spiral and elliptical galaxies across a vast range in mass, surface brightness and gas fraction (Lelli et al. 2017 and references therein). It is worth emphasizing that MOND does all this based solely on the distribution of luminous matter. Given that most of these rotation curves were obtained in the decades after the MOND field equation was first published (Bekenstein & Milgrom 1984), it is clear that these achievements are successful a priori predictions. These predictions work due to underlying regularities in galaxy rotation curves that are difficult to reconcile with the collisionless dark matter haloes of the ΛCDM paradigm (Desmond 2017a,b; Salucci & Turini 2017). In the LG, MOND implies a much stronger MW–M31 mutual attraction than does ΛCDM. Acting on their almost radial orbit (van der Marel et al. 2012), this leads to a close encounter 9 ± 2 Gyr ago (Zhao et al. 2013). Due to the high MW–M31 relative velocity around the time of their flyby, they would likely have gravitationally slingshot several LG dwarfs out at high speed. The dwarfs flung out almost parallel to the motion of the perturber would have reached the fastest speeds. We used a MOND model of the LG to demonstrate this effect, showing that the HVGs ought to define the MW–M31 orbital plane (Banik & Zhao 2018b, section 3). Observationally, the HVGs do define a rather thin plane, with the MW–M31 line only 16° out of this plane (see their table 4). Thus, we argued that the HVGs may preserve evidence of a past close MW–M31 flyby and their fast relative motion at that time. During this flyby, a tidal tail would likely have formed and might later have condensed into satellite galaxies of the MW and M31. This phenomenon occurs in some observed galactic interactions (Mirabel, Dottori & Lutz 1992) and in MOND simulations of them (Tiret & Combes 2008; Renaud, Famaey & Kroupa 2016). Due to the way in which such tidal dwarf galaxies form out of a thin tidal tail, they would end up lying close to a plane and corotating within that plane (Wetzstein, Naab & Burkert 2007), though a small fraction might well end up counter-rotating depending on the exact details (Pawlowski, Kroupa & de Boer 2011). This might be how the MW and M31 satellite planes formed (Zhao et al. 2013). In principle, the LG satellite planes could form in some other galactic interaction(s) than a MW–M31 flyby (e.g. Hammer et al. 2013), an infeasible scenario in ΛCDM as dynamical friction between their overlapping dark matter haloes would cause a rapid merger (Privon et al. 2013). However, any second-generation origin of the MW and M31 satellite planes runs into the issue that such satellites would be free of dark matter (Barnes & Hernquist 1992; Wetzstein et al. 2007). This is due to the dissipationless nature of dark matter and its initial distribution in a dispersion-supported near-spherical halo. During a tidal interaction, dark matter of this form is clearly incapable of forming into a thin dense tidal tail out of which dwarf galaxies might condense. As a result, any tidal dwarf galaxies would be purely baryonic and thus have only a very small escape velocity, preventing them from subsequently accreting dark matter out of the MW halo. For this reason, ΛCDM struggles to explain the high observed internal velocity dispersions of the MW satellites coherently rotating in a thin plane (McGaugh & Wolf 2010). That work paid careful attention to the issue of whether tides were responsible for the high velocity dispersions, eventually ruling out this explanation in a ΛCDM context (see their fig. 6). A similar situation is evident around M31 – its satellite galaxies also have rather high internal velocity dispersions (McGaugh & Milgrom 2013) but mostly lie in a thin plane (Ibata et al. 2013). The velocity dispersions could be explained if the M31 satellites formed primordially and thus retained their dark matter haloes. However, this is difficult to reconcile with the anisotropy of the M31 satellite system as a whole (Ibata et al. 2014), which is suggestive of a tidal origin. After careful consideration of several proposed ΛCDM scenarios for how primordial satellites might end up in a thin plane, Pawlowski et al. (2014) concluded that none of them agreed with observations for either the MW or M31. It was later shown that baryonic effects are unlikely to provide the necessary anisotropy if one sticks to a primordial origin for the satellites (Pawlowski et al. 2015b). This issue was revisited by performing a high-resolution ΛCDM hydrodynamical simulation of a MW analogue in a cosmological context (Maji et al. 2017). Although this unpublished article claimed that the results were consistent with observations, it has recently been shown that this is not the case (Pawlowski et al. 2017). Those authors showed that the satellite galaxy distribution of Maji et al. (2017) was consistent with isotropy. However, the actual MW satellite system is inconsistent with isotropy at more than 5σ once the survey footprint is taken into account (Pawlowski 2016). More recent hydrodynamical ΛCDM simulations also fail to yield highly flattened satellite systems like those observed in the LG (Ahmed et al. 2017). The mild amount of flattening in these simulations might not even be related to baryonic effects as similar results arise in dark matter-only simulations (Garaldi et al. 2018). In any case, it is difficult to see how baryonic effects can help to form a 200 kpc-wideplane of primordial satellites (Shao et al. 2018). On the other hand, tidally formed satellites would lack dark matter and thus have only a very small internal velocity dispersion in disagreement with observations. For a recent review of the satellite plane problem, we refer the reader to Pawlowski (2018). This considers both LG satellite planes and the recently discovered one around Cen A (Müller et al. 2018). Given how difficult it is for primordial satellites to explain the observed satellite planes, it is important to note that tidal dwarf galaxies ought to form in any viable cosmological model. In ΛCDM, these would have rather different properties to primordial dwarfs due to the former lacking significant quantities of dark matter (Kroupa 2012, section 4). However, observational evidence exists only for the latter (Dabringhausen et al. 2016). In MOND, this is to be expected because no galaxies have dark matter haloes, whether formed primordially or as tidal dwarfs. MOND predicts that there will be a one-to-one correspondence between the baryonic distribution of a galaxy and its internal dynamics regardless of how it formed. In the LG, this modification to gravity raises the expected internal velocity dispersions of purely baryonic MW and M31 satellites enough to match observations (McGaugh & Wolf 2010; McGaugh & Milgrom 2013, respectively). In addition to explaining the internal dynamics of LG satellites, MOND inevitably provides a past major tidal interaction in the LG out of which the MW and M31 satellite planes could have formed (Zhao et al. 2013). Recent MOND-based N-body simulations of this interaction suggest that it could produce structures resembling the LG satellite planes (Bílek et al. 2017). However, it is not clear if this scenario can explain the well-observed orientations of these planes (Pawlowski & Kroupa 2013). Here, we address this question by searching the parameter space much more thoroughly using restricted N-body models of the MW and M31 over a Hubble time. The galaxies are treated as point masses surrounded by test particle discs. We use a grid method to explore a range of MW–M31 relative tangential velocities consistent with the small observed proper motion of M31 (van der Marel et al. 2012) and the unknown external gravitational field (EF) on the LG (Section 2.2). The grid of models we try out thus contains three parameters – the magnitude and direction of the MW–M31 relative tangential velocity and the EF from large-scale structure in which the LG is embedded. In each case, we look at the orbital pole distribution of the tidal debris around each major LG galaxy. Ideally, we would find a clear clustering of orbital poles around the actually observed satellite plane angular momentum direction for that galaxy. We also consider how much the MW disc is dynamically heated by its interaction with M31. This is to check the feasibility of our scenario for producing the MW thick disc (Gilmore & Reid 1983), a structure which seems to have formed fairly rapidly from its thin disc (Hayden et al. 2015). The thick disc formed 9 ± 1 Gyr ago (Quillen & Garnett 2001), a constraint we consider on the time of the MW–M31 flyby. More recent investigations suggest that there was a burst of star formation at that time (Snaith et al. 2014, fig. 2). An enhanced star formation rate could have further thickened the MW disc through the popping star cluster scenario (Kroupa 2002). The star formation rate of M31 also appears to have been much higher for lookback times |${\gtrsim 9}$| Gyr (Williams et al. 2017b, fig. 12). The disc heating which likely formed the Galactic thick disc appears to have been stronger in the outer parts of the MW, characteristic of a tidal effect (Banik 2014). This may explain why its thick disc has a longer scale length than its thin disc (Jurić et al. 2008; Jayaraman et al. 2013). Around each galaxy, there is usually a clear clustering of tidal debris orbital poles, though its direction generally differs for the MW and M31 as their discs are oriented differently. Sometimes, there is also an orbital pole peak in the opposite direction. Such a counter-rotating peak often contains a significant amount of material, especially for the MW. This is good to some extent because Sculptor counter-rotates within the plane preferred by most MW satellites (Piatek et al. 2006; Pawlowski et al. 2011; Sohn et al. 2017). However, some models clearly have too high a counter-rotating fraction. Thus, we try to find models that have a reasonable fraction of counter-rotating material in addition to yielding orbital pole peaks that match the observed spin vectors of both LG satellite planes. We also expect the MW–M31 orbital plane to align with the HVG plane (Banik & Zhao 2018b, co-ordinates in table 4). In this contribution, our objective is to test whether any single MOND model can match all these constraints and thereby provide a reasonable origin scenario for the MW and M31 satellite planes. Having introduced the problem in this section, we explain how we obtain the MW–M31 trajectory (Section 2) and use it to advance test particles (Section 3). We compare each model to observations using the method described in Section 4. Using a χ2 statistic (equation 35), we identify the best-fitting model. The results of this model are described in Section 5, where we argue that it matches observations as well as can be expected given the deficiencies of our model and observational uncertainties. Our conclusions are given in Section 6. 2 MW–M31 TRAJECTORY 2.1 The two-body force in an external field 2.1.1 Deep-MOND limit We begin by determining the past MW–M31 trajectory in MOND. We use the quasi-linear formulation of MOND (QUMOND; Milgrom 2010) because it is easier to handle numerically. In spherical symmetry, this yields identical results to the more traditional aquadratic Lagrangian version of MOND (AQUAL; Bekenstein & Milgrom 1984) as long as the same interpolating function is used to transition from the Newtonian to the low-acceleration (deep-MOND) regimes. However, even in more complicated situations, AQUAL and QUMOND yield forces that only differ by |$\lesssim 5\hbox{ per cent}$| (Banik & Zhao 2015). N-body simulations of both MOND formulations yield quite similar results (Candlish 2016). QUMOND uses the Newtonian gravitational field |$\boldsymbol {g}_{N}$| to determine the true gravitational field |$\boldsymbol {g}$|⁠. Thus, the first step is to find |$\boldsymbol {g}_{N}$| by solving the usual Poisson equation sourced by the baryonic density ρb. \begin{eqnarray*} \nabla \cdot \boldsymbol {g}_{\rm N} \ =\ -4\pi G \rho _{b}. \end{eqnarray*} (3) To determine |$\boldsymbol {g}$|⁠, this Poisson equation must be solved again with a forcing analytically dependent on |$\boldsymbol {g}_{\rm N}$|⁠. Thus, QUMOND does not require a non-linear grid relaxation stage, a major improvement on AQUAL. \begin{eqnarray*} \overbrace{\nabla \cdot {g}}^{\propto \rho _{{\rm {PDM}}} + \rho _{b}} \ &=&\ \nabla \cdot \left[{g}_{\rm N} \nu \left( g_{\rm N} \right) \right], \ \ \text{ where} \end{eqnarray*} (4) \begin{eqnarray*} \nu \left( g_{\rm N} \right) &=& \frac{1}{2} \ +\ \sqrt{\frac{1}{4} + \frac{a_{0}}{g_{\rm N}}}. \end{eqnarray*} (5) Here, ν(gN) is the interpolating function used to transition between the Newtonian and deep-MOND regimes. We use the ‘simple’ form of this function (Famaey & Binney 2005). This works very well at explaining the rotation curve of the MW (Iocco, Pato & Bertone 2015, fig. 3) and the kinematics of ≈6000 elliptical galaxies probing accelerations up to ≈30a0 (Chae, Bernardi & Sheth 2017). The source term for the gravitational field is |$\nabla \cdot \left( \nu {g}_{\rm N} \right)$|⁠, which can be thought of as an ‘effective’ density ρ composed of the actual density ρb and an extra term which we define to be the phantom dark matter density ρPDM. This is the distribution of dark matter that would be necessary in Newtonian gravity to generate the same total gravitational field as QUMOND yields from the baryons alone. To get the gravity between the MW and M31 at some separation dM31, we assume that they are point masses in the deep-MOND regime, though we will later relax this assumption when their relative acceleration |$\ddot{d}_{\rm{M31}}$| is a significant fraction of a0 (Section 2.1.2). If the MW and M31 were the only objects in the Universe, then Zhao, Li & Bienaymé (2010) showed that we could approximate their mutual gravitational acceleration as \begin{eqnarray*} g_{\text {iso}} \ &=&\ \frac{Q\sqrt{GM_{{\text {LG}}}a_0}}{d_{{\rm M31}}}, \end{eqnarray*} (6) \begin{eqnarray*} Q \ &=&\ \frac{2 \left( 1 - {q_{1}}^\frac{3}{2} - {q_{2}}^\frac{3}{2}\right)}{3 q_{1} q_{2}}. \end{eqnarray*} (7) Here, the total MW and M31 mass is MLG, of which a fraction q1 is in the MW and the remaining q2 resides in M31. Equation (6) suggests that |$\ddot{d}_{{\rm M31}} {\approx 0.017\ a_{0}}$|⁠. However, the MOND field equations imply that |$\ddot{d}_{{\text M31}}$| is affected by the EF in which the LG is embedded (Milgrom 1986). This is true even if the EF is due to objects so far outside the LG that they affect the MW and M31 equally (i.e. there are no tides). Including this effect, the deep-MOND limit gravity gM between the MW and M31 can be written as \begin{eqnarray*} g_{M} &=& \frac{\sqrt{GMa_0}}{r} f\left( \tilde{r} \right), \end{eqnarray*} (8) \begin{eqnarray*} \tilde{r} &\equiv & \frac{d_{{\text M31}}}{r_t}, \end{eqnarray*} (9) \begin{eqnarray*} r_t &=& \frac{\sqrt{GMa_0}\left( 3 + \cos ^2 \phi \right)}{4Qg_{\text{ext}}}. \end{eqnarray*} (10) We have assumed that the form of |$f\left( \tilde{r} \right)$| is independent of the angle ϕ between the MW–M31 axis and |$\boldsymbol{\widehat{{g}}}_{{\text {ext}}}$|⁠, with the angular dependence introduced by the EF accounted for by making rt angle dependent. The function |$f\left( \tilde{r} \right)$| should be Q when |$\tilde{r} \ll 1$| or equivalently when gext ≪ giso, representing the isolated regime (Zhao et al. 2010). In the opposite EF-dominated case, |$f \rightarrow \frac{Q}{\tilde{r}}$| but the dependence on Q is cancelled by how this also affects rt (Banik & Zhao 2015, equation 37). It is in this limit that f becomes angle dependent as the MW and M31 start to feel the preferred direction introduced by the EF. The EF on the whole LG is gext ≈ 0.03 a0 (Famaey, Bruneton & Zhao 2007), suggesting that it may have a significant effect. This is especially true as the MW and M31 spend a significant amount of time near apocentre, where the EF is expected to be particularly important. However, the EF is clearly not dominant, i.e. much stronger than giso (equation 6). To include the EF and determine |$f\left( \tilde{r} \right)$|⁠, we thus need a numerical scheme. For this purpose, we use a method similar to that outlined in section 2.2 of Banik & Zhao (2018a). We first obtain |${\boldsymbol{ g}}_{\text N}$| using the superposition principle (i.e. the force from a combination of masses is the sum of the forces due to each one on its own). This is done by adding the Newtonian-equivalent external field |${\boldsymbol{ g}}_{{\text {N, ext}}}$| to the forces from the two point masses. |${\boldsymbol{ g}}_{{\text {N, ext}}}$| is what the EF on the LG would have been if the Universe had the same matter distribution as at present but was governed by Newtonian gravity rather than QUMOND. Assuming the EF is caused by a distant point mass (something we justify in Section 2.3), we can use the spherically symmetric relation between |${\boldsymbol{ g}}_{{\text {N, ext}}}$| and the true EF of |${\boldsymbol{ g}}_{{\text {ext}}}$|⁠. \begin{eqnarray*} \overbrace{\nu \left(\frac{\left| {\boldsymbol{ g}}_{{\text {N, ext}}} \right|}{a_{0}} \right)}^{\nu _{\text {ext}}} {\boldsymbol{ g}}_{\text{N, ext}} &=& {\boldsymbol{ g}}_{\text {ext}}. \end{eqnarray*} (11) In the deep-MOND limit, this reduces to \begin{eqnarray*} g_{\text {N,ext}} \ =\ \frac{{g_{\text {ext}}}^2}{a_{0}}. \end{eqnarray*} (12) Having obtained |${\boldsymbol{ g}}_{\text N}$|⁠, we can use equation (4) to obtain |$\nabla \cdot {\boldsymbol{ g}}$|⁠, the effective matter distribution. We apply a direct summation procedure to this in order to obtain |${\boldsymbol{ g}}$|⁠. \begin{eqnarray*} {\boldsymbol{ g }\left({\boldsymbol{ r}} \right)} \ =\ \int \nabla \cdot {\boldsymbol{ g}} \left( {\boldsymbol{ r}^{\prime }}\right) \frac{\left({\boldsymbol{ r}} - {\boldsymbol{ r}^{\prime }} \right)}{4 \pi |{\boldsymbol{ r}} - {\boldsymbol{ r}^{\prime }}|^3} \ {\rm d}^3 {\boldsymbol{ r}^{\prime }}. \end{eqnarray*} (13) Due to the axisymmetry of the problem, the phantom dark matter can be thought of as a large number of uniform density rings. To find |$\boldsymbol {g}$| on the MW and M31 and thus their relative acceleration, we sum the Newtonian gravity contributed by each of these rings. This is simple to do as the MW and M31 both lie on the symmetry axis. We use equation (50) to estimate that the mass ratio between the MW and M31 is such that q1 = 0.3 based on the flatline level of the MW rotation curve being 180 km s−1 (Kafle et al. 2012) and that of M31 being 225 km s−1 (Carignan et al. 2006). With q1 fixed as we vary gext, our numerical results can be summarized as \begin{eqnarray*} f\left( \tilde{r} \right) = \frac{Q}{\sqrt{1 + \tilde{r}^2}}\left[ 1 + \overbrace{0.032\ {\rm {exp}} \left({-\frac{\left( {\rm ln}\ \tilde{r} - 0.42\right)^2}{2\times 0.65^2}}\right)}^u\right]. \end{eqnarray*} (14) We determine |$f\left( \tilde{r} \right)$| by performing 18 force calculations which vary |$\tilde{r}$| over the range (0−12.5), thus probing a similar range in gext/giso. We fit |$f\left( \tilde{r} \right)$| to these numerical results using three degrees of freedom corresponding to the maximum fractional deviation of f from our initial guess of |$\frac{Q}{\sqrt{1 + \tilde{r}^2}}$|⁠, the value of |$\tilde{r}$| where this occurs and how quickly f approaches our initial guess on either side. It is evident that we accurately recover the expected behaviour in the asymptotic limits (Fig. 1). With 18 force calculations and 3 degrees of freedom in our fitting function, we obtain a total χ2 of 16.1. Therefore, equation (14) provides a good fit to our numerical results. Figure 1. Open in new tabDownload slide Natural logarithm (ln) of |$f\left( \tilde{r} \right)$| (equation 8) against |${\rm {ln}}\ \tilde{r}$|⁠, for two point masses with a 3:7 mass ratio and a total mass M embedded in an external field aligned with the masses. Equivalently, the graph shows the deep-MOND limit force between the masses relative to the naive expectation of equation (1), which is off by a constant factor Q = 0.794 for |${\tilde{r} \ll 1}$| (weak external field) but becomes much less accurate in a strong external field. Our numerical results (red dots with uncertainties) reveal the precise form of the transition between these analytic limits, with the case of no external field |$\left( \tilde{r} = 0 \right)$| arbitrarily plotted at |$\text {ln}\ \tilde{r} = -2.5$| instead of the correct value of −∞. All forces shown here are determined to within 1 per cent. Indistinguishable results are obtained if we set |${g}_{{\rm {ext}}} \rightarrow -{g}_{\rm{ext}}$|⁠. We fit these results using equation (14) (blue line). The total χ2 is 16.1 for 18 points and 3 fit parameters, indicating an acceptable fit. Figure 1. Open in new tabDownload slide Natural logarithm (ln) of |$f\left( \tilde{r} \right)$| (equation 8) against |${\rm {ln}}\ \tilde{r}$|⁠, for two point masses with a 3:7 mass ratio and a total mass M embedded in an external field aligned with the masses. Equivalently, the graph shows the deep-MOND limit force between the masses relative to the naive expectation of equation (1), which is off by a constant factor Q = 0.794 for |${\tilde{r} \ll 1}$| (weak external field) but becomes much less accurate in a strong external field. Our numerical results (red dots with uncertainties) reveal the precise form of the transition between these analytic limits, with the case of no external field |$\left( \tilde{r} = 0 \right)$| arbitrarily plotted at |$\text {ln}\ \tilde{r} = -2.5$| instead of the correct value of −∞. All forces shown here are determined to within 1 per cent. Indistinguishable results are obtained if we set |${g}_{{\rm {ext}}} \rightarrow -{g}_{\rm{ext}}$|⁠. We fit these results using equation (14) (blue line). The total χ2 is 16.1 for 18 points and 3 fit parameters, indicating an acceptable fit. 2.1.2 Going beyond the deep-MOND limit The mutual acceleration between the MW and M31 can be estimated using equation (14), a fit to numerical results in the deep-MOND limit. Although this is valid most of the time, it becomes inaccurate when the MW and M31 are near pericentre. At this time, their larger mutual attraction simplifies the problem in another way – it becomes more accurate to neglect the EF on the LG. Thus, we first estimate the Newtonian attraction between the MW and M31 with total mass MLG. \begin{eqnarray*} g_{\rm N} \ =\ \frac{GM_{\rm {LG}}}{{d_{\rm {M31}}}^2}. \end{eqnarray*} (15) We need to determine whether it is more important to consider the EF or the violation of the deep-MOND limit. We do this by comparing dM31 to the geometric means of the MW–M31 separations at which gN = a0 and where gN = gN, ext, the Newtonian-equivalent EF on the LG (equation 11). For a weak EF, this distinction of cases can be achieved simply by determining whether gN < gext. If it is, then the EF is more important. As gext ≪ a0, it is not as important to consider violation of the deep-MOND limit. In this case, we therefore use our previously derived result for gM including the EF but assuming the deep-MOND limit (equation 8). However, if gN > gext, it is not as important to consider the EF because gext ≫ gN, ext for a weak EF (equation 12). In this case, we neglect the EF but do not assume the deep-MOND limit. This requires a way to transition between that limit and the high-acceleration Newtonian limit. We do this using the ‘simple’ interpolating function (Famaey & Binney 2005) in another sequence of numerical MW–M31 mutual gravity calculations, each of which involves solving equation (4). This extra step is necessary to better determine the MW–M31 trajectory in the period around the MW–M31 encounter as it is not reasonable to assume the deep-MOND limit then. As the EF is time-dependent (Section 2.2), it is difficult to obtain |$\ddot{d}_{\rm {M31}}$| including both the EF and a variable MW–M31 separation. This can only be done by having a 2D grid of numerical simulations, as there are two physical scales in the problem – the MOND radius (equation 1) and the EF radius governing the transition from isolated to EF dominated (where gN = gN, ext). However, we feel that it is unnecessary to consider both scales simultaneously because they are well separated, a consequence of the EF being much weaker than a0. Thus, if the EF is important, then by definition we must have that |$g \lesssim g_{\rm {ext}}$|⁠. This implies that g ≪ a0, making it acceptable to assume the deep-MOND limit. Conversely, a high-acceleration situation will only mildly be affected by a weak EF. 2.2 The external field on the Local Group A viable EF history gext(t) must explain the present peculiar velocity |${\boldsymbol{ v}}_{\rm {pec}}$| of the LG, i.e. its deviation from the Hubble flow. \begin{eqnarray*} {\boldsymbol{ v}}_{\rm {pec}} \ \equiv \ \dot{{\boldsymbol{ r}}}_{\rm {LG}} - H{\boldsymbol{ r}}_{\rm {LG}}. \end{eqnarray*} (16) Here, |${\boldsymbol{ r}}_{\rm {LG}}$| is the position of the LG relative to a comoving observer. We use an overlying dot to indicate a time derivative (e.g. |$\dot{q} \equiv \frac{\mathrm{\partial} q}{\mathrm{\partial} t}$| for any quantity q). The Hubble parameter |$H \equiv \frac{\dot{a}}{a}$|⁠, where a(t) is the cosmic scale factor relative to its present value (i.e. a ≡ 1 at the current time t = tf). For the LG, |${v}_{\rm {pec}}$| has a value close to 630 km s−1 towards Galactic co-ordinates (276°, 30°) when measured with respect to the rest frame of the cosmic microwave background (CMB; Kogut et al. 1993). For reasons discussed in Section 2.3.1, we assume that this has always been the EF direction |${\boldsymbol{\widehat{g}}}_{\rm {ext}} \equiv {\boldsymbol{g}}_{ext} \div \left| {\boldsymbol{g}}_{\rm {ext}} \right|.$| Peculiar velocities must be caused by the time-integrated effect of |${\boldsymbol{ g}}_{\rm {ext}} \left( t \right)$|⁠. Assuming that the EF arises due to very distant objects that have always been in the same direction as viewed from the LG, we get that \begin{eqnarray*} \int _{0}^{t_f} a\ g_{\rm {ext}}\ {\rm{d}}t \ =\ v_{\rm {pec}}. \end{eqnarray*} (17) The integrating factor a(t) accounts for the effect of Hubble drag, which tends to reduce peculiar velocities in the absence of further forcing (as derived in e.g. Banik & Zhao 2016, section 2.1).2 There is no unique way to solve for gext(t), but neither is it wholly unconstrained. To make further progress, we assume that gextis due to a very distant point mass. Thus, \begin{eqnarray*} g_{\rm {ext}} \propto \left\lbrace \begin{array}{ll}\frac{1}{a} &\quad \mbox{if } a \ge a_{\rm {peak}} \\\left( a - a_c \right) &\quad \mbox{if } a_c < a < a_{\rm {peak}}. \\0 &\quad \mbox{if } a \le a_c \end{array} \right. \end{eqnarray*} (18) It is not possible that gext ∝ a−1 indefinitely into the past because the mass must correspond to an inhomogeneity in the Universe which was much smaller in the past. Thus, we make the further assumption that at very early times (a < ac = 0.1), gext = 0. To avoid a discontinuous jump in gext, we assume that gext  ∝  (a − ac) when ac < a < apeak. To transform gext(a) into gext(t) and thus solve equation (17), we need a specific form for the background cosmology a(t). We use a standard flat3 dark energy-dominated cosmology with parameters given in Table 1. \begin{eqnarray*} \frac{\overset{..}{a}}{a} &=& -\frac{4 \rm {\pi } G}{3} \left( \rho _{\rm m} - 2 \rho _\Lambda \right) \end{eqnarray*} (19) \begin{eqnarray*} &=& {H_{0}}^2 \left(- \frac{1}{2} \Omega _{{\rm{m}},0}\ a^{-3} + \Omega _{\Lambda , 0} \right). \end{eqnarray*} (20) Table 1. Cosmological parameters used in our work (Planck Collaboration XIII 2016). The critical density of the Universe at any epoch is |$\frac{3H^2}{8 \pi G}$|⁠. Present values are indicated with 0 subscripts. Our MW–M31 trajectory calculations start when a = 0.05, which occurs 193 Myr after the big bang (equation 21). Name Meaning and units Value Ωm, 0 Fraction of critical density in matter 0.315 ΩΛ, 0 Fraction of critical density in dark energy 0.685 H0 Hubble constant (km s−1 Mpc−1) 67.3 Name Meaning and units Value Ωm, 0 Fraction of critical density in matter 0.315 ΩΛ, 0 Fraction of critical density in dark energy 0.685 H0 Hubble constant (km s−1 Mpc−1) 67.3 Open in new tab Table 1. Cosmological parameters used in our work (Planck Collaboration XIII 2016). The critical density of the Universe at any epoch is |$\frac{3H^2}{8 \pi G}$|⁠. Present values are indicated with 0 subscripts. Our MW–M31 trajectory calculations start when a = 0.05, which occurs 193 Myr after the big bang (equation 21). Name Meaning and units Value Ωm, 0 Fraction of critical density in matter 0.315 ΩΛ, 0 Fraction of critical density in dark energy 0.685 H0 Hubble constant (km s−1 Mpc−1) 67.3 Name Meaning and units Value Ωm, 0 Fraction of critical density in matter 0.315 ΩΛ, 0 Fraction of critical density in dark energy 0.685 H0 Hubble constant (km s−1 Mpc−1) 67.3 Open in new tab Here, ρm and |$\rho _\Lambda$| are the cosmic mean densities of matter and dark energy, respectively. Currently, these components comprise fractions Ωm, 0 and ΩΛ, 0 of the cosmic critical density |$\frac{3{H_{0}}^2}{8 \pi G}$|⁠. Defining time t to start when a = 0 and requiring that |${\overset{.}{a} = H_{0}}$| when a = 1, we obtain the unique solution \begin{eqnarray*} a \left( t \right) &=& {{\left( \frac{{{\Omega }_{{\rm {m}},0}}}{{{\Omega }_{\Lambda ,0}}} \right)}^{\frac{1}{3}}}{{\sinh }^{\frac{2}{3}}}\left( \frac{3}{2}\sqrt{{{\Omega }_{\Lambda ,0}}}{{H}_{0}}t \right). \end{eqnarray*} (21) The present values of H0 and Ωm, 0 uniquely determine the present age of the Universe tf via inversion of this equation to solve for when a = 1. We also use it to determine when a = 0.05, thereby fixing the start time of our simulations. Equation (21) is valid for ΛCDM but may be inaccurate for a MOND model. However, we expect only small deviations from our adopted cosmology because this works fairly well observationally (e.g. Joudaki et al. 2017). Moreover, even a purely baryonic universe could have an expansion history very similar to ΛCDM (Zhao 2007). In our model, gextpeaks when a = apeak at a value higher than the present value (gext, 0) by the factor |$\frac{1}{a_{\rm {peak}}}$|⁠. We consider models to be realistic if apeak ≈ 0.77, when the Universe stopped being matter dominated. This requires gext,0 ≈ 0.023 a0, though we try a range of values to see which works best (Section 5). 2.3 Tides on the Local Group 2.3.1 Perturber trajectories Any method for determining tides exerted by massive perturbers on the LG requires us to estimate where these perturbers were at earlier times. To do this, we treat the LG as a single point mass at the MW–M31 barycentre. We then use equation (8) to get the relative gravitational acceleration between the LG particle at position |${\boldsymbol{ r}}_{\rm{LG}}$| and a perturber at position |${\boldsymbol{ r}}_{\rm {p}}$|⁠, neglecting the other perturbers (i.e. applying the superposition principle to the effects of the perturbers, something we will justify shortly). As the perturbers we consider are all very distant (Table 2), we assume that each perturber has always been in the same direction as perceived from the LG barycentre |${\boldsymbol{ r}}_{\rm {LG}}$|⁠. The distance dp of the perturber from |${\boldsymbol{ r}}_{\rm{LG}}$| is thus the only dynamical quantity. We vary its initial value at some early time ti in order to match its presently observed value. This is a basic form of timing argument analysis (Kahn & Woltjer 1959), a technique we will use more extensively in Section 2.5. \begin{eqnarray*} \ddot{d}_{\rm p} \ &=&\ \frac{\ddot{a}}{a} d_{\rm p} \ +\ g_{{M}}, \end{eqnarray*} (22) \begin{eqnarray*} \dot{d}_{\rm p} \ &=&\ H d_{\rm p} \ \ \text{initially}. \end{eqnarray*} (23) Table 2. Properties of mass concentrations outside the Local Group which we consider. Sky positions are shown using Galactic co-ordinates. The Galactocentric distances dMW are given in Mpc for Cen A (Harris, Rejkuba & Harris 2010), M81 (Gerke et al. 2011), and IC 342 (Wu et al. 2014). Baryonic masses M in units of 1011 M⊙ are from Karachentsev (2005) for Cen A and IC 342, while Karachentsev & Kashibadze (2006) is used for M81. Name b l dMW M Centaurus A 19.4173° 309.5159° 3.8 6 M81 40.9001° 142.0918° 3.6 2.2 IC 342 10.5799° 138.1726° 3.45 2.2 Name b l dMW M Centaurus A 19.4173° 309.5159° 3.8 6 M81 40.9001° 142.0918° 3.6 2.2 IC 342 10.5799° 138.1726° 3.45 2.2 Open in new tab Table 2. Properties of mass concentrations outside the Local Group which we consider. Sky positions are shown using Galactic co-ordinates. The Galactocentric distances dMW are given in Mpc for Cen A (Harris, Rejkuba & Harris 2010), M81 (Gerke et al. 2011), and IC 342 (Wu et al. 2014). Baryonic masses M in units of 1011 M⊙ are from Karachentsev (2005) for Cen A and IC 342, while Karachentsev & Kashibadze (2006) is used for M81. Name b l dMW M Centaurus A 19.4173° 309.5159° 3.8 6 M81 40.9001° 142.0918° 3.6 2.2 IC 342 10.5799° 138.1726° 3.45 2.2 Name b l dMW M Centaurus A 19.4173° 309.5159° 3.8 6 M81 40.9001° 142.0918° 3.6 2.2 IC 342 10.5799° 138.1726° 3.45 2.2 Open in new tab Here, gM is the mutual gravitational acceleration between the perturber and the LG according to equation (8), which assumes the deep-MOND limit. The relevant mass M is now the total of the LG and the perturber. Their vector separation and mass ratio enter the problem in the same way as we used these quantities for the MW and M31 in Section 2.1.1. Because the Universe was expanding nearly homogeneously at very early times (Planck Collaboration XIII 2016), we assume that the LG and the perturber both started exactly on the Hubble flow when our simulations start (redshift 19 or a = 0.05). In this problem, the deep-MOND limit is a very good approximation due to the large distance between the LG and perturbers like Cen A. It should also be reasonable to neglect the angular momentum of the LG as a whole with respect to these rather distant (⁠|${\gtrsim 3}$| Mpc) perturbers. Although we numerically calibrate the form of gM at a particular mass ratio between the objects, we assume it is valid for the mass ratios between the LG and the perturbers listed in Table 2. This is partly because the mass ratios are similar to the 3:7 case we investigated in detail. Moreover, the greater distances imply that the EF is more significant. In the limit where the EF dominates, we do not need a numerical scheme to determine the force between the LG and a perturber. In this limit, the superposition principle also becomes valid as the phantom dark matter density is linear in the matter distribution if it is embedded in a dominant EF (Banik & Zhao 2015, equation 25). This is because a perturbation theory approach becomes valid, suppressing the non-linearity usually inherent to MOND. When determining gM, we need to know the EF acting on the LG–perturber system. In principle, some of the EF on the LG could arise due to the perturber itself such that it affects the MW–M31 mutual attraction but should not be considered when determining the LG–perturber attraction. For simplicity, we assume that the EF on the LG mostly arises due to very distant objects rather than the relatively nearby (within 5 Mpc) perturbers we consider (Table 2). This can be justified by considering the effects of Cen A, the perturber likely to impose the largest EF on the LG due to its high mass and relative proximity. Supposing that most of the EF on the LG arose due to Cen A, we can consider the LG–Cen A problem in isolation to put an upper limit on their mutual acceleration. In this case, the gravity on the LG is |$g_{{\rm {iso}}} \approx \frac{\sqrt{GMa_{0}}}{d_{\rm p}}$|⁠, where M is the total (LG + Cen A) mass. For a distance to Cen A of dp = 3.8 Mpc, this yields 0.0070a0. As the actual EF on the LG is very likely larger (Section 2.2), the LG–Cen A problem cannot be treated in isolation. If instead we assume that the EF from more distant objects dominates and set gext = 0.02a0, then we get that Cen A pulls on the LG by |$g_{\rm {EFE}} \approx \frac{GMa_{0}}{{d_{\rm p}}^2 g_{\rm {ext}}}$|⁠, i.e. by only 0.0024a0. As this is much less than gext, it is self-consistent to assume that the EF on the LG mostly arises from objects much further away than the perturbers we consider. However, it is not self-consistent to suppose that most of this EF is provided by nearby objects. As Cen A is the most massive perturber we consider, this justifies our approximation that each LG–perturber problem can be considered embedded in the dominant EF on the LG found in Section 2.2. As a result, it becomes valid to superpose the forces due to each perturber. For simplicity, we neglect how perturbers affect each other as they are all rather distant and in different sky directions (Table 2). When running our algorithm for the first time, we treat the LG as a single point mass because we do not yet know dM31(t). The perturber trajectories thus found allow us to obtain the MW–M31 trajectory. We take advantage of this by running the algorithm on a second pass, this time determining what gM would be if the LG was at the position of the MW. After repeating this calculation for the M31 position, we take a weighted average of these determinations according to the MW:M31 mass ratio. \begin{eqnarray*} &&g_{{M}} =\nonumber \\&& \frac{\sqrt{GMa_0}}{d_{\rm p}} f\left( \tilde{r} \right) \left( 1 + \frac{q_{\rm {MW}}q_{\rm {M31}}n\left(n\cos ^2 \theta + 1\right){d_{\rm {M31}}}^2}{2{d_{\rm p}}^2}\right), \end{eqnarray*} (24) \begin{eqnarray*} n \equiv 1 \ -\ \frac{\mathrm{\partial} \,{\rm {ln}}\ f \left( \tilde{r} \right)}{\mathrm{\partial} \tilde{r}} \ \ \left( g_{{M}} \propto_{_{_{\!\!\!\!\sim}}} {\tilde{r}}^{-n} \right) \end{eqnarray*} (25) \begin{eqnarray*} = 1 \ +\ \frac{\tilde{r}^2}{\tilde{r}^2 + 1} \ +\ \frac{\left( {\rm ln}\ \tilde{r} - 0.42 \right)u}{0.65^2\left(1+u\right)}, \end{eqnarray*} (26) \begin{eqnarray*} \tilde{r} &=& \frac{4d_{\rm p}Qg_{\rm {ext}}}{\sqrt{GMa_0}\left( 3 + \left( \widehat{\boldsymbol{r}} \cdot {\widehat{\boldsymbol{ g}}}_{\rm {ext}}\right)^2 \right)}. \end{eqnarray*} (27) Here, M is the total mass of the LG and the perturber, while Q is calculated for the ratio between their masses. qMW and qM31 retain their previous meanings and values. u is found according to equation (14) based on the scaled distance |$\tilde{r}$|⁠. The second-order correction to gM found here depends on the angle θ between the MW−M31 line and the direction from the LG barycentre to the perturber. Using the long bracketed term in equation (24) improves the accuracy of |$\ddot{d}_{\rm p}$| and allows us to refine their trajectory. We only do this once as it leads to fairly small corrections. The final perturber trajectories dp(t) are shown in Fig. 2 for our best-fitting model (Section 5). Figure 2. Open in new tabDownload slide Trajectories of the perturbers in Table 2, i.e. their distance from the LG barycentre |${r}_{\rm {LG}}$| over time. The MW is shown here for comparison (M31 is much closer to |${r}_{\rm {LG}}$| due to its higher mass). Figure 2. Open in new tabDownload slide Trajectories of the perturbers in Table 2, i.e. their distance from the LG barycentre |${r}_{\rm {LG}}$| over time. The MW is shown here for comparison (M31 is much closer to |${r}_{\rm {LG}}$| due to its higher mass). 2.3.2 Effect of perturbers on the MW–M31 trajectory We estimate the effect of tides raised by each perturber on the LG based on how much their mutual attraction alters if the LG particle is placed at the location of the MW instead of M31. This difference in forces is projected on to the MW–M31 line to estimate how each perturber affects the MW–M31 mutual acceleration. For simplicity, we neglect the fact that there would also be a tidal torque on the MW–M31 orbit which would change its angular momentum. In general, the gravitational force between two objects is not along the line connecting them (e.g. Banik & Zhao 2015, fig. 1). When handling the trajectory of the LG with respect to a perturber, we neglect this because we expect distant perturbers to be on nearly radial Hubble flow orbits. However, when determining the effect of each perturber on the LG, the important consideration is how this differs depending on whether |${r}_{\rm {LG}} = {r}_{\rm {MW}}$| or |${r}_{\rm {LG}} = {r}_{\rm {M31}}$|⁠. Thus, it may be important to consider the tangential force, i.e. the force on the LG orthogonal to the LG−perturber line. We find this by determining the gravity due to a point mass embedded in a constant EF. At long range, the EF dominates and so the gravitational field becomes analytic (Banik & Zhao 2015, equation 37). To get the tangential force at smaller distances, we construct a force library for a point mass embedded in an EF in the deep-MOND limit. In this way, we can find the tangential force due to the mass, something without a Newtonian analogue. Our results are only valid if one of the objects is much less massive than the other. At present, we have no way of handling intermediate mass ratios as we can only rigorously solve problems that are at least axisymmetric. Thus, we guess how the relative tangential gravity |${g}_{\rm {tan}}$| between two objects varies with their mass ratio (Appendix A). In this way, we estimate the gravity |$\boldsymbol {g}_{\rm{LG}}$| exerted by each perturber on the MW and M31, both directly towards the perturber and in the orthogonal direction (within the plane defined by the EF and the LG–perturber line). The difference between the gravity on the MW and M31 is then projected on to the line connecting them so as to estimate gtide, our symbol for how much each perturber tidally alters the MW–M31 mutual attraction |$\ddot{d}_{\rm {M31}}$|⁠. \begin{eqnarray*} g_{\rm {tide}} &=& \left( \left. \boldsymbol {g}_{\rm {LG}} \right|_{\boldsymbol {r}_{\rm {LG}} = \boldsymbol {r}_{\rm {M31}}} - \left. \boldsymbol {g}_{\rm {LG}} \right|_{\boldsymbol {r}_{\rm {LG}} = \boldsymbol {r}_{\rm {MW}}} \right) \cdot \widehat{\boldsymbol {d}}_{\rm {M31}} \nonumber, \\\widehat{\boldsymbol {d}}_{\rm {M31}} &=& \frac{\boldsymbol {r}_{\rm {M31}} - \boldsymbol {r}_{\rm {MW}}}{\left| \boldsymbol {r}_{\rm {M31}} - \boldsymbol {r}_{\rm {MW}} \right|}. \end{eqnarray*} (28) We assume that the EF on the LG arises due to a mass concentration similar to the GA (Mieske, Hilker & Infante 2005). Thus, we take the EF to be entirely caused by a point mass currently 84 Mpc away and always in the same direction as viewed from the LG barycentre (Section 2.2). Because of the large distance to the GA, we assume that its distance from the LG scaled directly with a(t), i.e. it was always on the Hubble flow. Although this may not have been true at very early times, forces then have only a very small effect on present peculiar velocities due to Hubble drag (equation 17). The actual effect is even more pronounced (Banik & Zhao 2016, fig. 4). At 84 Mpc, the GA is expected to raise only a small tide on the LG. However, for completeness, we include this using the distant tide approximation. In this scheme, we approximate the tidal effect of the EF as an extra contribution of \begin{eqnarray*} g_{\rm{tide,GA}} \ =\ \frac{\left( 2\cos ^2 \theta _{\rm {ext}} - 1 \right)d_{\rm{M31}}g_{\rm{ext}}}{d_{\rm{GA}}}. \end{eqnarray*} (29) Here, θext is the angle between the MW–M31 line and |$\boldsymbol {\widehat{g}}_{\rm {ext}}$|⁠. The factor of (2cos 2θext − 1) is appropriate for an inverse distance gravity law (equation 1). This applies to the LG–GA problem due to the very low accelerations involved and our assumption that it can be considered in isolation. 2.4 The angular momentum barrier The MW–M31 trajectory cannot be approximated as radial if they undergo a past close encounter. For simplicity, we assume that their mutual orbital angular momentum |${\boldsymbol {h}}_{\rm {LG}}$| has remained constant since their first turnaround when the Universe was  ∼ 3 Gyr old. This is because |${\boldsymbol {h}}_{\rm{LG}}$| must have been gained by tidal torques, which are most significant when the MW and M31 are furthest apart. Although we expect the MW and M31 to have gone through two apocentres (Banik & Zhao 2018b, fig. 4), any large-scale structures are much further away in the latter epoch, suggesting that tidal torques then are less significant than at the first MW–M31 turnaround. Moreover, we will see later that the main constraint on our simulations arises from |${\boldsymbol {h}}_{\rm{LG}}$| around the time of the MW–M31 flyby. Their present proper motion is not constrained all that well and is in any case dependent on other uncertain quantities like the masses of the LMC and M33 relative to the MW and M31, respectively (Section 5.1). Thus, our analysis would not be much affected if a dwarf galaxy recently merged with M31, totally changing |$\widehat{\boldsymbol {h}}_{\rm{LG}}$| as |${\boldsymbol {h}}_{\rm{LG}}$| is so small.4 In this case, our best-fitting model would have the same |${\boldsymbol {h}}_{\rm{LG}}$|⁠, the only difference being that this would be slightly less consistent with the observed proper motion of M31. Thus, it is beyond the scope of our analysis to consider whether |${\boldsymbol {h}}_{\rm{LG}}$| has changed since the MW–M31 flyby. Instead, we consider all possible orientations for |$\widehat{\boldsymbol {h}}_{\rm{LG}}$| and a wide range of values for |$\left| {\boldsymbol {h}}_{\rm{LG}} \right|$|⁠. This ensures that the plausible ranges are covered, thus giving MOND a fair chance of explaining the LG satellite planes. 2.5 Satisfying the timing argument Combining the contributions to the MW–M31 mutual attraction |$\ddot{d}_{\rm{M31}}$| from cosmological acceleration (⁠|$\frac{\ddot{a}}{a}d_{\rm{M31}}$|⁠), angular momentum (hLG2/dM313), self-gravity g (Section 2.1), and the tidal contribution gtide due to each perturber (equations 28 and 29), we get that \begin{eqnarray*} \ddot{d}_{\rm{M31}} \ =\ \frac{\ddot{a}}{a}d_{\rm{M31}} \ +\ \frac{{h_{\rm{LG}}}^2}{{d_{\rm{M31}}}^3} \ +\ g \ +\ \sum _{\rm{Perturbers}} g_{\rm{tide}}. \end{eqnarray*} (30) We use this to integrate the MW–M31 trajectory back from present conditions, assuming a particular MW–M31 specific angular momentum |$\boldsymbol {h}_{\rm{LG}}$| in the period after their first turnaround but a purely radial orbit before then. As the MW–M31 line |$\widehat{\boldsymbol {d}}_{\rm{M31}}$| rotates, it changes orientation relative to the EF and the directions towards nearby perturbers. Thus, the results of our simulations depend on |$\widehat{{\boldsymbol {h}}}_{\rm{LG}}$|⁠. \begin{eqnarray*} \dot{\widehat{\boldsymbol {d}}}_{\rm{M31}} \ =\ \frac{\boldsymbol {h}_{\rm{LG}} \times \widehat{\boldsymbol {d}}_{\rm{M31}}}{{{d}_{\rm{M31}}}^2}. \end{eqnarray*} (31) |$\widehat{\boldsymbol {d}}_{\rm{M31}}$| rotates substantially during our simulation, but almost all of this rotation occurs close to the time of the MW–M31 flyby as dM31 is rather small then. Thus, it is not too important to accurately know |$\boldsymbol {h}_{\rm{LG}}$| at other times. A critical constraint on our problem is provided by the timing argument (Kahn & Woltjer 1959). The near-uniformity of the CMB indicates that the Universe was expanding almost homogeneously at early times (Planck Collaboration XIII 2016). Thus, our backward integration should leave the MW–M31 separation satisfying \begin{eqnarray*} \dot{d}_{\rm{M31}} \ =\ H_{i}d_{\rm{M31}} \ \ \text{when } t = t_i. \end{eqnarray*} (32) Our simulations cover the period from t = tf (the present epoch) back to t = ti, the limit of our backward integration. We refer to this as the start time of our simulations, even though in a numerical sense they start at tf. We use equation (21) to determine ti by solving it for when a = 0.05 (redshift 19). This allows us to determine the Hubble parameter Hi at that epoch. In general, a particular set of model parameters will not satisfy equation (32). We vary the total LG mass using a Newton–Raphson procedure to ensure this boundary condition is satisfied. This occurs at a particular mass called the timing argument mass. In Fig. 3, we show how dM31(t) is affected by using a lower or higher LG mass. Visually, it is possible to see why these other trajectories are incorrect – loosely speaking, they do not ‘naturally’ (i.e. without the cosmic acceleration term) lead to an MW–M31 pericentre at the time of the big bang. Figure 3. Open in new tabDownload slide Three solutions for the MW–M31 separation dM31(t), all with the same present separation and relative velocity but using different total masses M. Only one of the trajectories (thick black) satisfies the timing argument (equation 32). If M is too low (red curve), |$v_{\rm {pec}} \propto \frac{1}{a} \rightarrow +\infty$|⁠; while if M is too high (blue curve), vpec → −∞ as a → 0. In neither case does dM31 smoothly approach 0 at early times. Figure 3. Open in new tabDownload slide Three solutions for the MW–M31 separation dM31(t), all with the same present separation and relative velocity but using different total masses M. Only one of the trajectories (thick black) satisfies the timing argument (equation 32). If M is too low (red curve), |$v_{\rm {pec}} \propto \frac{1}{a} \rightarrow +\infty$|⁠; while if M is too high (blue curve), vpec → −∞ as a → 0. In neither case does dM31 smoothly approach 0 at early times. 3 THE MW AND M31 DISCS To enable investigation of a wide range of model parameters governing a past flyby encounter between the MW and M31, we model them as point masses surrounded by test particle discs. For simplicity, we neglect the EF but do not assume the deep-MOND limit as this is not accurate close to the MW and M31, regions critical for our analysis. At any instant in time, the gravitational field |$\boldsymbol {g}$| governing our test particles is axisymmetric, though the symmetry axis |$\widehat{\boldsymbol {d}}_{\rm{M31}}$| changes with time according to equation (31). We can thus obtain |$\boldsymbol {g}$| by numerically solving equation (4) using the direct summation method (equation 13). This is straightforward for points on the symmetry axis of the problem. At other points, we reduce the computational burden by utilizing a ‘ring library’ which stores |$\boldsymbol {g}_{N}$| due to a unit radius ring. Thus, instead of having to further split each ring into a finite number of elements, we can simply interpolate within our densely allocated ring library to find the gravity exerted by the ring at the point where we wish to know its contribution to |$\boldsymbol {g}$|⁠. As our numerical scheme is inaccurate very close to the MW and M31, we terminate the trajectories of test particles that get within a disc scale length of either galaxy (Table 3; stellar disc used for MW). This should not affect our results too much because we are mostly interested in their satellites. These lie several disc scale lengths from their hosts, making a point mass approximation reasonably accurate for a theory where galaxies lack extended dark matter haloes. Table 3. Parameters of the MW and M31 discs. We use kpc for the exponential scale lengths of the MW stellar disc (Bovy & Rix 2013), its gas disc (McMillan 2017), and the M31 disc (Courteau et al. 2011). We model the MW disc surface density as a double exponential with 0.1764 of its mass in its gas disc and the rest in its stellar disc (Banik & Zhao 2018a, table 1). M31 is also treated as an exponential disc, though its mass is always |$\frac{7}{3} \times$| the total MW mass. Disc orientations are discussed in Table 5. Disc scale length of Value (kpc) MW stars 2.15 MW gas 7 M31 5.3 Disc scale length of Value (kpc) MW stars 2.15 MW gas 7 M31 5.3 Open in new tab Table 3. Parameters of the MW and M31 discs. We use kpc for the exponential scale lengths of the MW stellar disc (Bovy & Rix 2013), its gas disc (McMillan 2017), and the M31 disc (Courteau et al. 2011). We model the MW disc surface density as a double exponential with 0.1764 of its mass in its gas disc and the rest in its stellar disc (Banik & Zhao 2018a, table 1). M31 is also treated as an exponential disc, though its mass is always |$\frac{7}{3} \times$| the total MW mass. Disc orientations are discussed in Table 5. Disc scale length of Value (kpc) MW stars 2.15 MW gas 7 M31 5.3 Disc scale length of Value (kpc) MW stars 2.15 MW gas 7 M31 5.3 Open in new tab Each test particle is initially set up at position |$\boldsymbol {r}$| on a circular orbit around galaxy i of mass Mi and position |$\boldsymbol {r}_{i}$|⁠. The circular speed is based on a gravitational field directly towards galaxy i of strength \begin{eqnarray*} g &=& \frac{g_{\rm N}}{2} \ +\ \sqrt{\left( \frac{g_{\rm N}}{2} \right)^2 + g_{\rm N}a_{0}}, \end{eqnarray*} (33) \begin{eqnarray*} g_{\rm N} &=& \frac{GM_i}{\left| \boldsymbol {r} - \boldsymbol {r}_{i} \right|^2}. \end{eqnarray*} (34) We add the systemic velocity of galaxy i to the circular velocity |$v_{c} = \sqrt{\left| \boldsymbol {r} - \boldsymbol {r}_{i} \right|g}$|⁠. Our estimates of vc are valid if we neglect the EF and the other galaxy. For a particle very close to one galaxy, this should be a reasonable assumption. However, starting when a = 0.05 leads to a rather small MW–M31 separation, causing some MW particles to actually start closer to M31. Despite starting quite close, the MW and M31 subsequently spend several Gyr at rather large separations. During this time, we expect their discs to dissipatively settle down as they should have been gas rich at early times. However, our simulations lack dissipative processes. To help address this issue, we start with completely dynamically cold MW and M31 discs at a fairly large separation and a relative velocity consistent with the MW–M31 trajectory from Section 2. Thus, we start the test particle integrations 2 Gyr after the big bang. This increases the MW–M31 separation from 200  to 700 kpc, thereby preventing any issues due to M31 significantly perturbing MW particles before the flyby. Our results are not much different if we start integrating the test particle trajectories when a = 0.05, though we consider such an early start less realistic. Our grids of initial test particle positions are uniform in distance from the nearest galaxy and azimuthal angle with respect to it (Table 4). Although our particles do not have any mass, they none the less represent a certain amount of mass in the real world. To account for this, we integrate the assumed surface density law (Table 3) to obtain the mass represented by each particle, later using these masses as statistical weights when analysing our simulations (Section 4). We treat M31 as an exponential disc but use a double exponential model for the MW, whose adopted mass distribution is the same as that used in Banik & Zhao (2018a) to calculate its rotation and escape velocity curves in MOND. Table 4. Resolutions and extents of the test particle grids used to simulate the MW and M31. We terminate trajectories of particles that get within a disc scale length of each galaxy (using the stellar disc for the MW). Parameter Value for MW Value for M31 Radial resolution (kpc) 0.215 0.53 Azimuthal resolution 2.8125° 2.8125° Inner radial limit (kpc) 2.15 5.3 Outer radial limit (kpc) 107.5 106 Parameter Value for MW Value for M31 Radial resolution (kpc) 0.215 0.53 Azimuthal resolution 2.8125° 2.8125° Inner radial limit (kpc) 2.15 5.3 Outer radial limit (kpc) 107.5 106 Open in new tab Table 4. Resolutions and extents of the test particle grids used to simulate the MW and M31. We terminate trajectories of particles that get within a disc scale length of each galaxy (using the stellar disc for the MW). Parameter Value for MW Value for M31 Radial resolution (kpc) 0.215 0.53 Azimuthal resolution 2.8125° 2.8125° Inner radial limit (kpc) 2.15 5.3 Outer radial limit (kpc) 107.5 106 Parameter Value for MW Value for M31 Radial resolution (kpc) 0.215 0.53 Azimuthal resolution 2.8125° 2.8125° Inner radial limit (kpc) 2.15 5.3 Outer radial limit (kpc) 107.5 106 Open in new tab Using our numerically determined gravitational field for the MW–M31 trajectory found in Section 2, we advance test particles using a fourth-order Runge–Kutta method. The time-step is quantized and varied by integer powers of 2 so that there are at least 70 steps per dynamical time, which we define as the ratio of distance to velocity relative to each major LG galaxy. Naturally, we use whichever galaxy returns a smaller dynamical time, regardless of whether this is the galaxy in which the particle started. At early times, cosmic expansion becomes important. Thus, we also require our time-step to be below |$\frac{1}{200}$| of the Universe’s age. If necessary, our algorithm can reduce the time-step to 213 kyr. This should be sufficient to handle very early times and/or particles quite close to the centre of the MW or M31. For more distant particles, a longer time-step is generally used, though this never exceeds 6.8 Myr. 3.1 Disc orientations Due to the MW–M31 flyby, their discs end up in a slightly different orientation to what they start with. To account for this, we adjust the initial disc spin vectors of the MW and M31 so as to match their presently observed orientations. Because our simulations use test particle discs, it is possible to adjust the discs independently. We use a 2D Newton–Raphson procedure to match the presently observed orientation of each disc.5 We determine the final MW disc spin vector in our simulations by fitting a plane to the particles within a cylindrical radius of 10.75 kpc and within 21.5 kpc of the disc plane (using the present MW disc orientation – we are only interested in broadly selecting the relevant region). For M31, we use 21.2 kpc radially and 15.9 kpc vertically. Our plane-fitting procedure is not only similar to that used in section 5.1 of Banik & Zhao (2018b) but also considers the unequal statistical weights of the particles and applies outlier rejection with a threshold of 2.58σ, corresponding to the 99 per cent confidence interval of a Gaussian distribution. As rejecting outliers reduces the root mean square (rms) plane thickness, we repeat the procedure iteratively until the algorithm no longer finds any outliers. Our results do not depend much on the details of the outlier rejection system. To match the presently observed orientations of the MW and M31 discs, we need their initial orientations to be as shown in Table 5. Fortunately, both discs precess by only a few degrees. Thus, we are always sure which hemisphere contains the total angular momentum of each disc. This allows us to get by with a plane-fitting algorithm, though we do also check the angular momenta of the particles to which we are fitting a plane. Table 5 . Initial spin vectors of the MW and M31 discs in Galactic co-ordinates. These differ slightly from their presently observed orientations (last column) because the discs precess during their close encounter. The present M31 disc orientation is obtained from Banik & Zhao (2018b, table 1) and changes by 5.37° while that of the MW changes by only 2.52°. Open in new tab Table 5 . Initial spin vectors of the MW and M31 discs in Galactic co-ordinates. These differ slightly from their presently observed orientations (last column) because the discs precess during their close encounter. The present M31 disc orientation is obtained from Banik & Zhao (2018b, table 1) and changes by 5.37° while that of the MW changes by only 2.52°. Open in new tab We assume that the MW and M31 have a constant mass throughout the entire simulation. We will see later that only a small fraction of their mass ends up in their satellite regions, so the galaxies lose rather little mass in our best-fitting model (Section 5.2.3). In general, including a small amount of mass-loss at the time of the MW–M31 flyby further reduces their total timing argument mass while pushing their flyby further into the past. This makes the results slightly more consistent with observations. However, we neglect this effect to avoid having too many model parameters and because there would also be some accretion on to the MW and M31. 4 COMPARISON WITH OBSERVATIONS We compare each model with observations using a χ2 statistic that includes contributions for the satellite plane orientations |$\chi ^2_{{{\rm{SP}},i}}$| (Section 4.1) and counter-rotating fractions |$\chi ^2_{{{\rm{CR}},i}}$| (Section 4.2) as well as |$\chi ^2_{\rm{TD}}$| for how the MW disc is heated by the flyby (Section 4.3) and |$\chi ^2_{\rm{HVG}}$| (Section 4.4) for the alignment of the MW–M31 orbital plane with the plane of HVGs identified by Banik & Zhao (2018b). \begin{eqnarray*} \chi ^2 = \sum _{i = \text{MW, M31}} \left( \chi ^2_{\rm{SP,i}} + \chi ^2_{\rm{CR,i}} \right) \ +\ \chi ^2_{\rm{TD}} \ +\ \chi ^2_{\rm{HVG}}. \end{eqnarray*} (35) Although the flyby must have affected both the MW and M31 discs, this is much more difficult to test for M31 than for the MW because M31 is an external galaxy viewed nearly face-on (Chemin, Carignan & Foster 2009). The accretion history of M31 is also more complicated and includes recent interactions with other galaxies, significantly complicating any analysis (Sadoun, Mohayaee & Colin 2014). For example, the M31 disc is likely being perturbed by the close proximity to M32 (Weisz et al. 2014, table 4). In future, accurate ages and kinematics of M31 stars could help distinguish the effects of more recent interactions. The work of Williams et al. (2017b) may be an important step in this regard. 4.1 Satellite plane angular momenta The main objective of this investigation is to compare our simulations with observations of the MW and M31 satellite planes. To facilitate this comparison, we first need to select which simulated particles correspond to satellites of e.g. the MW. We define a ‘satellite region’ around it by considering MW particles within an outer cut-off radius Rmax of the MW, as this encompasses the region usually considered in the satellite plane problem. In addition, we discard particles too close to its disc by excluding the region with rplane < rmin and |z| < zmin in cylindrical polar co-ordinates (rplane, ϕ, z) centred on the MW with axis along the observed MW disc spin vector. Our adopted values for these parameters are given in Table 6. We exclude everything within 15 disc scale lengths of its disc plane out to a large distance, completely removing the MW disc remnant. However, most of its simulated satellite plane should still remain outside our excluded region if it is similar to the corresponding observed structure. This is because the MW satellite plane is much larger than zmin and almost polar relative to the Galactic disc (Pawlowski 2018, fig. 1). For M31, we use a similar procedure to the MW, though with a slightly higher zmin to compensate for its larger disc. Table 6. Our adopted filters for selecting simulation particles corresponding to the satellite galaxy regions around the MW and M31. Variables are defined in Section 4.1. Although we could have used rmin = Rmax for both galaxies, in practice the shorter scale length of the MW meant a lower rmin was sufficient. Quantity Value for MW (kpc) Value for M31 (kpc) Rmax 250 250 rmin 172 250 zmin 32.25 53 Quantity Value for MW (kpc) Value for M31 (kpc) Rmax 250 250 rmin 172 250 zmin 32.25 53 Open in new tab Table 6. Our adopted filters for selecting simulation particles corresponding to the satellite galaxy regions around the MW and M31. Variables are defined in Section 4.1. Although we could have used rmin = Rmax for both galaxies, in practice the shorter scale length of the MW meant a lower rmin was sufficient. Quantity Value for MW (kpc) Value for M31 (kpc) Rmax 250 250 rmin 172 250 zmin 32.25 53 Quantity Value for MW (kpc) Value for M31 (kpc) Rmax 250 250 rmin 172 250 zmin 32.25 53 Open in new tab For simplicity, we neglect particles starting in one galaxy and ending up within the satellite region of the other galaxy. There is little tendency for particles to end up very far from their original galaxy (Figs 6 and 7). After selecting ‘satellite’ particles in this way, we create a 2D histogram of the orbital poles |$\widehat{\boldsymbol {h}} \equiv \frac{\boldsymbol {h}}{\left| \boldsymbol {h} \right|}$| of these particles relative to their nearest galaxy, where \begin{eqnarray*} \boldsymbol {h} \ =\ \left(\boldsymbol {r} - \boldsymbol {r}_{i} \right) \times \left(\boldsymbol {v} - \boldsymbol {v}_{i} \right). \end{eqnarray*} (36) Here, the nearest galaxy to the particle is galaxy i with position |$\boldsymbol {r}_{i}$| and velocity |$\boldsymbol {v}_{i}$|⁠. We assign the statistical weight of the particle to the bin in Galactic latitude and longitude defined by its orbital pole |$\widehat{\boldsymbol {h}}$|⁠. This results in a 2D histogram for the MW and another one for M31. Each of these histograms almost always has a clear peak indicating a preferred orbital pole, though its direction varies with model parameters and between the galaxies. Based on the bin that is assigned the most mass, we identify the peak |$\widehat{\boldsymbol {h}}_{\rm{peak}}$| in each orbital pole histogram. We then quantify the angular dispersion of this peak, considering only particles with |$\widehat{\boldsymbol {h}}$| in the same hemisphere as the peak, i.e. with |$\widehat{\boldsymbol {h}} \cdot \widehat{\boldsymbol {h}}_{\rm{peak}} >0$|⁠. We also quantify the angular dispersion of the remaining particles’ orbital poles about the direction |$-\widehat{\boldsymbol {h}}_{\rm{peak}}$|⁠. This allows for the possibility of a counter-rotating peak, i.e. a secondary peak in the orbital pole distribution at |$-\widehat{\boldsymbol {h}}_{\rm{peak}}$|⁠. To make our results more robust, we apply a σ-clipping procedure to reject outliers from the primary (≡ corotating) and counter-rotating peaks. The threshold we use is 2.58σ, corresponding to the 99 per cent confidence interval of a Gaussian distribution.6 The angular dispersion of each peak is defined as the weighted rms angle to |$\pm \widehat{\boldsymbol {h}}_{\rm{peak}}$| of particles whose orbital pole lies in the appropriate hemisphere. Considering the statistical weight Wj of each particle j, the angular dispersion of the counter-rotating peak is thus \begin{eqnarray*} \sigma _{-} \ =\ \sqrt{\frac{\sum _j W_j \left[ \overbrace{\cos ^{-1} \left( -\widehat{\boldsymbol {h}}_{j} \cdot \widehat{\boldsymbol {h}}_{\rm{peak}} \right)}^{\,<\, 90^\circ }\right]^2}{\sum _j W_j}} \end{eqnarray*} (37) After this outlier rejection stage, we obtain the total angular momentum |$\boldsymbol {h}_{{{\rm{tot}}, \pm }}$| of the remaining particles in each peak, taking care of their statistical weights and using + to denote the corotating side (− for the counter-rotating side). We use this to refine our estimate of |$\widehat{\boldsymbol {h}}_{\rm{peak}}$|⁠, considering particles corresponding to both orbital pole peaks. \begin{eqnarray*} \widehat{\boldsymbol {h}}_{\rm{peak}} \ \propto \ \boldsymbol {h}_{{{\rm{tot}}, +}} - \boldsymbol {h}_{{{\rm{tot}}, -}}. \end{eqnarray*} (38) This leads to a revised measure of each peak’s angular dispersion. Using our newly estimated |${\widehat{\boldsymbol {h}}_{\rm{peak}}}$| and σ±, we repeat the σ-clipping procedure and thus update the peak locations and dispersions. This is continued iteratively until the algorithm converges (finds no more outliers), which usually takes  ∼ 10 steps. In this way, we obtain the directions of the corotating and counter-rotating peaks as well as their angular dispersions. We use |${\widehat{\boldsymbol {h}}_{\rm{peak}}}$| to increment the χ2 of our model by \begin{eqnarray*} \Delta \chi ^2 \ &=&\ \sum _{i = {\rm{MW,M31}} \chi ^2_{\rm{{SP}},i}}, \end{eqnarray*} (39) \begin{eqnarray*} \chi ^2_{{\rm{SP}},i} \ &=&\ \left[ \frac{\cos ^{-1} \left( \widehat{\boldsymbol {h}}_{\rm{peak}} \cdot \widehat{\boldsymbol {h}}_{{{\rm{SP}},i}} \right)}{\sigma _{{{\rm{SP}},i}}} \right]^2. \end{eqnarray*} (40) Here, |$\widehat{\boldsymbol {h}}_{{\rm{SP}},i}$| is the observed satellite plane spin vector of galaxy i (Table 7). This table also shows σSP, i, how well we expect our model to reproduce the observed satellite plane orbital pole of galaxy i. We include an allowance for both observational and modelling uncertainties. Compared to the case of the MW, the M31 satellite plane is observed to be thinner (Pawlowski et al. 2014, tables 2 and 4) and also has a smaller angular spread in most of our simulations. Thus, we allow a smaller angular uncertainty for the M31 satellite plane (i.e. σSP, M31 < σSP, MW). Table 7. Observed directions of the MW and M31 satellite plane angular momenta, shown in Galactic co-ordinates. Both point roughly in the direction of the Galactic anticentre. The MW satellite plane orientation is from Pawlowski & Kroupa (2013, section 3), while that of M31 is from their section 4. The error budgets used in our analysis (last column) allow for both observational and modelling uncertainties. Satellite plane spin vector of Direction Error Milky Way (MW) (176.4°, −15.0°) 25° Andromeda (M31) (206.2°, 7.8°) 15° Satellite plane spin vector of Direction Error Milky Way (MW) (176.4°, −15.0°) 25° Andromeda (M31) (206.2°, 7.8°) 15° Open in new tab Table 7. Observed directions of the MW and M31 satellite plane angular momenta, shown in Galactic co-ordinates. Both point roughly in the direction of the Galactic anticentre. The MW satellite plane orientation is from Pawlowski & Kroupa (2013, section 3), while that of M31 is from their section 4. The error budgets used in our analysis (last column) allow for both observational and modelling uncertainties. Satellite plane spin vector of Direction Error Milky Way (MW) (176.4°, −15.0°) 25° Andromeda (M31) (206.2°, 7.8°) 15° Satellite plane spin vector of Direction Error Milky Way (MW) (176.4°, −15.0°) 25° Andromeda (M31) (206.2°, 7.8°) 15° Open in new tab 4.2 Fraction of tidal debris that counter-rotate The method described in Section 4.1 yields the amount of material in the corotating and counter-rotating peaks, which we call W+ and W−, respectively. This can be used to obtain a counter-rotating fraction \begin{eqnarray*} f_{\rm {counter}} \ \equiv \ \frac{W_-}{W_+ \ +\ W_-}. \end{eqnarray*} (41) The values of fcounter can be compared with the proportion of material in each satellite plane on a counter-rotating orbit. To penalize models that have a different fcounter to what we expect, we add an appropriate contribution to the total χ2 of the model. \begin{eqnarray*} \Delta \chi ^2 \ &=&\ \sum _{i = \rm{MW,M31}} \chi ^2_{{\rm{CR}},i}, \end{eqnarray*} (42) \begin{eqnarray*} \chi ^2_{\rm{{CR}},i} \ &=&\ \left( \frac{f_{{\rm{counter}},i} - f_{{\rm{counter,obs}},i}}{\sigma _{f_{{\rm{counter,}}i}}} \right)^2. \end{eqnarray*} (43) Around M31, only |${\,{\sim }\,\frac{1}{2}}$| of its satellites lie in its satellite plane (Ibata et al. 2013). 13 of these 15 planar satellites share the same radial velocity trend (i.e. those on different sides of M31 have a different sign of radial velocity, once the systemic motion of M31 is subtracted). The remaining two satellites could well be interlopers or genuine counter-rotators. As there is no clear evidence for the latter, we assume that fcounter, obs, M31 = 0 but allow an uncertainty of |$\sigma _{f_{\rm{counter,M31}}} = 0.05$|⁠. This could be refined with proper motion measurements of M31 satellites, especially those which break the radial velocity trend (And XIII and And XXVII). For the MW, the observational picture is clearer due to the availability of proper motions for most of its classical satellites (Pawlowski et al. 2013, and references therein). Sculptor is counter-rotating within the MW satellite plane (Sohn et al. 2017). As it is one object out of ≈10 Galactic satellites with known proper motions, we require our models to yield fcounter, obs, MW = 0.10 ± 0.05. 4.3 MW thick disc 4.3.1 Age Our simulations of a past MW–M31 flyby must involve significant tidal torques on their discs in order to explain their significant misalignments with their satellite planes ( ∼ 50° for M31 and  ∼ 75° for the MW). As a result, the MW and M31 discs would have been dynamically heated. It is difficult to test this for M31 because it is an external galaxy viewed nearly face-on (Chemin et al. 2009). M31 has also experienced more recent interactions with other galaxies (Sadoun et al. 2014), significantly complicating any analysis. The merger history of the MW has been more quiescent (Wyse 2009). It is also much easier to get detailed ages, chemistry and kinematics for MW stars due to their proximity. Thus, one might expect our Galaxy to retain some evidence of a past interaction with M31. We suggest that its thick disc was formed in this event. Originally, the Galactic thick disc was identified based on the density of stars having a double-exponential law with distance out of the Galactic disc plane (Gilmore & Reid 1983). The lower mass thick disc component may well have arisen from dynamical heating of the MW thin disc during its close encounter with M31 (Zhao et al. 2013). Indeed, the thick disc seems to have formed fairly rapidly from its thin disc (Hayden et al. 2015). The associated burst of star formation (Snaith et al. 2014, fig. 2) could have further thickened the MW disc (Kroupa 2002). If this scenario is correct, then there ought to be a sudden jump in the age–velocity dispersion relation of nearby MW stars. This is in fact observed 9 ± 1 Gyr ago (Quillen & Garnett 2001), so the MW–M31 flyby almost certainly occurred then if at all. We consider this constraint on their trajectory by adding a contribution to the model χ2 of \begin{eqnarray*} \chi ^2_{\rm{TD, age}} \ =\ \left(\frac{t_f - t_{\rm{flyby} - 9\, Gyr}}{1 \,\rm{Gyr}} \right)^2. \end{eqnarray*} (44) 4.3.2 Thickness In addition to its age, another important property of the thick disc is its thickness or scale height. As our simulation lacks dissipation, it is most analogous to that part of the MW stellar disc which existed before the flyby. The gas disc would likely have settled down after the flyby, leading to a thin component to the MW stellar disc that would not be accurately included in our model. Thus, we expect to roughly match the observed scale height of the old MW thick disc. This was recently measured by Ferguson, Gardner & Yanny (2017), who fit its vertical density profile using \begin{eqnarray*} \rho \left( z \right) \ \propto \ \operatorname{sech}^2 \left( \frac{z}{2H_{2}}\right). \end{eqnarray*} (45) Based on their fig. 9, we take a value of H2 = 0.7 ± 0.1 kpc, though we inflate the uncertainty to 0.15 kpc for comparison with our simulations as these also have some uncertainty. This is partly because they lack the self-gravity required to obtain a |$\operatorname{sech}^2$| profile. To facilitate a comparison, we obtain the rms thickness of the observed and simulated thick discs. \begin{eqnarray*} z_{\rm{rms,obs}} \ =\ \frac{\pi }{\sqrt{3}}H_{2}. \end{eqnarray*} (46) The observed thickness zrms, obs must be at least partly due to processes internal to the MW such as interactions with stars, molecular clouds and spiral arms (e.g. Carlberg & Sellwood 1985; Kroupa 2002). Our simulations lack such processes, so any disc heating arises almost entirely due to the effect of M31. To account for this, we need to estimate how thick the Galactic disc might have been without such secular heating effects. We estimate this using Quillen & Garnett (2001), who showed that very young stars have a rather small vertical velocity dispersion σz. This is presumably because stars form out of dissipational gas which has settled itself into a rather thin disc (Nakanishi & Sofue 2016, fig. 6). Older stars have higher σz, but this quickly saturates at  ∼ 20 km s−1. We assume that this is what σz would have been without the M31 flyby. However, the thick disc has σz ≈ 40 km s−1 (Minchev et al. 2014, fig. 1). This suggests that internal processes were responsible for |${\,{\sim }\,\frac{1}{4}}$| of the disc heating required to explain the Galactic thick disc, with the rest arising from the close passage of M31 or the massive MW star clusters which formed in the associated starburst. Thus, our models should yield a thick disc which is thinner than observed. To account for this, we only require our models to reach \begin{eqnarray*} z_{\rm{TD}} \ =\ \frac{3}{4} z_{\rm{rms,obs}}. \end{eqnarray*} (47) This can be compared with the simulated value of zrms. Following Ferguson et al. (2017, fig. 12), we obtain this based on particles having a cylindrical Galactocentric radius within 1 kpc of the Solar circle, which we take to be 8.2 kpc from the Galactic Centre (McMillan 2017). We neglect effects due to non-axisymmetry, thus using the entire annulus corresponding to Galactocentric cylindrical radii of 7.2−9.2 kpc. This maximizes the number of particles used to determine the zrms of our simulated thick disc. Vertically, we allow particles up to 21.5 kpc either side of the MW disc, ensuring this is properly included. We then find zrms using a plane fitting procedure similar to that described in section 5.1 of Banik & Zhao (2018b) augmented by σ-clipping with a threshold of 2.58σ, corresponding to the 99 per cent confidence interval of a Gaussian distribution. We use as many σ-clipping stages as we need for the algorithm to converge, without reintroducing particles excluded at a previous stage. Usually, this process requires  ∼  10 iterations. Once it converges, we can determine |$\chi ^2_{\rm{TD, height}}$| based on the appropriately scaled uncertainty. This allows us to add the thick disc contribution to the total χ2 of each model. \begin{eqnarray*} \chi ^2_{\rm{TD}} \ =\ \chi ^2_{\rm{TD, age}} + \chi ^2_{\rm{TD, height}}. \end{eqnarray*} (48) 4.4 Alignment with the HVG plane So far, we have focused on how the MW and M31 may have affected each other. Their high relative velocity around the time of their flyby (668 km s−1, Fig. 4) would have enabled them to gravitationally slingshot any passing LG dwarf galaxy out at a high speed. Pawlowski & McGaugh (2014) previously identified several LG dwarfs whose radial velocities are indeed much higher than would be expected in ΛCDM. Similar results were subsequently obtained using more detailed dynamical models of the LG (Banik & Zhao 2016, 2017; Peebles 2017). These HVGs mostly lie very close to a well-defined plane which goes near both the MW and M31 (Banik & Zhao 2018b). Section 3 of that work shows how this naturally arises in an MOND simulation of the LG because the MW and M31 would be most efficient at scattering LG dwarfs parallel to their velocity. We therefore expect the HVG plane to coincide with the MW–M31 orbital plane. Figure 4. Open in new tabDownload slide Top: The MW–M31 separation d(t) in our best-fitting model, using the cosmological parameters in Table 1. The parameters of this trajectory are given in Table 9. We start integrating test particle trajectories 2 Gyr after the big bang (Section 3). Bottom: The MW–M31 relative velocity, which is very high at early times due to cosmic expansion. The discontinuity at 2.6 Gyr (first turnaround) arises because we assume the orbit is radial prior to this point but has its present angular momentum at later times (Section 2.4). The MW–M31 flyby occurs at 6.17 Gyr with a relative velocity of 668.1 km s−1. Notice the very slow relative velocity at the second apocentre – even a tiny perturbation can significantly reorient the MW–M31 orbit. Thus, the present direction of |$\widehat{\boldsymbol {h}}_{{LG}}$| is not a good guide to what |$\widehat{\boldsymbol {h}}_{{LG}}$| was several Gyr ago. Figure 4. Open in new tabDownload slide Top: The MW–M31 separation d(t) in our best-fitting model, using the cosmological parameters in Table 1. The parameters of this trajectory are given in Table 9. We start integrating test particle trajectories 2 Gyr after the big bang (Section 3). Bottom: The MW–M31 relative velocity, which is very high at early times due to cosmic expansion. The discontinuity at 2.6 Gyr (first turnaround) arises because we assume the orbit is radial prior to this point but has its present angular momentum at later times (Section 2.4). The MW–M31 flyby occurs at 6.17 Gyr with a relative velocity of 668.1 km s−1. Notice the very slow relative velocity at the second apocentre – even a tiny perturbation can significantly reorient the MW–M31 orbit. Thus, the present direction of |$\widehat{\boldsymbol {h}}_{{LG}}$| is not a good guide to what |$\widehat{\boldsymbol {h}}_{{LG}}$| was several Gyr ago. Unfortunately, there are some serious difficulties with testing if this is indeed the case. The main problem is the very small number of independent HVGs and thus the uncertain orientation of what plane they prefer. Most of the HVGs we identified belong to the NGC 3109 association. In ΛCDM, it is unlikely to have been gravitationally bound as this would require a rather large and massive dark matter halo (Bellazzini et al. 2013; Kourkchi & Tully 2017). Given the high radial velocities of these galaxies, they need to have been rather close to the MW at some point in the past. Dynamical friction with its dark matter halo probably rules out this scenario, which is why we treated these galaxies as independent HVGs in Banik & Zhao (2018b). However, this consideration would not arise in MOND as dynamical friction would only operate on baryons. Thus, the NGC 3109 association might have been gravitationally bound and still have passed rather close to the MW without merging. Indeed, this seems rather likely given the filamentary nature of the association (Bellazzini et al. 2013). It may be an unbound structure somewhat similar to a tidal stream, defined by dwarf galaxies instead of stars. In this scenario, we would have only a very small number of independent HVGs with which to reliably determine the HVG plane normal. The angular dispersion of HVGs with respect to the MW–M31 orbital plane would be ∼cos −10.9 ≈ 25° (Banik & Zhao 2018b, fig. 7).7 An uncertainty of this order seems reasonable given the 16° angle of the MW–M31 line out of the HVG plane. As there are not very many clearly independent HVGs and our simulations also have some uncertainty, we allow an angular mismatch of σHVG = 25° between the HVG plane normal |$\widehat{\boldsymbol {n}}_{\rm{HVG}}$| and the MW–M31 orbital pole |$\widehat{\boldsymbol {h}}_{\rm{LG}}$|⁠. Also allowing for the unknown sense in which the HVGs rotate within their preferred plane, we thus add an extra contribution to the model χ2 of \begin{eqnarray*} \chi ^2_{\rm{HVG}} \ =\ \left( \frac{\cos ^{-1} \left| \widehat{\boldsymbol {n}}_{\rm{HVG}} \cdot \widehat{\boldsymbol {h}}_{\rm{LG}} \right|}{\sigma _{\rm{HVG}}}\right)^2. \end{eqnarray*} (49) 5 RESULTS AND DISCUSSION To see if MOND can explain the observed properties of the LG, we run a grid of models in the present values of the EF on the LG, MW–M31 orbital pole, and Galactocentric tangential speed of M31. We then obtain a χ2 statistic for each model (equation 35) and thereby select the best-fitting model. We improve its accuracy by using the Newton–Raphson method to adjust the initial orientations of the MW and M31 discs so as to best match their presently observed orientations (Section 3.1). We then use this in a new grid of models, focusing on the most promising region of parameter space in the previous grid. After alternating between the Newton–Raphson and grid searches a few times, we improve the resolution in parameter space to several times better than the expected modelling and observational uncertainties. Our best-fitting model is compared with LG observations in Table 8 and discussed next. Table 8. Comparison of our best-fitting model with observational constraints on the LG. The uncertainties σ are shown in the central column e.g. σHVG = 25°. 0.0045 of the mass in the MW ends up in its satellite plane, including both the corotating and counter-rotating peaks in the orbital pole distribution of material in its satellite region (Table 6). For M31, the equivalent quantity is 6.5 × 10−6. As discussed in Section 4.4, we expect the MW–M31 orbital plane to align with the plane of high-velocity galaxies discovered by Banik & Zhao (2018b). Open in new tab Table 8. Comparison of our best-fitting model with observational constraints on the LG. The uncertainties σ are shown in the central column e.g. σHVG = 25°. 0.0045 of the mass in the MW ends up in its satellite plane, including both the corotating and counter-rotating peaks in the orbital pole distribution of material in its satellite region (Table 6). For M31, the equivalent quantity is 6.5 × 10−6. As discussed in Section 4.4, we expect the MW–M31 orbital plane to align with the plane of high-velocity galaxies discovered by Banik & Zhao (2018b). Open in new tab The total χ2 of our best model is 14.08. There are three model parameters and seven contributions to χ2. However, the orientation of a plane is a 2D concept such that the satellite plane and HVG constraints should really count for six constraints rather than three. In this case, there are 10 constraints – 3 model parameters = 7 ‘degrees of freedom’, for which the chance of an even higher χ2 is |$P = 4.98\hbox{ per cent}$| based on the chi2cdf function of matlab®. However, even if we count the satellite and HVG planes as one constraint each, this still leaves 4 degrees of freedom for which |$P = 0.7\hbox{ per cent}$|⁠. Given the crudeness of our model, we believe that it matches a wide array of LG observations well enough to form the basis for more detailed modelling (e.g. Bílek et al. 2017). To get an idea of how the different constraints affect the parameters of our best-fitting model, we repeat our χ2 analysis neglecting one contribution to χ2 at a time. The best-fitting model is unaffected by the constraints from the thick disc age, the thickness of the MW disc, the orientation of the HVG plane found by Banik & Zhao (2018b) and the fraction of counter-rotating material around M31 (presumably because our models had no such material). Neglecting the constraint from the MW satellite plane orientation, the preferred model has a slightly smaller closest approach distance of 17.90 kpc. Ignoring the counter-rotating fraction of tidal debris around the MW, the best-fit orbital pole tilts 5° further south and the model prefers a weaker present gext of 0.02a0. Perhaps the most important single constraint is the accurately known M31 satellite plane orientation. Without this information, the best-fitting model prefers an orbital pole 15° further north and a smaller closest approach distance of 17.90 kpc. Even these changes are very small, suggesting that no individual observational constraint greatly affects our analysis. 5.1 Implied MW–M31 trajectory 5.1.1 The timing argument Fig. 4 shows the MW–M31 trajectory in our best-fitting model, with the main model parameters summarized in Table 9. The flyby occurs 7.65 Gyr ago, consistent with the thick disc age of 9 ± 1 Gyr (Quillen & Garnett 2001). To satisfy the timing argument (Section 2.5), our model requires the MW and M31 to have a total mass of 2.81 × 1011 M⊙. This is slightly higher than our estimate of the total baryonic mass in their discs, assuming the MW rotation curve flatlines at vf = 180 km s−1 (Kafle et al. 2012) while that of M31 flattens at 225 km s−1 (Carignan et al. 2006). Taken at face value, these rotation speeds imply a total mass of 2.3 × 1011 M⊙ because the MOND dynamical mass of a galaxy is \begin{eqnarray*} M_{\rm{dyn}} \ \ =\ \ \frac{{v_{f}}^4}{Ga_{0}}. \end{eqnarray*} (50) Table 9. Data concerning the MW–M31 trajectory in our best-fitting model, with the present separation obtained from McConnachie (2012). Our adopted cosmology (Table 1) implies the Universe is 13.82 Gyr old, so the last MW–M31 pericentre was 7.65 Gyr ago. The direction of the EF on the LG is obtained from Kogut et al. (1993) and is the same in all models (Section 2.2). The model has a peak EF strength of 0.029a0 when the cosmic scale factor a = 0.75. Although the simulated Galactocentric tangential velocity of M31 is currently in almost the opposite direction to that suggested by van der Marel et al. (2012), it is consistent within uncertainties as the orbit is nearly radial (see their fig. 3). The grid of models we run has a step of 0.001a0 in the present EF on the LG, 1.15 km s−1 in the tangential velocity of M31 (≈1.7 kpc in closest approach distance), and 5° in the MW–M31 orbital pole. Quantity Value Orbital pole of MW–M31 orbit (24.6°, −16.3°) MW–M31 separation 783 kpc MW–M31 tangential velocity 16.70 km s−1 MW+M31 mass 2.81 × 1011 M⊙ Flyby time (after big bang) 6.17 Gyr Closest approach distance 19.6 kpc Present external field on LG 0.022 a0 External field direction (276°, 30°) Quantity Value Orbital pole of MW–M31 orbit (24.6°, −16.3°) MW–M31 separation 783 kpc MW–M31 tangential velocity 16.70 km s−1 MW+M31 mass 2.81 × 1011 M⊙ Flyby time (after big bang) 6.17 Gyr Closest approach distance 19.6 kpc Present external field on LG 0.022 a0 External field direction (276°, 30°) Open in new tab Table 9. Data concerning the MW–M31 trajectory in our best-fitting model, with the present separation obtained from McConnachie (2012). Our adopted cosmology (Table 1) implies the Universe is 13.82 Gyr old, so the last MW–M31 pericentre was 7.65 Gyr ago. The direction of the EF on the LG is obtained from Kogut et al. (1993) and is the same in all models (Section 2.2). The model has a peak EF strength of 0.029a0 when the cosmic scale factor a = 0.75. Although the simulated Galactocentric tangential velocity of M31 is currently in almost the opposite direction to that suggested by van der Marel et al. (2012), it is consistent within uncertainties as the orbit is nearly radial (see their fig. 3). The grid of models we run has a step of 0.001a0 in the present EF on the LG, 1.15 km s−1 in the tangential velocity of M31 (≈1.7 kpc in closest approach distance), and 5° in the MW–M31 orbital pole. Quantity Value Orbital pole of MW–M31 orbit (24.6°, −16.3°) MW–M31 separation 783 kpc MW–M31 tangential velocity 16.70 km s−1 MW+M31 mass 2.81 × 1011 M⊙ Flyby time (after big bang) 6.17 Gyr Closest approach distance 19.6 kpc Present external field on LG 0.022 a0 External field direction (276°, 30°) Quantity Value Orbital pole of MW–M31 orbit (24.6°, −16.3°) MW–M31 separation 783 kpc MW–M31 tangential velocity 16.70 km s−1 MW+M31 mass 2.81 × 1011 M⊙ Flyby time (after big bang) 6.17 Gyr Closest approach distance 19.6 kpc Present external field on LG 0.022 a0 External field direction (276°, 30°) Open in new tab When it comes to the MW–M31 timing argument, their large separation requires us to consider any massive satellites as these also contribute to the MW–M31 attraction. The LMC has vf ≈ 70 km s−1 (Alves & Nelson 2000; van der Marel & Kallivayalil 2014), while vf ≈ 120 km s−1 for M33 (Corbelli 2003; Corbelli & Salucci 2007). Thus, the most luminous LG satellites only contribute a mass of ≈1.5 × 1010 M⊙ without much affecting the MW:M31 mass ratio. This can only explain |${\approx \frac{1}{3}}$| of the discrepancy between the timing argument mass of the LG and that inferred from rotation curves of LG objects. The remaining discrepancy disappears if actual rotation speeds are 3 per cent higher than we assume. This is quite possible given expected observational uncertainties for the dominant LG galaxy M31 (e.g. Chemin et al. 2009, table 4). The high timing argument mass might also be a sign of hot gas haloes surrounding the MW and M31. Nicastro et al. (2016) conducted observations at a range of Galactic latitudes and found that the best-fitting model has a hot gas corona around the MW with a mass of 2 × 1010 M⊙ (see their table 2 model A). This is consistent with the estimate of 2.7 ± 1.4 × 1010 M⊙ obtained by Salem et al. (2015) based on how the LMC gas disc is truncated at a smaller radius than its stellar disc, a phenomenon which they attribute to ram pressure stripping by the MW corona. If this was much more massive, the MOND-predicted MW escape velocity curve would be higher and flatter than inferred from observations (Banik & Zhao 2018a, fig. 5). A very massive Galactic hot gas corona would also reduce the apocentre of the Sagittarius tidal stream (Thomas et al. 2017, section 4.3), making it difficult to explain the large measured distances to some parts of the stream (Belokurov et al. 2014). A corona has also been detected around M31 based on absorption features in spectra of background quasars (Lehner, Howk & Wakker 2015). Thus, a small amount of hot gas around the MW and M31 could easily explain why their total MOND timing argument mass slightly exceeds what is implied by their rotation curves and those of their most luminous satellites. 5.1.2 The apocentre asymmetry An interesting feature of our MW–M31 trajectory is that their second apocentre is larger than their first. This apocentre asymmetry is very important because it is closely related to the timing of the MW–M31 flyby. Naively, one might expect that the MW–M31 system is currently near apocentre so the flyby must have been relatively recent. For example, if we take the most recent apocentre to be 12 Gyr after the big bang and consider this to be |$\frac{3}{2}$| orbital periods, then the flyby occurred after 1 orbital period i.e. at a time of 8 Gyr, only  ∼ 6 Gyr ago. This is much more recent than the observed 9 ± 1 Gyr age of the Galactic thick disc which we expect formed due to the flyby. To match this constraint, our model will try to find some way of getting a small first apocentre so the MW and M31 can quickly turn around and undergo their flyby. Afterwards, some extra effect is needed to push the galaxies further out, thus matching their present separation. Part of the necessary repulsion must come from the cosmological acceleration term in equation (30). However, this would cause only a small effect because the MW–M31 dynamics is largely decoupled from cosmic expansion. Indeed, a quick look at the EF-included trajectory in fig. 1 of Zhao et al. (2013) shows that the two apocentres are expected to be very similar despite those authors including the cosmological acceleration. Although the apocentre asymmetry could be due to our assumed EF history (Section 2.2), this is not the whole story because the asymmetry is much weaker if |$\widehat{\boldsymbol {h}}_{\rm{LG}}$| is changed while not adjusting |$\boldsymbol {g}_{\rm{ext}}$|⁠. Instead, our investigations reveal that the apocentre asymmetry is mainly due to tides raised by perturbers on the LG (Section 2.3). To investigate this possibility in more detail, we perform some analytic calculations of how much tides might be expected to influence the MW–M31 trajectory. We treat the MW and M31 as test particles in the EF-dominated gravitational field due to each perturber, whose tidal stress tensor is found in Appendix B. The EF-dominated assumption is likely to hold true at the first MW–M31 apocentre and should be very accurate at their second apocentre. Table 10 gives our analytic estimates for how much tides from each perturber would have affected the MW–M31 relative acceleration at the time of their first apocentre and at the present time. For reference, the present MW–M31 attraction may be estimated as |$\frac{\sqrt{GM_{\rm{LG}} a_{0}}}{d_{\rm{M31}}} = 2.77 \times 10^{-12}$| m s−2, though this would be weakened |${\approx 30\hbox{ per cent}}$| due to the EF (Fig. 1). It is evident that the perturbers we consider would likely have created extra attraction between the MW and M31 in the period before their flyby, thus limiting the distance at which they turned around and causing them to have an earlier flyby. However, the perturbers would not have continued ‘compressing’ the MW–M31 orbit after the flyby. This is no doubt partly because of cosmic expansion, which affects the LG-perturber distances far more than it affects the much more bound MW–M31 system. Table 10. Our estimates for how much each perturber affects the MW–M31 relative acceleration at the times shown. We obtain these estimates by applying the analytic method described in Appendix B. For the GA, we use equation (29) instead. The last row shows the contribution of the cosmological acceleration term. The MW and M31 are expected to have a relative acceleration of ≈2 × 10−12 m s−2once the external field is considered (Section 5.1.2). Perturber Estimated gtide at Estimated first apocentre (m s−2) gtide now (m s−2) Cen A −1.85 × 10−13 +6.22 × 10−14 M81 −6.80 × 10−14 −4.11 × 10−15 IC 342 −5.33 × 10−13 +9.75 × 10−15 GA −8.34 × 10−15 +1.64 × 10−14 |$\frac{\ddot{a}}{a} d_{\rm{M31}}$| −6.47 × 10−13 +6.06 × 10−14 Perturber Estimated gtide at Estimated first apocentre (m s−2) gtide now (m s−2) Cen A −1.85 × 10−13 +6.22 × 10−14 M81 −6.80 × 10−14 −4.11 × 10−15 IC 342 −5.33 × 10−13 +9.75 × 10−15 GA −8.34 × 10−15 +1.64 × 10−14 |$\frac{\ddot{a}}{a} d_{\rm{M31}}$| −6.47 × 10−13 +6.06 × 10−14 Open in new tab Table 10. Our estimates for how much each perturber affects the MW–M31 relative acceleration at the times shown. We obtain these estimates by applying the analytic method described in Appendix B. For the GA, we use equation (29) instead. The last row shows the contribution of the cosmological acceleration term. The MW and M31 are expected to have a relative acceleration of ≈2 × 10−12 m s−2once the external field is considered (Section 5.1.2). Perturber Estimated gtide at Estimated first apocentre (m s−2) gtide now (m s−2) Cen A −1.85 × 10−13 +6.22 × 10−14 M81 −6.80 × 10−14 −4.11 × 10−15 IC 342 −5.33 × 10−13 +9.75 × 10−15 GA −8.34 × 10−15 +1.64 × 10−14 |$\frac{\ddot{a}}{a} d_{\rm{M31}}$| −6.47 × 10−13 +6.06 × 10−14 Perturber Estimated gtide at Estimated first apocentre (m s−2) gtide now (m s−2) Cen A −1.85 × 10−13 +6.22 × 10−14 M81 −6.80 × 10−14 −4.11 × 10−15 IC 342 −5.33 × 10−13 +9.75 × 10−15 GA −8.34 × 10−15 +1.64 × 10−14 |$\frac{\ddot{a}}{a} d_{\rm{M31}}$| −6.47 × 10−13 +6.06 × 10−14 Open in new tab However, cosmic expansion cannot be the entire explanation because the tidal forcing changes sign for Cen A and IC 342, something that expansion alone cannot do. The sign change is an indication that the geometry of the problem is different after the flyby than before it. This is unsurprising given that the MW–M31 line must rotate substantially around the time of their close encounter. In fact, our model has |$\widehat{\boldsymbol {d}}_{\rm{M31}}$| rotate 233° between the first MW–M31 apocentre and today, significantly changing the relative orientations of the MW, M31, and external perturbers. Combined with other effects, this explains the apocentre asymmetry in our best model and thus how the MW–M31 flyby could be as long ago as the formation of the Galactic thick disc. 5.1.3 MW–M31 orbital pole Our best model has a very different MW–M31 orbital pole to that indicated by the toy model of Banik & Zhao (2018b, section 2.2). That model relied on many assumptions, almost certainly oversimplifying the problem. Of particular concern is the impulse approximation, which is not likely to work, given the long duration of the flyby compared to the dynamical time of MW stars. For example, the MW–M31 separation is below twice its minimum value for  ∼ 120 Myr (Fig. 4), a significant fraction of the time taken for the Sun to orbit around the MW (McMillan 2017). It is also unrealistic to assume that M31 only significantly affects the MW at the time of closest approach as M31 gets within 20 kpc of the MW (Table 9). In any case, M31 would come closest to different parts of the MW at different times, an effect not included in the toy model. Moreover, the model does not contain any rigorous dynamical arguments and is a purely kinematic analysis. It is therefore quite possible that the toy model did not adequately capture how the MW–M31 flyby would have worked. The main goal of this contribution is to construct a more rigorous dynamical model of the flyby as a stepping stone to N-body simulations. An important constraint on our model is the present relative velocity between the MW and M31. We determine this by adjusting the observed value (van der Marel et al. 2012) for the motion of the LMC (Kallivayalil et al. 2013) and M33 (Brunthaler et al. 2005) as well as the motion of the Sun within the MW (McMillan 2011; Francis & Anderson 2014). Applying equation (50) suggests that the LMC has 2.25 per cent as much mass as the MW (Alves & Nelson 2000; Kafle et al. 2012), while M33 has 7.5 per cent as much mass as M31 (Corbelli 2003; Carignan et al. 2006). Combining this information yields the relative velocity between the MW–LMC and M31–M33 barycentres. We fix the radial component of this velocity at –93.4 km s−1 but let the tangential component vary so as to alter the MW–M31 orbital pole |$\widehat{\boldsymbol {h}}_{\rm{LG}}$| and their closest approach distance. In our best-fitting model, the tangential velocity of M31 is similar in magnitude to the observational estimate of van der Marel et al. (2012) but points in almost the opposite direction. Their fig. 3 shows that this is consistent within uncertainties because the MW–M31 orbit is nearly radial. Future high-precision proper motion measurements of M31 could provide some constraints on |$\widehat{\boldsymbol {h}}_{\rm{LG}}$| and thus whether it might be the same as in our best-fitting model. However, this requires a good idea of the LMC and M33 masses as one would need to consider the relative velocity between the MW–LMC barycentre and the M31–M33 barycentre. Both M33 and the LMC have velocities of ≈200 km s−1 orthogonal to the MW–M31 line, affecting their tangential velocity by |${\approx } \left( \frac{120}{225} \right)^4\times$| as much i.e. by 16 km s−1, even if we only consider the more massive M33. If we assume uncertainties of 10 per cent on its rotation curve, this implies a 6 km s−1 uncertainty on the present tangential velocity of M31. In addition to this, M33 also affects the position of M31 as the two orbit around their barycentre. Given a M31–M33 distance of ≈200 kpc (Patel, Besla & Sohn 2017, table 1), the reflex motion of M31 must be ≈16 kpc. This is a significant fraction of our estimated 20 kpc distance between the MW and M31 at the time of their closest approach (Table 9). To avoid this becoming a major source of uncertainty, we would need to know how the M31–M33 orbit behaved 8 Gyr into the past. Given the nearly radial MW–M31 orbit, determining |$\widehat{\boldsymbol {h}}_{\rm{LG}}$| also requires accurate knowledge of the motion of the Sun with respect to the MW, in particular the velocity of the Local Standard of Rest (McMillan 2017). Although these considerations have only a small effect on the MW–M31 orbit over the last few Gyr, it is so nearly radial that such effects could determine whether M31 flew past the MW on one side or on the other. This makes it very difficult to directly determine the geometry of the MW–M31 flyby through backward integration of present conditions. Fortunately, some insights can be gained from our MOND timing argument analysis. This indicates that the LG must be unrealistically massive if the MW–M31 orbital pole lies within  ∼ 75° of Galactic co-ordinates (227°, −34°). Other orbital poles can satisfy the timing argument with a LG mass more consistent with the observed baryonic masses of its constituents. As this still leaves a wide range of possibilities, we focus on other ways of testing our model. 5.2 The Local Group satellite planes 5.2.1 Orientation A strong test of our models is provided by the observed orientations of the MW and M31 satellite planes (Table 7). Although these may have precessed slightly due to the non-spherical nature of the MW and M31, significant precession is unlikely as all satellites are not at the same radial distance (Pawlowski & Kroupa 2013, fig. 2). Thus, different satellites would precess at different rates, thickening the satellite planes (Fernando et al. 2017). Their small observed thickness therefore strongly suggests that their orientations have not substantially changed since they were formed. To see how well we can fit these constraints, we show the orbital pole distribution of all satellite particles around the MW and M31 (Fig. 5). Both orbital pole distributions show a lack of material on a similar orbit to the disc of the host galaxy, a consequence of the filters we applied to our simulation in order to focus on the MW and M31 satellite regions (Table 6). Figure 5. Open in new tabDownload slide The distribution of orbital angular momentum directions (spin vectors) for tidal debris around the MW and M31 discs at the end of our best-fitting simulation (parameters in Table 9). Satellite particles are selected according to the criteria in Table 6, thus rejecting the MW and M31 discs themselves. The mass scale used is arbitrary (satellite plane masses are discussed in Section 5.2.3). For model comparison, the method described in Section 4.1 is applied to a low-resolution version of these histograms prior to further refinement. Top: Results for the MW. Its disc spin vector points at the South Galactic Pole, while the open pink circle shows the spin vector of its satellite plane (Table 7). Bottom: Results for M31. We use open pink circles to show the observed spin vectors of its disc (lower left) and satellite plane (upper right). Figure 5. Open in new tabDownload slide The distribution of orbital angular momentum directions (spin vectors) for tidal debris around the MW and M31 discs at the end of our best-fitting simulation (parameters in Table 9). Satellite particles are selected according to the criteria in Table 6, thus rejecting the MW and M31 discs themselves. The mass scale used is arbitrary (satellite plane masses are discussed in Section 5.2.3). For model comparison, the method described in Section 4.1 is applied to a low-resolution version of these histograms prior to further refinement. Top: Results for the MW. Its disc spin vector points at the South Galactic Pole, while the open pink circle shows the spin vector of its satellite plane (Table 7). Bottom: Results for M31. We use open pink circles to show the observed spin vectors of its disc (lower left) and satellite plane (upper right). For the MW, there is a broad peak in the orbital pole distribution near a Galactic longitude of 180°, similar to the observed spin vector of the MW satellite plane. The Galactic latitude of the simulated satellite plane is  ∼ 25° lower than observed, indicating that there is not quite enough tidal torque on the MW in our model. However, a discrepancy at this level is reasonable giving expected observational and modelling errors. For example, treating M31 as a disc rather than a point mass could allow parts of it to get much closer to the MW, thus exerting a larger pull on it. Moreover, the orbital poles of MW satellites have an angular dispersion of ≈25° (Pawlowski & Kroupa 2013, section 4). The simulated M31 satellite particles prefer a much narrower range of orbital poles – our orbital pole analysis (Section 4.1) yields an angular dispersion of 41° for the MW satellite plane but only 23° for M31. This may explain why the MW satellite plane indeed has nearly twice the thickness of the M31 satellite plane (Pawlowski et al. 2014, tables 2 and 4). However, both structures are observed to be much thinner than in our best-fitting model. This is to be expected given its lack of dissipative processes. Due to reduced observational and modelling uncertainties in its orientation, the M31 satellite plane provides a strong test of our models. Thus, it is interesting that our best-fitting model matches the observed spin vector of this structure to within 10° (Table 8). 5.2.2 Counter-rotating fraction Around the MW, a secondary orbital pole peak is also evident near a Galactic longitude of 350° (Fig. 5). This is not quite where a true counter-rotating peak would be, though it is in roughly the right place. Due to dissipation, we expect that the actual tidal debris around the MW would settle into a plane with the counter-rotating material also forced to orbit within this plane but in the opposite sense. Thus, our model strongly suggests that some MW satellites ought to be counter-rotating. Indeed, Sculptor is just such an object (Sohn et al. 2017). Taken at face value, our model predicts that 22 per cent of the MW satellite plane mass ought to be in similar objects, much more than observed. If there is a significant amount of dissipation before the tidal debris coalesce into a few satellites, then the counter-rotating material could fall to much lower radii and perhaps get re-accreted on to the MW. Thus, the counter-rotating fraction in our model only provides an upper limit on how much counter-rotating material might actually be present in the MW satellite plane. Unlike the MW, our model yields no counter-rotating material around M31. Although it lacks dissipation, it is very difficult to see how this could explain future observations that reveal several counter-rotating satellites within the M31 satellite plane. Thus, our flyby model may be tested with proper motions of the 2 M31 satellites in its satellite plane that do not share the coherent radial velocity trend of the remaining 13 M31 satellites in this structure (Ibata et al. 2013). If both of these anomalous satellites (And XIII and And XXVII) have proper motions indicating orbits within the M31 satellite plane but opposite to the majority, then this implies a significant counter-rotating fraction, thus challenging our model. However, if these satellites have orbits that are not well aligned with the M31 satellite plane, then by definition they are just interlopers. In this situation, the planar M31 satellites would all be corotating, much more consistent with our model. To be sure of this conclusion, one would need to obtain proper motions for any subsequently discovered members of the M31 satellite plane in addition to its already known members. For example, Cassiopeia III (Martin et al. 2013) is another example of a M31 satellite spatially within the M31 satellite plane but whose radial velocity indicates that it is counter-rotating with respect to this structure (Martin et al. 2014). 5.2.3 Mass As well as explaining their orientations, our simulations should broadly match the observed masses of the LG satellite planes. In principle, masses can be estimated from luminosities as galaxies lack dark matter haloes in MOND. However, the MW and M31 satellite planes consist of rather faint galaxies, many of which were discovered rather recently. Due to the large associated uncertainties, we do not use satellite plane masses to select our best-fitting model. We can therefore use this as an independent check on whichever model best fits the other constraints (Section 4). In this model, 0.0045 of the MW mass ends up in its satellite plane, including both the corotating and counter-rotating peaks in the orbital pole distribution. We can estimate the actual mass of the MW satellite plane by considering its most massive satellite, the LMC. If we assume that it has vf = 70 km s−1 (Alves & Nelson 2000) while the MW has vf = 180 km s−1 (Kafle et al. 2012), then we expect the satellite plane mass fraction to be |$\,{\sim }\,\left( \frac{70}{180} \right)^4 = 0.02$| (equation 50). For M31, this quantity is only 6.5 × 10−6 in our model. Observationally, M31 satellite plane galaxies are each ≈4−5 orders of magnitude fainter than M31 itself (Martin et al. 2016, table 2). Although M33 is much brighter, it must have an independent origin as it is not part of the M31 satellite plane (Ibata et al. 2013, fig. 1). If this structure also contains M32, this would make the M31 satellite plane much more massive (McConnachie 2012, table 4). However, it is unclear whether M32 is really part of this structure due to the very small distance between M32 and M31 (Weisz et al. 2014, table 4) – any M31 satellite this close in would lie within the M31 satellite plane. Moreover, M32 has likely been accreted recently by M31 (Dierickx, Blecha & Loeb 2014). Although that work considered ΛCDM, this is also possible in MOND because M32 would experience some dynamical friction against the baryonic disc of M31 if the galaxies had a sufficiently close encounter. As this process would eventually cause M32 to merge with M31, the current detectability of M32 and M31 as separate objects makes it unlikely that M32 has maintained its present tight orbit around M31 over the ≈8 Gyr since the MW–M31 flyby. Thus, we expect that M32 is unrelated to the M31 satellite plane. Our model correctly predicts that the M31 satellite plane should be ≈3 orders of magnitude fainter than the MW satellite plane, though both satellite plane masses are underestimated by about an order of magnitude. The lower mass of the M31 satellite plane is almost certainly related to the MW having a lower mass than M31. This makes it more difficult for the MW to tidally disrupt M31 than vice versa. 5.2.4 Initial radial distribution Our low predicted satellite plane masses may be related to our assumed exponential disc law, which could well be inaccurate in the outer parts of each galaxy that go on to form their satellite planes. The importance of this can be seen in Figs 6 and 7, whose left-hand sides (left of vertical green lines) show the initial radial distribution of material ending up in the MW and M31 satellite planes, respectively. It is evident that a significant fraction of the MW satellite particles were initially ≈50 kpc from the MW, even though its stellar disc has a scale length of only 2.15 kpc. Thus, our simulated satellite plane masses are sensitive to assumptions concerning the MW disc surface density at large radii. Partly to deal with this, we also include a gas disc with a longer scale length of 7 kpc (Table 3). Thus, we have to consider test particles starting at quite large radii. Our adopted outer limit of 107.5 kpc (Table 4) seems to be sufficient for our investigation as the initial radial distribution of MW satellite plane particles tails off by this point. Figure 6. Open in new tabDownload slide Histogram of the initial and final Galactocentric distances of MW particles in its satellite region (defined in Section 4.1). Each particle contributes the mass it represents within the MW disc (Section 3). The particles are divided according to whether they are currently in the corotating peak (top) or the counter-rotating peak (bottom). The initial radial distribution of the particles is shown to the left of the vertical green lines, while the parts of the figures to the right of these lines show the present-day distribution of the same material. The dashed horizontal lines show distances to the indicated MW satellite. The proper motion of Leo II indicates it is in the MW satellite plane (Piatek, Pryor & Olszewski 2016), but most of this structure is defined by less distant satellites (Pawlowski & Kroupa 2013, fig. 2). The only known counter-rotating MW satellite is Sculptor (Sohn et al. 2017). The observed radial extents of both satellite planes are nicely illustrated in Pawlowski (2018, fig. 1). Figure 6. Open in new tabDownload slide Histogram of the initial and final Galactocentric distances of MW particles in its satellite region (defined in Section 4.1). Each particle contributes the mass it represents within the MW disc (Section 3). The particles are divided according to whether they are currently in the corotating peak (top) or the counter-rotating peak (bottom). The initial radial distribution of the particles is shown to the left of the vertical green lines, while the parts of the figures to the right of these lines show the present-day distribution of the same material. The dashed horizontal lines show distances to the indicated MW satellite. The proper motion of Leo II indicates it is in the MW satellite plane (Piatek, Pryor & Olszewski 2016), but most of this structure is defined by less distant satellites (Pawlowski & Kroupa 2013, fig. 2). The only known counter-rotating MW satellite is Sculptor (Sohn et al. 2017). The observed radial extents of both satellite planes are nicely illustrated in Pawlowski (2018, fig. 1). Figure 7. Open in new tabDownload slide Similar to Fig. 6, but for M31. The dashed horizontal line at 192 kpc is the rms radial extent of the observed M31 satellite plane (Ibata et al. 2014, section 3). Our model lacks counter-rotating particles in the M31 satellite region (Fig. 5). Figure 7. Open in new tabDownload slide Similar to Fig. 6, but for M31. The dashed horizontal line at 192 kpc is the rms radial extent of the observed M31 satellite plane (Ibata et al. 2014, section 3). Our model lacks counter-rotating particles in the M31 satellite region (Fig. 5). For M31, we use a similar initial radial range – our simulations start with M31 particles out to 106 kpc from M31 (20 disc scale lengths, Courteau et al. 2011). This seems to be enough for our purposes (Fig. 7). However, this is such a large distance that our estimated M31 satellite plane mass is rather sensitive to its disc surface density at large radii. Even if M31 can be treated as an exponential disc, there is a |${10\hbox{ per cent}}$| uncertainty in its scale length (Courteau et al. 2011). To estimate how this might affect our results, we consider the M31 mass fraction exterior to 60 kpc, the approximate initial radius of most M31 satellite plane particles (Fig. 7). This fraction would be 2.4 × higher if we assume the M31 disc scale length is 5.8 kpc instead of 5.3 kpc. Thus, it is unlikely that our model will provide a very good match to the masses of the LG satellite planes, though it should yield the correct order of magnitude. Although these uncertainties clearly affect the mass of the LG satellite planes, their orientations should be less sensitive – we weighted all MW and M31 particles equally and found very little change to the appearance of Fig. 5. This is because there is a clear peak in the orbital pole distribution for each galaxy such that its location is not very sensitive to details. For the MW, there is a clear secondary orbital pole peak which is approximately counter-rotating relative to the primary peak (Fig. 5). The particles in the secondary peak start at similar radii to the corotating particles (Fig. 6). This result is different to that obtained by Pawlowski et al. (2011), whose Fig. 8 indicates that the counter-rotating particles start further out. There are many possible reasons for this difference, e.g. using a different mass ratio, geometry, and gravity law. Figure 8. Open in new tabDownload slide Left: Histogram of |$\widehat{\boldsymbol {r}} \cdot \widehat{\boldsymbol {v}}$|⁠, where |$\boldsymbol {r}$| and |$\boldsymbol {v}$| are, respectively, the position and velocity of each particle relative to the MW. This is a basic measure of orbital eccentricity. We only show particles corresponding to the MW satellite plane according to the criteria in Table 6. The dashed blue line shows the expected distribution for particles on random orbits around the MW. Right: Analogous results for M31. All particles shown are in the corotating peak around M31 as there are no counter-rotating particles (Fig. 5). Figure 8. Open in new tabDownload slide Left: Histogram of |$\widehat{\boldsymbol {r}} \cdot \widehat{\boldsymbol {v}}$|⁠, where |$\boldsymbol {r}$| and |$\boldsymbol {v}$| are, respectively, the position and velocity of each particle relative to the MW. This is a basic measure of orbital eccentricity. We only show particles corresponding to the MW satellite plane according to the criteria in Table 6. The dashed blue line shows the expected distribution for particles on random orbits around the MW. Right: Analogous results for M31. All particles shown are in the corotating peak around M31 as there are no counter-rotating particles (Fig. 5). 5.2.5 Present radial distribution The right-hand sides of the panels in Fig. 6 (the right of vertical green lines) show the present radial distribution of MW satellite particles in our best-fitting model. This is similar to the range of Galactocentric distances at which we observe MW satellites in its satellite plane (Pawlowski & Kroupa 2013, fig. 2). The most luminous MW satellite is the LMC, whose Galactocentric distance is 50 kpc (Pietrzyński et al. 2013). This is certainly a distance at which our simulations indicate there should be a significant amount of tidal debris from the MW–M31 flyby. Our model indicates that the counter-rotating particles do not go as far out as the corotating particles. This is similar to the result obtained by Pawlowski et al. (2011), as can be seen from their fig. 6. The observed counter-rotator Sculptor (Sohn et al. 2017) has a Galactocentric distance of 86 kpc (Martínez-Vázquez et al. 2015), well within the simulated counter-rotating material’s final radial distribution. Thus, both the corotating and counter-rotating MW satellite particles in our simulation have a radial distribution broadly similar to what we observe for actual MW satellites orbiting within its satellite plane. Fig. 7 shows the initial and final radial distribution of M31 satellite particles. We do not show results for counter-rotating M31 satellite particles as our simulation lacks such particles (the bottom panel of Fig. 5). At the present time, our best-fitting model indicates that M31 satellite particles should have a broad range of distances from M31. This is similar to the range of observed distances between M31 and its satellites within its satellite plane (Ibata et al. 2013; fig. 1).8 Indeed, our simulated M31 satellite particle radial distribution peaks at ≈200 kpc, very similar to the observed 192 kpc rms radial extent of the M31 satellite plane (Ibata et al. 2014, section 3). Their results suggest that this structure does not extend much further out, so we impose a maximum distance limit of 250 kpc from M31 when selecting particles belonging to its simulated satellite plane. 5.2.6 Orbital eccentricity To better understand the orbits of satellite plane particles in our simulations, we consider how eccentric their orbits are. This is difficult to quantify rigorously due to the non-Keplerian potential (equation 52). As a proxy, we show the cosine distribution of angles between the position and velocity of each particle with respect to its host galaxy (Fig. 8), with weights assigned as in Section 3. For the MW, our results are similar to what they would be if its satellite particles had randomly distributed velocities (dashed blue line). The only major difference is a deficit of particles on extremely radial orbits i.e. whose velocity is within cos −10.9 = 26° of the radial direction. This is unlike the actual MW satellite plane as that is corotating and mostly rotation-supported (Pawlowski & Kroupa 2013). The discrepancy may be due to the lack of dissipative processes in our simulation. Even so, the nearly symmetric distribution of |$\widehat{\boldsymbol {r}} \cdot \widehat{\boldsymbol {v}}$| suggests that the MW satellite system has reached some sort of equilibrium in the ≈8 Gyr since its interaction with M31. We obtain rather different results for the satellite particles around M31 (the right-hand panel of Fig. 8). Despite the lack of dissipation, we get much fewer particles on radial orbits than would be expected for randomly distributed velocities. The clear peak near 0 (tangential motion) implies near-circular orbits. This provides a natural explanation for the clear radial velocity gradient across the M31 satellite plane, whereby satellites on different sides of M31 almost always have a different sign of radial velocity relative to that of M31 (Ibata et al. 2013). Our models suggest that the satellite plane members of M31 ought to be on more nearly circular orbits than those of the MW. Observationally, even the MW satellite plane is rotationally supported (Pawlowski & Kroupa 2013). Thus, we would certainly expect the M31 satellite plane to be predominantly rotation-supported. It will be interesting to test this prediction once proper motions of M31 satellites become available. 5.3 The MW thick disc An important constraint on our models is the extent to which they disrupt the well-observed MW disc in which we live (Section 4.3). For comparison with observations, we extract the rms thickness of the plane best-fitting MW particles at Galactocentric cylindrical radii rplane within 1 kpc of the Solar value (Ferguson et al. 2017, fig. 12). In our best model, this agrees fairly well with observations (Table 8). To check this result, we show the vertical distribution of these particles (Fig. 9). Their distance from the MW disc follows a nearly Gaussian distribution whose width is very similar to the rms thickness inferred by our plane-fitting algorithm (Banik & Zhao 2018b, section 5.1). Figure 9. Open in new tabDownload slide Weighted histogram of the vertical distances of MW particles from its disc. Here, we only show particles with cylindrical radius between 7.2 and 9.2 kpc, the region constrained by the observations of Ferguson et al. (2017). The solid blue line shows a Gaussian whose width is the thickness of the plane best-fitting these particles (Table 8). Figure 9. Open in new tabDownload slide Weighted histogram of the vertical distances of MW particles from its disc. Here, we only show particles with cylindrical radius between 7.2 and 9.2 kpc, the region constrained by the observations of Ferguson et al. (2017). The solid blue line shows a Gaussian whose width is the thickness of the plane best-fitting these particles (Table 8). The lower scale height than observed may be due to the Galactic disc having been thickened by processes other than tides raised during the M31 flyby (Section 4.3.2). Although we allow for this to some extent (equation 47), the MW disc may have experienced a burst of star formation due to the close passage of M31 (Snaith et al. 2014). This may have heated the Galactic disc via ‘popping’ star clusters (Kroupa 2002; Assmann et al. 2011). As our model does not include this process, it should find a dynamically colder MW than observed, even if secular heating is accounted for. In the future, it will become possible to obtain detailed observations of a larger region of the thick disc. We thus use Fig. 10 to show the rms height of our simulated MW disc for several bins in rplane, using a similar method to that used in Section 4.3.2. This reveals that the thick disc ought to be even thicker further from the MW, a natural consequence of its tidal origin. Thus, the inner regions of the MW should suffer less disruption than its more weakly bound outer parts. Figure 10. Open in new tabDownload slide We use red diamonds to show the rms thickness of MW disc material in various annuli (bins in Galactocentric radii), found using the plane-fitting method described in section 5.1 of Banik & Zhao (2018b) augmented by σ-clipping. The blue dots show the mean height (displacement towards the North Galactic Pole) of particles in each annulus relative to the position of the MW particle, our proxy for the position of the Galactic Centre. The vertical dashed line at 17.7 kpc is the radius which encloses 95 per cent of the MW mass in our best model. We have likely underestimated the thickness of the innermost annulus shown here as we terminate the trajectories of all particles that get within 2.15 kpc of the MW. Figure 10. Open in new tabDownload slide We use red diamonds to show the rms thickness of MW disc material in various annuli (bins in Galactocentric radii), found using the plane-fitting method described in section 5.1 of Banik & Zhao (2018b) augmented by σ-clipping. The blue dots show the mean height (displacement towards the North Galactic Pole) of particles in each annulus relative to the position of the MW particle, our proxy for the position of the Galactic Centre. The vertical dashed line at 17.7 kpc is the radius which encloses 95 per cent of the MW mass in our best model. We have likely underestimated the thickness of the innermost annulus shown here as we terminate the trajectories of all particles that get within 2.15 kpc of the MW. However, the trend may well be opposite this if the MW disc formed through internal processes e.g. interactions between stars, whose spatial density is higher closer to the Galactic Centre. Thus, it may be possible to determine if the thick disc originated in a tidal heating event. If it did not, this would cast doubt on our scenario. Interestingly, the thick disc does appear to have a higher velocity dispersion further from the MW (Banik 2014). The apparent lack of accreted stars suggests that a flyby may be a better explanation than a minor merger of the sort typical in ΛCDM (Ruchti et al. 2015). The expected 7.65 Gyr age of the thick disc in our model is lower than the 9 ± 1 Gyr estimated by (Quillen & Garnett 2001). Although this discrepancy is not significant, more recent observations suggest that the thick disc has an age closer to 7 Gyr (Yu & Liu 2018), more in line with our model. 5.4 Analogues to the LMC If our model yields a realistic satellite plane around the MW, it ought to yield particles on orbits similar to its most massive member, the LMC. To look for analogues of the LMC, we search for particles whose orbit around the MW takes them to the LMC’s Galactocentric distance rLMC with at least as much radial kinetic energy and specific angular momentum h as the LMC (measuring distances r and velocities v with respect to the MW). This allows us to consider particles on orbits similar to that of the LMC even if they are currently at a different point on their orbit, e.g. near apocentre rather than pericentre. We thus make use of the adjusted energy \begin{eqnarray*} \tilde{E} \ \equiv \ \overbrace{\Phi \left( r \right) + \frac{1}{2}v^2}^E \ -\ \overbrace{\frac{h^2}{2{r_{\rm{LMC}}}^2}}^{\begin{array}{cc}\text{Tangential kinetic} \\\text{energy at $r = r_{\rm{LMC}}$} \end{array}}. \end{eqnarray*} (51) For particles close to the MW, their motion is mainly governed by their total energy E and angular momentum h with respect to the MW as M31 is very distant. Up to an additive constant, |$\tilde{E}$| is thus the specific radial kinetic energy of a particle when it has a Galactocentric distance of rLMC (if |$\tilde{E}$| is too low, then the particle is not capable of reaching this distance). We define an analogue of the LMC as a particle whose |$\tilde{E}$| and h both exceed the LMC value. The constraint on |$\tilde{E}$| ensures that analogues can reach r = rLMC with its observed radial velocity while the constraint on h ensures analogues would do so at the correct tangential velocity. A LMC analogue selected in this way can have a very different present position and velocity to the LMC, as long as the particle belongs to the simulated MW satellite plane. If we neglect M31 and the EF, the form of the Galactic potential Φ(r) follows from integrating the simple interpolating function (Famaey & Binney 2005) for a single point mass using a hyperbolic substitution. \begin{eqnarray*} \Phi \left( r \right) &=& \sqrt{GMa_{0}} \left[ {\rm{ln}} \left(1 + \sqrt{1 + \tilde{r}^2} \right) - \frac{1}{\tilde{r}} - \sqrt{\frac{1}{\tilde{r}^2} + 1} \right], \nonumber \\\tilde{r} &=& \frac{2r}{r_{M}} \ \ \text{, where } r_{M} = \frac{\sqrt{GM}}{a_{0}}. \end{eqnarray*} (52) Here, M is the mass of the MW while rM is its ‘MOND radius’, beyond which the effects of MOND become significant. For the MW, rM ≈ 9 kpc, not much smaller than the LMC distance of rLMC = 50 kpc (Pietrzyński et al. 2013). This makes it inaccurate to assume that the LMC is in the deep-MOND limit. In Fig. 11, we show the values of |$\tilde{E}$| and h for particles in the satellite region of the MW. It is evident that a large number of particles have orbits analogous to that of the LMC using our definition (above and to the right of the black dot representing the LMC). In fact, the LMC is nearly in the middle of the sea of red dots representing simulated particles. Thus, our best-fitting model suggests that the LMC may well have condensed out of material in the outer part of the MW disc that was raised on to a higher and more inclined orbit by tides from M31. If our simulation had included self-gravity, it is reasonable to suppose that these particles on LMC-like orbits would have formed into a LMC-like object. This is similar to the scenario proposed by Yang et al. (2014) in which the LMC is a tidal dwarf galaxy. Figure 11. Open in new tabDownload slide Values of |$\tilde{E}$| (equation 51) against h (equation 36) relative to the LMC value for MW satellite particles (red dots, only part of graph shown). The blue curve shows the relation for particles on circular orbits, with text labels along this curve (next to blue dots) giving the corresponding Galactocentric distance in kpc. We consider particles analogous to the LMC if their |$\tilde{E}$| and h both exceed the LMC value (black dot, uncertainty indicated with dashed line). We do not show the uncertainty in |$\tilde{E}$| for the LMC because we neglect the rather small uncertainty in its radial velocity, whereas its tangential velocity is mainly a function of its less accurately known proper motion (Kallivayalil et al. 2013, table 5). Our adopted distance to the LMC is 49.39 kpc based on a heliocentric distance of 49.97 kpc (Pietrzyński et al. 2013) and a Sun-Galactic Centre distance of 8.20 kpc (McMillan 2017). Figure 11. Open in new tabDownload slide Values of |$\tilde{E}$| (equation 51) against h (equation 36) relative to the LMC value for MW satellite particles (red dots, only part of graph shown). The blue curve shows the relation for particles on circular orbits, with text labels along this curve (next to blue dots) giving the corresponding Galactocentric distance in kpc. We consider particles analogous to the LMC if their |$\tilde{E}$| and h both exceed the LMC value (black dot, uncertainty indicated with dashed line). We do not show the uncertainty in |$\tilde{E}$| for the LMC because we neglect the rather small uncertainty in its radial velocity, whereas its tangential velocity is mainly a function of its less accurately known proper motion (Kallivayalil et al. 2013, table 5). Our adopted distance to the LMC is 49.39 kpc based on a heliocentric distance of 49.97 kpc (Pietrzyński et al. 2013) and a Sun-Galactic Centre distance of 8.20 kpc (McMillan 2017). In this case, the internal dynamics of the LMC need to be explained using its baryons alone as tidal dwarf galaxies are expected to be purely baryonic (Wetzstein et al. 2007). This is very difficult to do in Newtonian gravity (van der Marel et al. 2002), a result that was confirmed recently with proper motion-based measurements of the LMC rotation curve (van der Marel & Kallivayalil 2014; van der Marel & Sahlmann 2016). On the other hand, fig. 8 of van der Marel & Kallivayalil (2014) shows that the LMC falls on the baryonic Tully–Fisher relation (equation 50), implying consistency with MOND regardless of how the LMC formed. In our model, the LMC is bound to the MW. This is also evident if we look at fig. 2 of Patel et al. (2017) and put in a MW mass of 1.5 × 1012 M⊙ corresponding to what Newtonian gravity requires in order to sustain a flat rotation curve of 180 km s−1 out to 200 kpc. The bound nature of the LMC is unsurprising given that the observed MW escape velocity at the LMC position (380 ± 30 km s−1, Williams et al. 2017a) exceeds its Galactocentric velocity (321 ± 24 km s−1, Kallivayalil et al. 2013, table 5). In MOND, a bound LMC must have been orbiting the MW for ≈8 Gyr as there is almost no dynamical friction 40 kpc from a purely baryonic MW. Assuming the Small Magellanic Cloud (SMC) formed on a similar orbit as part of the MW satellite plane, it would be quite feasible for these two objects to closely interact, helping to explain the observed properties of the Magellanic Stream (Hammer et al. 2015). A common origin during a past MW–M31 flyby provides a natural explanation for the alignment between the Galactic orbit of the LMC, the LMC-SMC orbit and the MW satellite plane (Pawlowski, McGaugh & Jerjen 2015a, fig. 1). It has long been assumed that the LMC and SMC have been bound to the MW for several Gyr (Diaz & Bekki 2011), consistent with our expectation. If instead they were recently accreted by the MW due to dynamical friction with its dark matter halo (Besla et al. 2007), then these geometrical alignments would be fortuitous. After all, the MW satellite plane orientation would be a priori unrelated to the LMC-SMC binary orbit which could itself be oriented rather differently to the motion of their barycentre about the MW (Kroupa 2015, section 5.3.1). The competing scenarios might be distinguished using more detailed observations of the Magellanic Stream (Fox 2016). 6 CONCLUSIONS Using a MOND model of the LG, we showed that the MW–M31 trajectory could be made consistent with cosmological initial conditions, a technique known as the timing argument. Such trajectories always involve a past close MW–M31 flyby (Fig. 4). Assuming a standard ΛCDM cosmology (equation 21), the MW–M31 trajectory satisfies the timing argument with a total implied mass close to the sum of their observed baryonic masses and a close pericentre around the time when the Galactic thick disc formed (Table 9). This conclusion would be affected somewhat if the actual cosmology is non-standard, as seems likely if the MW and M31 indeed had a past close flyby. However, any viable model should yield an expansion history similar to the standard one given that this works well observationally (e.g. Joudaki et al. 2017). Using a restricted N-body model of the MW and M31 discs advanced according to our MW–M31 trajectory, we simulated how their close flyby would have affected their internal structure. By performing a grid of models covering the parameter space, we obtained a reasonably good match to the observed properties of the MW and M31 (Table 8). In particular, we found clear peaks in their tidal debris orbital pole distributions at orientations similar to those observed for their satellite planes (Fig. 5). Our model correctly predicts a much smaller angular dispersion for the M31 satellite plane compared to the MW (Pawlowski et al. 2014, tables 2 and 4). It also matches the observation that the M31 satellite plane is ≈3 orders of magnitude less massive than that of the MW. However, it underestimates the masses of both satellite planes by about an order of magnitude. This may be related to the rather strong sensitivity of our models to the disc surface density of the MW and M31 at radii of  ∼  50 kpc (Section 5.2.4). Our model naturally yields counter-rotating material in the satellite plane around the MW, though not around M31. Thus, we can naturally explain why Sculptor counter-rotates within the MW satellite plane (Pawlowski et al. 2011; Sohn et al. 2017). However, our simulated MW satellite plane has a rather high fraction of counter-rotators. Interestingly, the observed MW satellite plane has a smaller radial extent than in our simulations, as can be seen by comparing our Fig. 6 with fig. 2 of Pawlowski & Kroupa (2013). Both these issues may be related to the lack of dissipative processes in our simulations. We hope to conduct more detailed investigations like those of Bílek et al. (2017) using the Phantom of ramses algorithm (Lüghausen, Famaey & Kroupa 2015) which adapts the gravity solver of the ramsescode (Teyssier 2002) widely used by astronomers to conduct N-body simulations that include hydrodynamic effects. Assuming these yield similar results to our flyby model, it may be tested based on proper motions of M31 satellites in its satellite plane. In particular, two of these (And XIII & And XXVII) do not share the coherent radial velocity trend of the remaining 13 M31 satellites in this structure (Ibata et al. 2013). If they turn out to orbit within the M31 satellite plane, then this implies a significant value of fcounter for M31, thus challenging our model (Table 8). We expect instead that these two dwarfs are just interlopers whose orbits take them far from the M31 satellite plane, as suggested by Ibata et al. (2013) in their supplementary information. This seems quite feasible given that half the M31 satellites do not belong to its satellite plane and have a more isotropic distribution. It has proven difficult to reconcile the ΛCDM paradigm with the highly flattened satellite galaxy planes of the MW and Andromeda while also explaining the high internal velocity dispersions of their constituent satellites. Using the alternative paradigm of MOND, we obtained a reasonable match to several important properties of these structures, in particular their 3D orientations. However, our restricted N-body model misses dissipative and self-gravitating effects which are likely important for accurate models of galactic interactions. Thus, more work will be required to test this scenario. Footnotes 1 The particle mass would have to be 10 − 80 MeV/c2 and it would need to carry a small but non-zero fraction of the electric charge. The cross-section for scattering with gas also needs to decline steeply with velocity. 2 |$v_{\rm{pec}} \propto \frac{1}{a \left( t \right)}$| for a test particle in an otherwise homogeneous universe, making avpec constant. 3 Ωm, 0 + ΩΛ, 0 = 1. 4 Mergers are possible in MOND as dynamical friction still exists when the baryonic parts of two galaxies overlap. As these are much smaller than galactic dark matter haloes in ΛCDM, mergers would be much rarer in MOND, though its stronger long-range gravity would lead to more gravitational focusing. 5 We accelerate this process by assuming a change to the initial disc spin of either galaxy causes an equal change to its final disc spin (e.g. tilting the MW 5° further north initially causes its present spin to also tilt 5° further north). This greatly increases the reliability of our algorithm as well as removing the burden of estimating derivatives through finite differencing, reducing the computational cost 3×. 6 If this threshold angle exceeds 90°, we do not perform outlier rejection – this often happens for the MW. 7 We use ∼ for order of magnitude estimates and ≈ when making more precise statements. 8 At the 783 kpc distance of M31 (McConnachie 2012), each square in this figure has a side length of 55 kpc. ACKNOWLEDGEMENTS The authors are grateful to the referee for her/his constructive comments. IB was supported by Science and Technology Facilities Council studentship 1506672. The work of DOR was enabled by Royal Astronomical Society undergraduate research bursary 6390/302/001. The algorithms were set up using matlab®. REFERENCES Ahmed S. H. , Brooks A. M. , Christensen C. R. , 2017 , MNRAS , 466 , 3119 https://doi.org/10.1093/mnras/stw3271 Crossref Search ADS Alves D. R. , Nelson C. A. , 2000 , ApJ , 542 , 789 https://doi.org/10.1086/317023 Crossref Search ADS Aragon-Calvo M. A. , Silk J. , Szalay A. S. , 2011 , MNRAS , 415 , L16 https://doi.org/10.1111/j.1745-3933.2011.01071.x Crossref Search ADS Assmann P. , Fellhauer M. , Kroupa P. , Brüns R. C. , Smith R. , 2011 , MNRAS , 415 , 1280 https://doi.org/10.1111/j.1365-2966.2011.18773.x Crossref Search ADS Banik I. , 2014 , preprint (arXiv:1406.4538) Banik I. , Zhao H. , 2015 , preprint (arXiv:1509.08457) Banik I. , Zhao H. , 2016 , MNRAS , 459 , 2237 https://doi.org/10.1093/mnras/stw787 Crossref Search ADS Banik I. , Zhao H. , 2017 , MNRAS , 467 , 2180 Banik I. , Zhao H. , 2018a , MNRAS , 473 , 419 https://doi.org/10.1093/mnras/stx2350 Crossref Search ADS Banik I. , Zhao H. , 2018b , MNRAS , 473 , 4033 https://doi.org/10.1093/mnras/stx2596 Crossref Search ADS Barnes J. E. , Hernquist L. , 1992 , Nature , 360 , 715 https://doi.org/10.1038/360715a0 Crossref Search ADS Bekenstein J. , Milgrom M. , 1984 , ApJ , 286 , 7 https://doi.org/10.1086/162570 Crossref Search ADS Bellazzini M. , Oosterloo T. , Fraternali F. , Beccari G. , 2013 , A&A , 559 , L11 Crossref Search ADS Belokurov V. et al. , 2014 , MNRAS , 437 , 116 https://doi.org/10.1093/mnras/stt1862 Crossref Search ADS Bengaly C. A. P. , Novaes C. P. , Xavier H. S. , Bilicki M. , Bernui A. , Alcaniz J. S. , 2018 , MNRAS , 475 , L106 https://doi.org/10.1093/mnrasl/sly002 Crossref Search ADS Berlin A. , Hooper D. , Krnjaic G. , McDermott S. D. , 2018 , preprint (arXiv:1803.02804) Besla G. , Kallivayalil N. , Hernquist L. , Robertson B. , Cox T. J. , van der Marel R. P. , Alcock C. , 2007 , ApJ , 668 , 949 https://doi.org/10.1086/521385 Crossref Search ADS Bílek M. , Thies I. , Kroupa P. , Famaey B. , 2017 , A&A , forthcoming Bovy J. , Rix H.-W. , 2013 , ApJ , 779 , 115 https://doi.org/10.1088/0004-637X/779/2/115 Crossref Search ADS Bowman J. D. , Rogers A. E. E. , Monsalve R. A. , Mozdzen T. J. , Mahesh N. , 2018 , Nature , 555 , 67 https://doi.org/10.1038/nature25792 Crossref Search ADS PubMed Brunthaler A. , Reid M. J. , Falcke H. , Greenhill L. J. , Henkel C. , 2005 , Science , 307 , 1440 https://doi.org/10.1126/science.1108342 Crossref Search ADS PubMed Butsky I. et al. , 2016 , MNRAS , 462 , 663 https://doi.org/10.1093/mnras/stw1688 Crossref Search ADS Candlish G. N. , 2016 , MNRAS , 460 , 2571 https://doi.org/10.1093/mnras/stw1130 Crossref Search ADS Carignan C. , Chemin L. , Huchtmeier W. K. , Lockman F. J. , 2006 , ApJ , 641 , L109 https://doi.org/10.1086/503869 Crossref Search ADS Carlberg R. G. , Sellwood J. A. , 1985 , ApJ , 292 , 79 https://doi.org/10.1086/163134 Crossref Search ADS Chae K.-H. , Bernardi M. , Sheth R. K. , 2017 , preprint (arXiv:1707.08280) Chemin L. , Carignan C. , Foster T. , 2009 , ApJ , 705 , 1395 https://doi.org/10.1088/0004-637X/705/2/1395 Crossref Search ADS Chiboucas K. , Jacobs B. A. , Tully R. B. , Karachentsev I. D. , 2013 , AJ , 146 , 126 https://doi.org/10.1088/0004-6256/146/5/126 Crossref Search ADS Corbelli E. , 2003 , MNRAS , 342 , 199 https://doi.org/10.1046/j.1365-8711.2003.06531.x Crossref Search ADS Corbelli E. , Salucci P. , 2007 , MNRAS , 374 , 1051 https://doi.org/10.1111/j.1365-2966.2006.11219.x Crossref Search ADS Courteau S. , Widrow L. M. , McDonald M. , Guhathakurta P. , Gilbert K. M. , Zhu Y. , Beaton R. L. , Majewski S. R. , 2011 , ApJ , 739 , 20 https://doi.org/10.1088/0004-637X/739/1/20 Crossref Search ADS Dabringhausen J. , Kroupa P. , Famaey B. , Fellhauer M. , 2016 , MNRAS , 463 , 1865 https://doi.org/10.1093/mnras/stw2001 Crossref Search ADS Desmond H. , 2017a , MNRAS , 464 , 4160 https://doi.org/10.1093/mnras/stw2571 Crossref Search ADS Desmond H. , 2017b , MNRAS , 472 , L35 https://doi.org/10.1093/mnrasl/slx134 Crossref Search ADS Diaz J. , Bekki K. , 2011 , MNRAS , 413 , 2015 https://doi.org/10.1111/j.1365-2966.2011.18289.x Crossref Search ADS Dierickx M. , Blecha L. , Loeb A. , 2014 , ApJ , 788 , L38 https://doi.org/10.1088/2041-8205/788/2/L38 Crossref Search ADS Einasto J. , Lynden-Bell D. , 1982 , MNRAS , 199 , 67 https://doi.org/10.1093/mnras/199.1.67 Crossref Search ADS Ewall-Wice A. , Chang T.-C. , Lazio J. , Dore O. , Seiffert M. , Monsalve R. A. , 2018 , ApJ , preprint (arXiv:1803.01815) Famaey B. , Binney J. , 2005 , MNRAS , 363 , 603 https://doi.org/10.1111/j.1365-2966.2005.09474.x Crossref Search ADS Famaey B. , McGaugh S. S. , 2012 , Living Rev. Relativ. , 15 , 10 https://doi.org/10.12942/lrr-2012-10 Crossref Search ADS PubMed Famaey B. , Bruneton J.-P. , Zhao H. , 2007 , MNRAS , 377 , L79 https://doi.org/10.1111/j.1745-3933.2007.00308.x Crossref Search ADS Ferguson D. , Gardner S. , Yanny B. , 2017 , ApJ , 843 , 141 https://doi.org/10.3847/1538-4357/aa77fd Crossref Search ADS Fernando N. , Arias V. , Guglielmo M. , Lewis G. F. , Ibata R. A. , Power C. , 2017 , MNRAS , 465 , 641 https://doi.org/10.1093/mnras/stw2694 Crossref Search ADS Fox A. , 2016 , The Origin of the Leading Arm of the Magellanic Stream, HST Proposal . Google Preview WorldCat COPAC Francis C. , Anderson E. , 2014 , Celest. Mech. Dyn. Astron. , 118 , 399 https://doi.org/10.1007/s10569-014-9541-z Crossref Search ADS Garaldi E. , Romano-Díaz E. , Borzyszkowski M. , Porciani C. , 2018 , MNRAS , 473 , 2234 https://doi.org/10.1093/mnras/stx2489 Crossref Search ADS Gerke J. R. , Kochanek C. S. , Prieto J. L. , Stanek K. Z. , Macri L. M. , 2011 , ApJ , 743 , 176 https://doi.org/10.1088/0004-637X/743/2/176 Crossref Search ADS Gilmore G. , Reid N. , 1983 , MNRAS , 202 , 1025 https://doi.org/10.1093/mnras/202.4.1025 Crossref Search ADS González R. E. , Padilla N. D. , 2010 , MNRAS , 407 , 1449 https://doi.org/10.1111/j.1365-2966.2010.17015.x Crossref Search ADS Hammer F. , Yang Y. , Fouquet S. , Pawlowski M. S. , Kroupa P. , Puech M. , Flores H. , Wang J. , 2013 , MNRAS , 431 , 3543 https://doi.org/10.1093/mnras/stt435 Crossref Search ADS Hammer F. , Yang Y. B. , Flores H. , Puech M. , Fouquet S. , 2015 , ApJ , 813 , 110 https://doi.org/10.1088/0004-637X/813/2/110 Crossref Search ADS Harris G. L. H. , Rejkuba M. , Harris W. E. , 2010 , PASA , 27 , 457 https://doi.org/10.1071/AS09061 Crossref Search ADS Hayden M. R. et al. , 2015 , ApJ , 808 , 132 https://doi.org/10.1088/0004-637X/808/2/132 Crossref Search ADS Heinrich C. , Hu W. , 2018 , preprint (arXiv:1802.00791) Ibata R. A. et al. , 2013 , Nature , 493 , 62 https://doi.org/10.1038/nature11717 Crossref Search ADS PubMed Ibata R. A. , Ibata N. G. , Lewis G. F. , Martin N. F. , Conn A. , Elahi P. , Arias V. , Fernando N. , 2014 , ApJ , 784 , L6 https://doi.org/10.1088/2041-8205/784/1/L6 Crossref Search ADS Iocco F. , Pato M. , Bertone G. , 2015 , Phys. Rev. D , 92 , 084046 https://doi.org/10.1103/PhysRevD.92.084046 Crossref Search ADS Javanmardi B. , Kroupa P. , 2017 , A&A , 597 , A120 Crossref Search ADS Javanmardi B. , Porciani C. , Kroupa P. , Pflamm-Altenburg J. , 2015 , ApJ , 810 , 47 https://doi.org/10.1088/0004-637X/810/1/47 Crossref Search ADS Jayaraman A. , Gilmore G. , Wyse R. F. G. , Norris J. E. , Belokurov V. , 2013 , MNRAS , 431 , 930 https://doi.org/10.1093/mnras/stt221 Crossref Search ADS Joudaki S. , Kaplinghat M. , Keeley R. , Kirkby D. , 2017 , preprint (arXiv:1710.04236) Jurić M. et al. , 2008 , ApJ , 673 , 864 https://doi.org/10.1086/523619 Crossref Search ADS Kafle P. R. , Sharma S. , Lewis G. F. , Bland-Hawthorn J. , 2012 , ApJ , 761 , 98 https://doi.org/10.1088/0004-637X/761/2/98 Crossref Search ADS Kahn F. D. , Woltjer L. , 1959 , ApJ , 130 , 705 https://doi.org/10.1086/146762 Crossref Search ADS Kallivayalil N. , van der Marel R. P. , Besla G. , Anderson J. , Alcock C. , 2013 , ApJ , 764 , 161 https://doi.org/10.1088/0004-637X/764/2/161 Crossref Search ADS Karachentsev I. D. , 2005 , AJ , 129 , 178 https://doi.org/10.1086/426368 Crossref Search ADS Karachentsev I. D. , Kashibadze O. G. , 2006 , Astrophysics , 49 , 3 https://doi.org/10.1007/s10511-006-0002-6 Crossref Search ADS Kogut A. et al. , 1993 , ApJ , 419 , 1 https://doi.org/10.1086/173453 Crossref Search ADS Kourkchi E. , Tully R. B. , 2017 , ApJ , 843 , 16 https://doi.org/10.3847/1538-4357/aa76db Crossref Search ADS Kroupa P. , 2002 , MNRAS , 330 , 707 https://doi.org/10.1046/j.1365-8711.2002.05128.x Crossref Search ADS Kroupa P. , 2012 , PASA , 29 , 395 https://doi.org/10.1071/AS12005 Crossref Search ADS Kroupa P. , 2015 , Can. J. Phys. , 93 , 169 https://doi.org/10.1139/cjp-2014-0179 Crossref Search ADS Kroupa P. , Theis C. , Boily C. M. , 2005 , A&A , 431 , 517 Crossref Search ADS Lehner N. , Howk J. C. , Wakker B. P. , 2015 , ApJ , 804 , 79 https://doi.org/10.1088/0004-637X/804/2/79 Crossref Search ADS Lelli F. , McGaugh S. S. , Schombert J. M. , Pawlowski M. S. , 2017 , ApJ , 836 , 152 https://doi.org/10.3847/1538-4357/836/2/152 Crossref Search ADS Lüghausen F. , Famaey B. , Kroupa P. , 2015 , Can. J. Phys. , 93 , 232 https://doi.org/10.1139/cjp-2014-0168 Crossref Search ADS Lynden-Bell D. , 1976 , MNRAS , 174 , 695 https://doi.org/10.1093/mnras/174.3.695 Crossref Search ADS Lynden-Bell D. , 1982 , Observatory , 102 , 202 Ma C. et al. , 1998 , AJ , 116 , 516 https://doi.org/10.1086/300408 Crossref Search ADS Maji M. , Zhu Q. , Marinacci F. , Li Y. , 2017 , preprint (arXiv:1702.00497) Martin N. F. et al. , 2013 , ApJ , 772 , 15 https://doi.org/10.1088/0004-637X/772/1/15 Crossref Search ADS Martin N. F. et al. , 2014 , ApJ , 793 , L14 https://doi.org/10.1088/2041-8205/793/1/L14 Crossref Search ADS Martin N. F. et al. , 2016 , ApJ , 833 , 167 https://doi.org/10.3847/1538-4357/833/2/167 Crossref Search ADS Martínez-Vázquez C. E. et al. , 2015 , MNRAS , 454 , 1509 https://doi.org/10.1093/mnras/stv2014 Crossref Search ADS McConnachie A. W. , 2012 , AJ , 144 , 4 https://doi.org/10.1088/0004-6256/144/1/4 Crossref Search ADS McConnachie A. W. et al. , 2009 , Nature , 461 , 66 https://doi.org/10.1038/nature08327 Crossref Search ADS PubMed McGaugh S. , 1999 , in Holt S. , Smith E. , eds, AIP Conf. Proc. Vol. 470, After the Dark Ages: When Galaxies Were Young (the Universe at 2 < z < 5). Am. Inst. Phys , New York , p. 72 Google Preview WorldCat COPAC McGaugh S. S. , 2011 , Phys. Rev. Lett. , 106 , 121303 https://doi.org/10.1103/PhysRevLett.106.121303 Crossref Search ADS PubMed McGaugh S. S. , 2018 , Res. Notes Am. Astron. Soc. , 2 , 37 Crossref Search ADS McGaugh S. , Milgrom M. , 2013 , ApJ , 775 , 139 https://doi.org/10.1088/0004-637X/775/2/139 Crossref Search ADS McGaugh S. S. , Wolf J. , 2010 , ApJ , 722 , 248 https://doi.org/10.1088/0004-637X/722/1/248 Crossref Search ADS McMillan P. J. , 2011 , MNRAS , 414 , 2446 https://doi.org/10.1111/j.1365-2966.2011.18564.x Crossref Search ADS McMillan P. J. , 2017 , MNRAS , 465 , 76 https://doi.org/10.1093/mnras/stw2759 Crossref Search ADS Metz M. , Kroupa P. , Jerjen H. , 2007 , MNRAS , 374 , 1125 https://doi.org/10.1111/j.1365-2966.2006.11228.x Crossref Search ADS Metz M. , Kroupa P. , Libeskind N. I. , 2008 , ApJ , 680 , 287 https://doi.org/10.1086/587833 Crossref Search ADS Metz M. , Kroupa P. , Jerjen H. , 2009 , MNRAS , 394 , 2223 https://doi.org/10.1111/j.1365-2966.2009.14489.x Crossref Search ADS Mieske S. , Hilker M. , Infante L. , 2005 , A&A , 438 , 103 Crossref Search ADS Milgrom M. , 1983 , ApJ , 270 , 365 https://doi.org/10.1086/161130 Crossref Search ADS Milgrom M. , 1986 , ApJ , 302 , 617 https://doi.org/10.1086/164021 Crossref Search ADS Milgrom M. , 1999 , Phys. Lett. A , 253 , 273 https://doi.org/10.1016/S0375-9601(99)00077-8 Crossref Search ADS Milgrom M. , 2010 , MNRAS , 403 , 886 https://doi.org/10.1111/j.1365-2966.2009.16184.x Crossref Search ADS Minchev I. et al. , 2014 , ApJ , 781 , L20 https://doi.org/10.1088/2041-8205/781/1/L20 Crossref Search ADS Mirabel I. F. , Dottori H. , Lutz D. , 1992 , A&A , 256 , L19 Müller O. , Scalera R. , Binggeli B. , Jerjen H. , 2017 , A&A , 602 , A119 Crossref Search ADS Müller O. , Pawlowski M. S. , Jerjen H. , Lelli F. , 2018 , Science , 359 , 534 https://doi.org/10.1126/science.aao1858 Crossref Search ADS PubMed Nakanishi H. , Sofue Y. , 2016 , PASJ , 68 , 5 https://doi.org/10.1093/pasj/psv108 Crossref Search ADS Nicastro F. , Senatore F. , Krongold Y. , Mathur S. , Elvis M. , 2016 , ApJ , 828 , L12 https://doi.org/10.3847/2041-8205/828/1/L12 Crossref Search ADS Ostriker J. P. , Steinhardt P. J. , 1995 , Nature , 377 , 600 https://doi.org/10.1038/377600a0 Crossref Search ADS Patel E. , Besla G. , Sohn S. T. , 2017 , MNRAS , 464 , 3825 https://doi.org/10.1093/mnras/stw2616 Crossref Search ADS Pawlowski M. S. , 2016 , MNRAS , 456 , 448 https://doi.org/10.1093/mnras/stv2673 Crossref Search ADS Pawlowski M. S. , 2018 , Mod. Phys. Lett. A , 33 , 1830004 https://doi.org/10.1142/S0217732318300045 Crossref Search ADS Pawlowski M. S. , Kroupa P. , 2013 , MNRAS , 435 , 2116 https://doi.org/10.1093/mnras/stt1429 Crossref Search ADS Pawlowski M. S. , Kroupa P. , 2014 , ApJ , 790 , 74 https://doi.org/10.1088/0004-637X/790/1/74 Crossref Search ADS Pawlowski M. S. , McGaugh S. S. , 2014 , MNRAS , 440 , 908 https://doi.org/10.1093/mnras/stu321 Crossref Search ADS Pawlowski M. S. , Kroupa P. , de Boer K. S. , 2011 , A&A , 532 , A118 Crossref Search ADS Pawlowski M. S. , Pflamm-Altenburg J. , Kroupa P. , 2012 , MNRAS , 423 , 1109 https://doi.org/10.1111/j.1365-2966.2012.20937.x Crossref Search ADS Pawlowski M. S. , Kroupa P. , Jerjen H. , 2013 , MNRAS , 435 , 1928 https://doi.org/10.1093/mnras/stt1384 Crossref Search ADS Pawlowski M. S. et al. , 2014 , MNRAS , 442 , 2362 https://doi.org/10.1093/mnras/stu1005 Crossref Search ADS Pawlowski M. S. , McGaugh S. S. , Jerjen H. , 2015a , MNRAS , 453 , 1047 https://doi.org/10.1093/mnras/stv1588 Crossref Search ADS Pawlowski M. S. , Famaey B. , Merritt D. , Kroupa P. , 2015b , ApJ , 815 , 19 https://doi.org/10.1088/0004-637X/815/1/19 Crossref Search ADS Pawlowski M. S. et al. , 2017 , Astron. Nachr. , 338 , 854 https://doi.org/10.1002/asna.201713366 Crossref Search ADS Pazy E. , 2013 , Phys. Rev. D , 87 , 084063 https://doi.org/10.1103/PhysRevD.87.084063 Crossref Search ADS Peebles P. J. E. , 2017 , preprint (arXiv:1705.10683) Peñarrubia J. , Ma Y.-Z. , Walker M. G. , McConnachie A. , 2014 , MNRAS , 443 , 2204 https://doi.org/10.1093/mnras/stu879 Crossref Search ADS Peters P. C. , 1981 , Am. J. Phys. , 49 , 564 https://doi.org/10.1119/1.12460 Crossref Search ADS Phelps S. , Nusser A. , Desjacques V. , 2013 , ApJ , 775 , 102 https://doi.org/10.1088/0004-637X/775/2/102 Crossref Search ADS Piatek S. , Pryor C. , Bristow P. , Olszewski E. W. , Harris H. C. , Mateo M. , Minniti D. , Tinney C. G. , 2006 , AJ , 131 , 1445 https://doi.org/10.1086/499526 Crossref Search ADS Piatek S. , Pryor C. , Olszewski E. W. , 2016 , AJ , 152 , 166 https://doi.org/10.3847/0004-6256/152/6/166 Crossref Search ADS Pietrzyński G. et al. , 2013 , Nature , 495 , 76 https://doi.org/10.1038/nature11878 Crossref Search ADS PubMed Planck Collaboration XIII , 2016 , A&A , 594 , A13 Crossref Search ADS Privon G. C. , Barnes J. E. , Evans A. S. , Hibbard J. E. , Yun M. S. , Mazzarella J. M. , Armus L. , Surace J. , 2013 , ApJ , 771 , 120 https://doi.org/10.1088/0004-637X/771/2/120 Crossref Search ADS Quillen A. C. , Garnett D. R. , 2001 , in Funes J. G. , Corsini E. M. , eds, ASP Conf. Ser. Vol. 230, Galaxy Disks and Disk Galaxies . Astron. Soc. Pac , San Francisco , p. 87 Google Preview WorldCat COPAC Renaud F. , Famaey B. , Kroupa P. , 2016 , MNRAS , 463 , 3637 https://doi.org/10.1093/mnras/stw2331 Crossref Search ADS Riess A. G. et al. , 1998 , AJ , 116 , 1009 https://doi.org/10.1086/300499 Crossref Search ADS Ruchti G. R. et al. , 2015 , MNRAS , 450 , 2874 https://doi.org/10.1093/mnras/stv807 Crossref Search ADS Sadoun R. , Mohayaee R. , Colin J. , 2014 , MNRAS , 442 , 160 https://doi.org/10.1093/mnras/stu850 Crossref Search ADS Salem M. , Besla G. , Bryan G. , Putman M. , van der Marel R. P. , Tonnesen S. , 2015 , ApJ , 815 , 77 https://doi.org/10.1088/0004-637X/815/1/77 Crossref Search ADS Salucci P. , Turini N. , 2017 , preprint (arXiv:1707.01059) Sandage A. , 1986 , ApJ , 307 , 1 https://doi.org/10.1086/164387 Crossref Search ADS Shao S. , Cautun M. , Frenk C. S. , Grand R. J. J. , Gómez F. A. , Marinacci F. , Simpson C. M. , 2018 , MNRAS , 476 , 1796 https://doi.org/10.1093/mnras/sty343 Crossref Search ADS Smolin L. , 2017 , Phys. Rev. D , 96 , 083523 Snaith O. N. , Haywood M. , Di Matteo P. , Lehnert M. D. , Combes F. , Katz D. , Gómez A. , 2014 , ApJ , 781 , L31 https://doi.org/10.1088/2041-8205/781/2/L31 Crossref Search ADS Sohn S. T. et al. , 2017 , ApJ , 849 , 93 https://doi.org/10.3847/1538-4357/aa917b Crossref Search ADS Teyssier R. , 2002 , A&A , 385 , 337 Crossref Search ADS Thomas G. F. , Famaey B. , Ibata R. , Lüghausen F. , Kroupa P. , 2017 , A&A , 603 , A65 Crossref Search ADS Tiret O. , Combes F. , 2008 , in Funes J. G. , Corsini E. M. , eds, ASP Conf. Ser. Vol. 396, Formation and Evolution of Galaxy Disks . Astron. Soc. Pac , San Francisco , p. 259 Google Preview WorldCat COPAC van der Marel R. P. , Kallivayalil N. , 2014 , ApJ , 781 , 121 https://doi.org/10.1088/0004-637X/781/2/121 Crossref Search ADS van der Marel R. P. , Sahlmann J. , 2016 , ApJ , 832 , L23 https://doi.org/10.3847/2041-8205/832/2/L23 Crossref Search ADS van der Marel R. P. , Alves D. R. , Hardy E. , Suntzeff N. B. , 2002 , AJ , 124 , 2639 https://doi.org/10.1086/343775 Crossref Search ADS van der Marel R. P. , Besla G. , Cox T. J. , Sohn S. T. , Anderson J. , 2012 , ApJ , 753 , 9 https://doi.org/10.1088/0004-637X/753/1/9 Crossref Search ADS Verlinde E. P. , 2017 , SciPost Physics , 2 , 016 Weisz D. R. , Dolphin A. E. , Skillman E. D. , Holtzman J. , Gilbert K. M. , Dalcanton J. J. , Williams B. F. , 2014 , ApJ , 789 , 147 https://doi.org/10.1088/0004-637X/789/2/147 Crossref Search ADS Wetzstein M. , Naab T. , Burkert A. , 2007 , MNRAS , 375 , 805 https://doi.org/10.1111/j.1365-2966.2006.11360.x Crossref Search ADS Williams A. A. , Belokurov V. , Casey A. R. , Evans N. W. , 2017a , MNRAS , 468 , 2359 https://doi.org/10.1093/mnras/stx508 Crossref Search ADS Williams B. F. et al. , 2017b , ApJ , 846 , 145 https://doi.org/10.3847/1538-4357/aa862a Crossref Search ADS Wu P.-F. , Tully R. B. , Rizzi L. , Dolphin A. E. , Jacobs B. A. , Karachentsev I. D. , 2014 , AJ , 148 , 7 https://doi.org/10.1088/0004-6256/148/1/7 Crossref Search ADS Wyse R. F. G. , 2009 , in Mamajek E. E. , Soderblom D. R. , Wyse R. F. G. , eds, Proc. IAU Symp. 258, The Ages of Stars , Cambridge Univ. Press , Cambridge . p. 11 Google Preview WorldCat COPAC Yang Y. , Hammer F. , Fouquet S. , Flores H. , Puech M. , Pawlowski M. S. , Kroupa P. , 2014 , MNRAS , 442 , 2419 https://doi.org/10.1093/mnras/stu931 Crossref Search ADS Yu J. , Liu C. , 2018 , MNRAS , 475 , 1093 https://doi.org/10.1093/mnras/stx3204 Crossref Search ADS Zhao H. , 2007 , ApJ , 671 , L1 https://doi.org/10.1086/524731 Crossref Search ADS Zhao H. , Li B. , Bienaymé O. , 2010 , Phys. Rev. D , 82 , 103001 https://doi.org/10.1103/PhysRevD.82.103001 Crossref Search ADS Zhao H. , Famaey B. , Lüghausen F. , Kroupa P. , 2013 , A&A , 557 , L3 Crossref Search ADS APPENDIX A THE TANGENTIAL FORCE BETWEEN TWO MASSES In Section 2.3.2, we used an approximation for the relative tangential gravity |$\boldsymbol {g}_{\rm{tan}}$| between two masses. This is the difference in the gravitational field experienced by the masses, less the component of this parallel to the line connecting them. In Newtonian gravity, |$\boldsymbol {g}_{\rm{tan}} = 0$|⁠. This is also true in MOND if the bodies are isolated. However, |$\boldsymbol {g}_{\rm{tan}} \ne 0$| in the presence of an external field, even if one of the objects can be considered a test particle (Banik & Zhao 2015, equation 39). Our estimate for |$\boldsymbol {g}_{\rm{tan}}$| is \begin{eqnarray*} \boldsymbol {g}_{\rm{tan}} = \left( \widehat{\boldsymbol {r}} \cos \theta - \boldsymbol {\widehat{g}}_{\rm{ext}} \right)k, \ \ \text{ where } k = \end{eqnarray*} (A1) \begin{eqnarray*} \left[ \cos \theta - \frac{4}{5} \sin ^2 \theta \cos \left( \pi q_{\rm{LG}} \right) \frac{\tilde{r}}{1 + \tilde{r}^2}\right] \frac{GM_pa_0r}{2g_{\rm{ext}}\left( r^3 + {r_t}^3\right)}, \end{eqnarray*} (A2) \begin{eqnarray*} \tilde{r} &\equiv & \frac{r}{r_t}, \end{eqnarray*} (A3) \begin{eqnarray*} r_t &=& \frac{\sqrt{GMa_0}}{Qg_{\rm{ext}}}, \end{eqnarray*} (A4) \begin{eqnarray*} q_{\rm{LG}} &=& \frac{M_{\rm{LG}}}{\underbrace{M_{\rm{LG}} + M_p}_M}, \end{eqnarray*} (A5) \begin{eqnarray*} \cos \theta &=& \widehat{\boldsymbol {r}} \cdot \boldsymbol {\widehat{g}}_{\rm{ext}}\ ,, \ \ \sin ^2 \theta = 1 - \cos ^2 \theta, \end{eqnarray*} (A6) \begin{eqnarray*} \boldsymbol {r} &=& \boldsymbol {r}_p - \boldsymbol {r}_i\ ,, \ \ \text{ where $i = $ MW or M31}. \end{eqnarray*} (A7) Here, M refers to the combination of the perturber mass Mp and the LG mass MLG. For the case qLG = 0 or 1, equation (A2) is a fit to numerical results obtained for a point mass embedded in a constant EF. This is a problem we solved previously using our ring library procedure (Banik & Zhao 2018a, section 4.3). The factor of cos (⁠|$\pi$|qLG) is our guess for how the tangential force behaves at other values of qLG, which we cannot solve directly due to the much higher computational cost associated with non-axisymmetric problems. Other functions could also work if they are normalized to be 1 for qLG = 0 and are antisymmetric with respect to qLG → 1 − qLG. APPENDIX B THE TIDAL STRESS TENSOR IN A DOMINANT EXTERNAL FIELD In Section 5.1.2, we discussed how tides raised by perturbers like Cen A could have influenced the MW–M31 trajectory. To understand this analytically, we consider the tidal stress tensor due to a point mass embedded in a dominant external field |$\boldsymbol {g}_{\rm{ext}}$|⁠. To facilitate this analysis, we define Cartesian co-ordinates x, y, z centred on the perturber of mass M and based on axes better suited to the geometry of the problem. \begin{eqnarray*} \widehat{\boldsymbol {z}} \ &=&\ \widehat{\boldsymbol {g}}_{\rm{ext}}, \end{eqnarray*} (B1) \begin{eqnarray*} \widehat{\boldsymbol {y}} \ &\propto &\ \boldsymbol {r} \times \boldsymbol {g}_{\rm{ext}}, \end{eqnarray*} (B2) \begin{eqnarray*} \widehat{\boldsymbol {x}} \ &=&\ \widehat{\boldsymbol {y}} \times \widehat{\boldsymbol {z}}. \end{eqnarray*} (B3) |$\boldsymbol {r}$| is the position of some point of interest relative to the perturber, whose potential is Φ. Due to axisymmetry and the curl-free nature of the gravitational field (∇ × ∇Φ = 0) required to prevent energy gain for a particle going around a closed loop, we have that \begin{eqnarray*} \Phi _{xy} \ &=&\ \Phi _{yx} \ =\ 0, \end{eqnarray*} (B4) \begin{eqnarray*} \Phi _{yz} \ &=&\ \Phi _{zy} \ =\ 0, \end{eqnarray*} (B5) \begin{eqnarray*} \Phi _{yy} \ &=&\ \frac{\Phi _{x}}{x}. \end{eqnarray*} (B6) Here, we use a common shorthand for derivatives according to which, e.g. |$\Phi _{yz} \equiv \frac{\mathrm{\partial} ^2 \Phi }{\mathrm{\partial} x \mathrm{\partial} y}$|⁠. The second spatial derivatives of Φ govern tides raised by a perturber. We approximate that it is sufficiently distant from the LG that the problem is EF dominated, reducing the potential to that given in equation (37) of Banik & Zhao (2015). \begin{eqnarray*} \Phi \ &=&\ -\frac{GM \nu _{\rm{ext}}}{r} \left( 1 + \frac{K_0}{2} \sin ^2 \theta \right), \end{eqnarray*} (B7) \begin{eqnarray*} K_0 \ &\equiv &\ \frac{\mathrm{\partial} {\rm{ln}} \ \nu }{\mathrm{\partial} {\rm{ln}} \ g_{\rm{N,ext}}}. \end{eqnarray*} (B8) This leads to the components of the tidal stress tensor given below, which are not obvious from symmetry arguments like equation (B6). \begin{eqnarray*} \frac{\Phi _{xx}}{\alpha } \ &=&\ \left(K_0 - 1 \right)\left(1 - 3\sin ^2 \theta \right) + \frac{3}{2}K_0 \sin ^2 \theta \left( 5 \sin ^2 \theta - 3\right), \nonumber \\\frac{\Phi _{zz}}{\alpha } \ &=&\ 3\cos ^2 \theta - 1 + \frac{3}{2}K_0 \sin ^2 \theta \left( 5 \cos ^2 \theta - 1\right), \nonumber \\\frac{\Phi _{xz}}{\alpha } \ &=&\ \Phi _{zx} \ =\ \frac{15}{2}K_0 \sin ^3 \theta \cos \theta - 3 \left( K_0 - 1 \right) \sin \theta \cos \theta, \nonumber \\\alpha &\equiv & -\frac{GM \nu _{\rm{ext}}}{r^3}. \end{eqnarray*} (B9) In general, diagonalizing the tidal stress tensor requires us to find the eigenvalues of a symmetric 3 × 3 matrix. However, the alignment of our coordinate system with the axisymmetry of the problem means that the stress tensor is already in block diagonal form. Thus, we only need to solve a 2 × 2 matrix and use equation (B6) to obtain the third eigenvalue λϕ corresponding to the eigenvector |$\widehat{\boldsymbol {y}}_{\rm{eff}}$|⁠. \begin{eqnarray*} \lambda _\phi \ &=&\ \frac{GM \nu _{\rm{ext}}}{r^3} \left( 1 - K_0 + \frac{3}{2} K_0 \sin ^2 \theta \right), \nonumber \\2 \lambda _r \ &=&\ \overbrace{\Phi _{xx} + \Phi _{zz}}^{\alpha \left(K_0 + 1 - \frac{3}{2}K_0 \sin ^2 \theta \right)} - \sqrt{4{\Phi _{xz}}^2 + \left(\Phi _{xx} - \Phi _{zz} \right)^2}, \nonumber \\2 \lambda _\theta \ &=&\ \Phi _{xx} + \Phi _{zz} + \sqrt{4{\Phi _{xz}}^2 + \left(\Phi _{xx} - \Phi _{zz} \right)^2}. \end{eqnarray*} (B10) The eigenvalues λi have been named according to the direction of the corresponding eigenvector |$\hat{\boldsymbol {e}}_i$| in the Newtonian limit (K0 = 0). Note that λr < 0 while the other eigenvalues are positive. \begin{eqnarray*} \hat{\boldsymbol {e}}_\phi &=& \widehat{\boldsymbol {y}}_{\rm {eff}}, \end{eqnarray*} (B11) \begin{eqnarray*} \hat{\boldsymbol {e}}_r &\propto & \widehat{\boldsymbol {x}}_{\rm {eff}} + \frac{\lambda _r - \Phi _{xx}}{\Phi _{xz}} \widehat{\boldsymbol {z}}_{\rm{eff}}, \end{eqnarray*} (B12) \begin{eqnarray*} \hat{\boldsymbol {e}}_\theta &\propto & \widehat{\boldsymbol {x}}_{\rm{eff}} + \frac{\lambda _\theta - \Phi _{xx}}{\Phi _{xz}} \widehat{\boldsymbol {z}}_{\rm{eff}}. \end{eqnarray*} (B13) To help visualize our results, we show the angle between the radial eigenvector |$\hat{\boldsymbol {e}}_r$| and the radial direction |$\hat{\boldsymbol {r}}$|⁠, with positive values indicating a larger angle to |$\widehat{\boldsymbol {g}}_{{ext}}$| (the top panel of Fig. B1). The bottom panel of this figure shows |$\beta \equiv -\frac{2\lambda _r}{\lambda _\theta + \lambda _\phi }$|to help give an idea about the tidal stress tensor’s shape. As the potential is symmetric with respect to |$\boldsymbol {g}_{\rm{ext}} \rightarrow -\boldsymbol {g}_{\rm{ext}}$|⁠, we only show results for |$\theta \le \frac{\pi }{2}$|⁠, i.e. |$\boldsymbol {r} \cdot \boldsymbol {g}_{\rm{ext}} >0$|⁠. Figure B1. Open in new tabDownload slide Top: Angle between the radial direction and the most radial eigenvector of the tidal stress tensor due to a point mass embedded in a dominant external field. Positive values indicate that |$\hat{\boldsymbol {e}}_r \cdot \widehat{\boldsymbol {g}}_{\rm{ext}} < \hat{\boldsymbol {r}} \cdot \widehat{\boldsymbol {g}}_{\rm{ext}}$| i.e. the radial eigenvector is tilted further away from |$\widehat{\boldsymbol {g}}_{\rm{ext}}$| than the radial direction. K0 is the logarithmic derivative of the ν function, with value 0 in Newtonian gravity and |$-\frac{1}{2}$| in the deep-MOND limit. Bottom: Values of |$-\frac{2\lambda _r}{\lambda _\theta + \lambda _\phi }$|⁠, giving an idea about the tidal stress tensor’s shape. Figure B1. Open in new tabDownload slide Top: Angle between the radial direction and the most radial eigenvector of the tidal stress tensor due to a point mass embedded in a dominant external field. Positive values indicate that |$\hat{\boldsymbol {e}}_r \cdot \widehat{\boldsymbol {g}}_{\rm{ext}} < \hat{\boldsymbol {r}} \cdot \widehat{\boldsymbol {g}}_{\rm{ext}}$| i.e. the radial eigenvector is tilted further away from |$\widehat{\boldsymbol {g}}_{\rm{ext}}$| than the radial direction. K0 is the logarithmic derivative of the ν function, with value 0 in Newtonian gravity and |$-\frac{1}{2}$| in the deep-MOND limit. Bottom: Values of |$-\frac{2\lambda _r}{\lambda _\theta + \lambda _\phi }$|⁠, giving an idea about the tidal stress tensor’s shape. In Newtonian gravity, β = 2 and the non-radial directions are equivalent (λθ = λϕ). In general, this is not true as there is an additional preferred direction due to the EF. Moreover, β can also differ from 2 as the gravitational field need not be divergence-free outside the matter distribution (equation 4). For two particles separated by a small amount |$\boldsymbol {s}$| compared to their distance r from the perturber, we can now find the tidal contribution to the acceleration of their separation |${\ddot{s}}$|⁠. This can be decomposed as a sum over contributions from the three eigenvalues λi which depend on how well |$\boldsymbol {s}$| aligns with the corresponding eigenvector |$\hat{\boldsymbol {e}}_{i}$|⁠. \begin{eqnarray*} \left. \frac{\ddot{s}}{s} \right|_{\rm {tides}} \ =\ -\sum _{i = 1}^3 \lambda _i \left( \hat{\boldsymbol {s}} \cdot \hat{\boldsymbol {e}}_{i} \right)^2 \end{eqnarray*} (B14) The actual value of |${\ddot{s}}$| will include contributions from other sources like the cosmological acceleration. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)
TI - Origin of the Local Group satellite planes
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty919
DA - 2018-07-11
UR - https://www.deepdyve.com/lp/oxford-university-press/origin-of-the-local-group-satellite-planes-TFYIWfJYLa
SP - 4768
VL - 477
IS - 4
DP - DeepDyve
ER -