pdftexdestination with the same hyperrefOption hyperrefToken pdftex(dest)

# Origin of the Local Group satellite planes

Indranil Banik1, David O’Ryan1 and Hongsheng Zhao1
1Scottish Universities Physics Alliance, University of St Andrews, North
Haugh, St Andrews, Fife, KY16 9SS, UK Email<EMAIL_ADDRESS>(Indranil
Banik)
<EMAIL_ADDRESS>(Hongsheng Zhao)

###### Abstract

We attempt to understand the planes of satellite galaxies orbiting the Milky
Way (MW) and M31 in the context of Modified Newtonian Dynamics (MOND), which
implies a close MW-M31 flyby occurred ${\approx 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 disks, 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 disk 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 disk. Thus, the MW thick disk may have formed
together with the MW and M31 satellite planes during a past MW-M31 flyby.

###### keywords:

galaxies: groups: individual: Local Group – Galaxy: kinematics and dynamics –
Dark Matter – methods: numerical – methods: data analysis – Galaxy: disc

††pubyear: 2018††pagerange: Origin of the Local Group satellite planes–B

##  1 Introduction

The standard cosmological paradigm ($\Lambda$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 et al. (2005) to argue that its satellite system is
inconsistent with $\Lambda$CDM as it is too anisotropic. Later proper motion
measurements showed that most of its satellites co-rotate within a well-
defined plane (Metz et al., 2008; Pawlowski & Kroupa, 2013). Recently
discovered ultra-faint satellites, globular clusters and tidal streams
independently prefer a similarly oriented plane (Pawlowski et al., 2012;
Pawlowski & Kroupa, 2014). Although some flattening is expected in
$\Lambda$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 et al., 2017).

An analogous situation was suspected around M31 (Metz et al., 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
$\Lambda$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 et al. (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 co-rotating
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 $\Lambda$CDM paradigm (González & Padilla, 2010, figure 8).

Although $\Lambda$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 (Kroupa, 2015;
Javanmardi et al., 2015; Javanmardi & Kroupa, 2017; Bengaly et al., 2018). The
model also has difficulty explaining the strength of the 21 cm absorption
feature corresponding to the hyperfine transition of cold neutral hydrogen at
redshift ${z\approx 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).111The particle mass would
have to be $10-80$ MeV/$c^{2}$ 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. 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
reionisation (Heinrich & Hu, 2018). For example, the black holes required by
their model would need to stop growing at ${z\approx 16}$.

Given these difficulties and the many a priori predictions of 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 $\Lambda$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,
$\Lambda$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 $\Lambda$CDM based on the
scatter about the Hubble flow in detailed $N$-body simulations of it (Aragon-
Calvo et al., 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, figure 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 and the Great Attractor.

To test whether the issue might be resolved using a three-dimensional (3D)
model of the LG, we borrowed an algorithm described in Phelps et al. (2013).
The typical mismatch between observed and predicted GRVs in the best-fitting
model is actually slightly higher than in the 2D case, with a clear tendency
persisting for faster outward motion than expected (Banik & Zhao, 2017,
figures 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 $\Lambda$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 ${\,{\sim}\,3}$ Mpc, with an isotropic
distribution ruled out at ${3\sigma}$ (see their table 7).

$\Lambda$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 Modified Newtonian
Dynamics (MOND, Milgrom, 1983). In MOND, galaxies are not surrounded by halos
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

$\displaystyle g~{}=~{}\frac{\sqrt{GMa_{{}_{0}}}}{r}~{}~{}~{}\text{for
}~{}r\gg\sqrt{\frac{GM}{a_{{}_{0}}}}$ (1)

MOND introduces Milgrom’s constant $a_{{}_{0}}$ as a fundamental acceleration
scale of nature below which the deviation from Newtonian dynamics becomes
significant. Empirically, $a_{{}_{0}}\approx 1.2\times{10}^{-10}$ m/s2 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_{{}_{\Lambda}}=\rho_{{}_{\Lambda}}c^{2}$ implied by the accelerating
expansion of the Universe (Riess et al., 1998). Thus,

$\displaystyle\frac{g^{2}}{8\rm{\pi}G}~{}<~{}u_{{}_{\Lambda}}~{}~{}\Leftrightarrow~{}~{}g~{}\la~{}2\rm{\pi}a_{{}_{0}}$
(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 emphasising 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 halos of the $\Lambda$CDM paradigm (Salucci &
Turini, 2017; Desmond, 2017a, b).

In the LG, MOND implies a much stronger MW-M31 mutual attraction than does
$\Lambda$CDM. Acting on their almost radial orbit (van der Marel et al.,
2012), this leads to a close encounter ${9\pm 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^{\circ}}$ 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 et al., 1992) and in MOND
simulations of them (Tiret & Combes, 2008; Renaud et al., 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 co-rotating within that plane
(Wetzstein et al., 2007), though a small fraction might well end up counter-
rotating depending on the exact details (Pawlowski et al., 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 $\Lambda$CDM as dynamical friction between their overlapping dark
matter halos 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,
$\Lambda$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 $\Lambda$CDM context (see their figure 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 halos. 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 $\Lambda$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
$\Lambda$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\sigma}$ once the survey footprint is taken into
account (Pawlowski, 2016).

More recent hydrodynamical $\Lambda$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-wide plane 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 $\Lambda$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 halos,
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
disks. 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 disk is dynamically heated by its interaction
with M31. This is to check the feasibility of our scenario for producing the
MW thick disk (Gilmore & Reid, 1983), a structure which seems to have formed
fairly rapidly from its thin disk (Hayden et al., 2015). The thick disk formed
${9\pm 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, figure 2). An
enhanced star formation rate could have further thickened the MW disk through
the popping star cluster scenario (Kroupa, 2002). The star formation rate of
M31 also appears to have been much higher for lookback times ${\ga 9}$ Gyr
(Williams et al., 2017b, figure 12). The disk heating which likely formed the
Galactic thick disk 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 disk has a longer scale length than its thin disk (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 disks 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 $\chi^{2}$
statistic (Equation 40), 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 Milky Way $-$ 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
quasilinear 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 $\la 5\%$ (Banik & Zhao, 2015). $N$-body simulations of
both MOND formulations yield quite similar results (Candlish, 2016).

QUMOND uses the Newtonian gravitational field $\bm{g}_{{}_{N}}$ to determine
the true gravitational field $\bm{g}$. Thus, the first step is to find
$\bm{g}_{{}_{N}}$ by solving the usual Poisson equation sourced by the
baryonic density $\rho_{{}_{b}}$.

$\displaystyle\nabla\cdot\bm{g}_{{}_{N}}~{}=~{}-4\pi G\rho_{{}_{b}}$ (3)

To determine $\bm{g}$, this Poisson equation must be solved again with a
forcing analytically dependent on $\bm{g}_{{}_{N}}$. Thus, QUMOND does not
require a non-linear grid relaxation stage, a major improvement on AQUAL.

$\displaystyle\overbrace{\nabla\cdot\bm{g}}^{\propto\rho_{{}_{PDM}}+\rho_{{}_{b}}}~{}$
$\displaystyle=$
$\displaystyle~{}\nabla\cdot\left[\bm{g}_{{}_{N}}\nu\left(g_{{}_{N}}\right)\right]~{}~{}\text{
where}$ (4) $\displaystyle\nu\left(g_{{}_{N}}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2}~{}+~{}\sqrt{\frac{1}{4}+\frac{a_{{}_{0}}}{g_{{}_{N}}}}$
(5)

$\nu\left(g_{{}_{N}}\right)$ is the interpolating function used to transition
between the Newtonian and deep-MOND regimes. Here, we use the ‘simple’ form of
this function (Famaey & Binney, 2005). This works very well at explaining the
rotation curve of the MW (Iocco et al., 2015, figure 3) and the kinematics of
${\approx 6000}$ elliptical galaxies probing accelerations up to ${\approx
30a_{{}_{0}}}$ (Chae et al., 2017). The source term for the gravitational
field is $\nabla\cdot\left(\nu\bm{g}_{{}_{N}}\right)$, which can be thought of
as an ‘effective’ density $\rho$ composed of the actual density
$\rho_{{}_{b}}$ and an extra term which we define to be the phantom dark
matter density $\rho_{{}_{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 $d_{{}_{M31}}$,
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}_{{}_{M31}}$ is a significant fraction of $a_{{}_{0}}$ (Section
2.1.2). If the MW and M31 were the only objects in the Universe, then Zhao et
al. (2010) showed that we could approximate their mutual gravitational
acceleration as

$\displaystyle g_{iso}~{}$ $\displaystyle=$
$\displaystyle~{}\frac{Q\sqrt{GM_{{}_{LG}}a_{0}}}{d_{{}_{M31}}}$ (6)
$\displaystyle Q~{}$ $\displaystyle=$
$\displaystyle~{}\frac{2\left(1-{q_{{}_{1}}}^{\frac{3}{2}}-{q_{{}_{2}}}^{\frac{3}{2}}\right)}{3q_{{}_{1}}q_{{}_{2}}}$
(7)

Here, the total MW and M31 mass is $M_{{}_{LG}}$, of which a fraction
$q_{{}_{1}}$ is in the MW and the remaining $q_{{}_{2}}$ resides in M31.
Equation 6 suggests that $\ddot{d}_{{}_{M31}}{\approx 0.017~{}a_{{}_{0}}}$.
However, the MOND field equations imply that $\ddot{d}_{{}_{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 $g_{{}_{M}}$ between the MW and M31 can be written as

$\displaystyle g_{{}_{M}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{GMa_{0}}}{r}f\left(\tilde{r}\right)$ (8)
$\displaystyle\tilde{r}$ $\displaystyle\equiv$
$\displaystyle\frac{d_{{}_{M31}}}{r_{t}}$ (9) $\displaystyle r_{t}$
$\displaystyle=$
$\displaystyle\frac{\sqrt{GMa_{0}}\left(3+\cos^{2}\phi\right)}{4Qg_{ext}}$
(10)

We have assumed that the form of $f\left(\tilde{r}\right)$ is independent of
the angle $\phi$ between the MW-M31 axis and $\widehat{\bm{g}}_{{}_{ext}}$,
with the angular dependence introduced by the EF accounted for by making
$r_{t}$ angle-dependent. The function $f\left(\tilde{r}\right)$ should be $Q$
when $\tilde{r}\ll 1$ or equivalently when $g_{{}_{ext}}\ll g_{iso}$,
representing the isolated regime (Zhao et al., 2010). In the opposite EF-
dominated case, $f\to\frac{Q}{\tilde{r}}$ but the dependence on $Q$ is
cancelled by how this also affects $r_{t}$ (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 $g_{ext}\approx 0.03~{}a_{{}_{0}}$ (Famaey et al.,
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 $g_{iso}$ (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 $\bm{g}_{{}_{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 $\bm{g}_{{}_{N,ext}}$ to the forces from
the two point masses. $\bm{g}_{{}_{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 $\bm{g}_{{}_{N,ext}}$ and the
true EF of $\bm{g}_{{}_{ext}}$.

$\displaystyle\overbrace{\nu\left(\frac{\left|\bm{g}_{{}_{N,ext}}\right|}{a_{{}_{0}}}\right)}^{\nu_{{}_{ext}}}\bm{g}_{{}_{N,ext}}$
$\displaystyle=$ $\displaystyle\bm{g}_{ext}$ (11)

In the deep-MOND limit, this reduces to

$\displaystyle g_{N,ext}~{}=~{}\frac{{g_{ext}}^{2}}{a_{{}_{0}}}$ (12)

Having obtained $\bm{g}_{{}_{N}}$, we can use Equation 4 to obtain
$\nabla\cdot\bm{g}$, the effective matter distribution. We apply a direct
summation procedure to this in order to obtain $\bm{g}$.

$\displaystyle\bm{g\left(\bm{r}\right)}~{}=~{}\int\nabla\cdot\bm{g}\left(\bm{r^{\prime}}\right)\frac{\left(\bm{r}-\bm{r^{\prime}}\right)}{4\pi|\bm{r}-\bm{r^{\prime}}|^{3}}~{}d^{3}\bm{r^{\prime}}$
(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 $\bm{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 55 to estimate that the mass ratio between the MW and M31 is
such that ${q_{{}_{1}}=0.3}$ based on the flatline level of the MW rotation
curve being ${180}$ km/s (Kafle et al., 2012) and that of M31 being ${225}$
km/s (Carignan et al., 2006). With ${q_{{}_{1}}}$ fixed as we vary
$g_{{}_{ext}}$, our numerical results can be summarised as

$\displaystyle
f\left(\tilde{r}\right)=\frac{Q}{\sqrt{1+\tilde{r}^{2}}}\left[1+\overbrace{0.032~{}\rm{exp}\left({-\frac{\left(Ln~{}\tilde{r}-0.42\right)^{2}}{2\times
0.65^{2}}}\right)}^{u}\right]$ (14)

We determine $f\left(\tilde{r}\right)$ by performing 18 force calculations
which vary $\tilde{r}$ over the range $\left(0-12.5\right)$, thus probing a
similar range in $g_{ext}/g_{iso}$. 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
(Figure 1). With 18 force calculations and 3 degrees of freedom in our fitting
function, we obtain a total $\chi^{2}$ of 16.1. Therefore, Equation 14
provides a good fit to our numerical results.

Figure 1: Natural logarithm (Ln) of $f\left(\tilde{r}\right)$ (Equation 8)
against $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 $Ln~{}\tilde{r}=-2.5$ instead of the correct value of
$-\infty$. All forces shown here are determined to within 1%.
Indistinguishable results are obtained if we set
$\bm{g}_{{}_{ext}}\to-\bm{g}_{{}_{ext}}$. We fit these results using Equation
14 (blue line). The total $\chi^{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 $M_{{}_{LG}}$.

$\displaystyle g_{{}_{N}}~{}=~{}\frac{GM_{{}_{LG}}}{{d_{{}_{M31}}}^{2}}$ (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 $d_{{}_{M31}}$ to
the geometric means of the MW-M31 separations at which $g_{{}_{N}}=a_{{}_{0}}$
and where $g_{{}_{N}}=g_{{}_{N,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 $g_{{}_{N}}<g_{{}_{ext}}$. If it is, then the EF is
more important. As $g_{{}_{ext}}\ll a_{{}_{0}}$, it is not as important to
consider violation of the deep-MOND limit. In this case, we therefore use our
previously derived result for $g_{{}_{M}}$ including the EF but assuming the
deep-MOND limit (Equation 8).

However, if $g_{{}_{N}}>g_{{}_{ext}}$, it is not as important to consider the
EF because $g_{{}_{ext}}\gg g_{{}_{N,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}_{{}_{M31}}$ including both the EF and a variable MW-M31 separation.
This can only be done by having a two-dimensional 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 $g_{{}_{N}}=g_{{}_{N,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 $a_{{}_{0}}$. Thus,
if the EF is important, then by definition we must have that $g\la
g_{{}_{ext}}$. This implies that $g\ll a_{{}_{0}}$, 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 $g_{ext}\left(t\right)$ must explain the present peculiar
velocity $\bm{v}_{pec}$ of the LG i.e. its deviation from the Hubble flow.

$\displaystyle\bm{v}_{pec}~{}\equiv~{}\dot{\bm{r}}_{{}_{LG}}-H\bm{r}_{{}_{LG}}$
(16)

Here, $\bm{r}_{{}_{LG}}$ is the position of the LG relative to a co-moving
observer. We use an overlying dot to indicate a time derivative (e.g.
$\dot{q}\equiv\frac{\partial q}{\partial t}$ for any quantity $q$). The Hubble
parameter $H\equiv\frac{\dot{a}}{a}$ where $a\left(t\right)$ is the cosmic
scale-factor relative to its present value (i.e. $a\equiv 1$ at the current
time $t=t_{f}$).

For the LG, $\bm{v}_{pec}$ has a value close to 630 km/s towards Galactic co-
ordinates $\left(276^{\circ},30^{\circ}\right)$ 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
$\bm{\widehat{g}}_{ext}\equiv\bm{{g}}_{ext}\div\left|\bm{{g}}_{ext}\right|$.

Peculiar velocities must be caused by the time-integrated effect of
$\bm{g}_{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

$\displaystyle\int_{0}^{t_{f}}a~{}g_{ext}~{}dt~{}=~{}v_{pec}$ (17)

The integrating factor $a\left(t\right)$ 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).222$v_{pec}\propto\frac{1}{a\left(t\right)}$ for a test particle in an
otherwise homogeneous universe, making $av_{pec}$ constant. There is no unique
way to solve for $g_{ext}\left(t\right)$, but neither is it wholly
unconstrained. To make further progress, we assume that $g_{ext}$ is due to a
very distant point mass. Thus,

$\displaystyle g_{ext}\propto\left\\{\begin{array}[]{ll}\frac{1}{a}&\mbox{if
}a\geq a_{peak}\\\ \left(a-a_{c}\right)&\mbox{if }a_{c}<a<a_{peak}\\\
0&\mbox{if }a\leq a_{c}\end{array}\right.$ (21)

It is not possible that $g_{ext}\propto 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<a_{c}=0.1$), $g_{ext}=0$. To avoid a discontinuous jump in
$g_{ext}$, we assume that $g_{ext}~{}\propto~{}\left(a-a_{c}\right)$ when
$a_{c}<a<a_{peak}$.

To transform $g_{ext}\left(a\right)$ into $g_{ext}\left(t\right)$ and thus
solve Equation 17, we need a specific form for the background cosmology
$a\left(t\right)$. We use a standard
flat333$\Omega_{m,0}+\Omega_{\Lambda,0}=1$ dark energy-dominated cosmology
with parameters given in Table 1.

$\displaystyle\frac{\overset{..}{a}}{a}$ $\displaystyle=$
$\displaystyle-\frac{4\rm{\pi}G}{3}\left(\rho_{m}-2\rho_{\Lambda}\right)$ (22)
$\displaystyle=$
$\displaystyle{H_{{}_{0}}}^{2}\left(-\frac{1}{2}\Omega_{m,0}~{}a^{-3}+\Omega_{\Lambda,0}\right)$
(23)

Here, $\rho_{m}$ and $\rho_{\Lambda}$ are the cosmic mean densities of matter
and dark energy, respectively. Currently, these components comprise fractions
$\Omega_{m,0}$ and $\Omega_{\Lambda,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

$\displaystyle a\left(t\right)$ $\displaystyle=$
$\displaystyle{{\left(\frac{{{\Omega}_{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)$
(24)

Name | Meaning and units | Value
---|---|---
$\Omega_{m,0}$ | Fraction of critical density in matter | 0.315
$\Omega_{\Lambda,0}$ | Fraction of critical density in dark energy | 0.685
$H_{{}_{0}}$ | Hubble constant, km/s/Mpc | 67.3
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 24).

The present values of $H_{{}_{0}}$ and $\Omega_{m,0}$ uniquely determine the
present age of the Universe $t_{f}$ 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 24 is valid for $\Lambda$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 $\Lambda$CDM (Zhao, 2007).

In our model, $g_{ext}$ peaks when $a=a_{peak}$ at a value higher than the
present value ($g_{ext,0}$) by the factor $\frac{1}{a_{peak}}$. We consider
models to be realistic if $a_{peak}\approx 0.77$, when the Universe stopped
being matter-dominated. This requires $g_{ext,0}\approx 0.023~{}a_{0}$, 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 $\bm{r}_{{}_{LG}}$ and a perturber at position
$\bm{r}_{{}_{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 $\bm{r}_{{}_{LG}}$. The distance $d_{p}$ of
the perturber from $\bm{r}_{{}_{LG}}$ is thus the only dynamical quantity. We
vary its initial value at some early time $t_{{}_{i}}$ 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.

$\displaystyle\ddot{d}_{p}~{}$ $\displaystyle=$
$\displaystyle~{}\frac{\ddot{a}}{a}d_{p}~{}+~{}g_{{}_{M}}$ (25)
$\displaystyle\dot{d}_{p}~{}$ $\displaystyle=$
$\displaystyle~{}Hd_{p}~{}~{}\text{initially}$ (26)

Here, $g_{{}_{M}}$ 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 (${\ga 3}$ Mpc) perturbers.

Name | $b$ | $l$ | $d_{MW}$ | $M$
---|---|---|---|---
Centaurus A | $19.4173^{\circ}$ | $309.5159^{\circ}$ | 3.8 | 6
M81 | $40.9001^{\circ}$ | $142.0918^{\circ}$ | 3.6 | 2.2
IC 342 | $10.5799^{\circ}$ | $138.1726^{\circ}$ | 3.45 | 2.2
Table 2: Properties of mass concentrations outside the Local Group which we
consider. Sky positions are shown using Galactic co-ordinates. The
Galactocentric distances $d_{MW}$ are given in Mpc for Cen A (Harris et al.,
2010), M81 (Gerke et al., 2011) and IC 342 (Wu et al., 2014). Baryonic masses
$M$ in units of $10^{11}M_{\odot}$ are from Karachentsev (2005) for Cen A and
IC 342 while Karachentsev & Kashibadze (2006) is used for M81.

Although we numerically calibrate the form of $g_{{}_{M}}$ 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 $g_{{}_{M}}$, 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_{{}_{iso}}\approx\frac{\sqrt{GMa_{{}_{0}}}}{d_{p}}$ where $M$ is the total
(LG $+$ Cen A) mass. For a distance to Cen A of $d_{p}=3.8$ Mpc, this yields
$0.0070a_{{}_{0}}$. As the actual EF on the LG is very likely larger (Section
2.2), the LG-Cen A problem can’t be treated in isolation.

If instead we assume that the EF from more distant objects dominates and set
$g_{{}_{ext}}=0.02a_{{}_{0}}$, then we get that Cen A pulls on the LG by
$g_{{}_{EFE}}\approx\frac{GMa_{{}_{0}}}{{d_{p}}^{2}g_{{}_{ext}}}$ i.e. by only
$0.0024a_{{}_{0}}$. As this is much less than $g_{{}_{ext}}$, 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 $d_{{}_{M31}}\left(t\right)$. 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 $g_{{}_{M}}$ 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.

$g_{{}_{M}}=$
$\displaystyle\frac{\sqrt{GMa_{0}}}{d_{p}}f\left(\tilde{r}\right)\left(1+\frac{q_{{}_{MW}}q_{{}_{M31}}n\left(n\cos^{2}\theta+1\right){d_{{}_{M31}}}^{2}}{2{d_{p}}^{2}}\right)$
(27) $\displaystyle n$ $\displaystyle\equiv$ $\displaystyle
1~{}-~{}\frac{\partial
Ln~{}f\left(\tilde{r}\right)}{\partial\tilde{r}}~{}~{}\left(g_{{}_{M}}\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}{\tilde{r}}^{-n}\right)$ (30) $\displaystyle=$
$\displaystyle
1~{}+~{}\frac{\tilde{r}^{2}}{\tilde{r}^{2}+1}~{}+~{}\frac{\left(Ln~{}\tilde{r}-0.42\right)u}{0.65^{2}\left(1+u\right)}$
(31) $\displaystyle\tilde{r}$ $\displaystyle=$
$\displaystyle\frac{4d_{p}Qg_{ext}}{\sqrt{GMa_{0}}\left(3+\left(\widehat{\bm{r}}\cdot\bm{\widehat{g}}_{ext}\right)^{2}\right)}$
(32)

Here, $M$ is the total mass of the LG and the perturber while $Q$ is
calculated for the ratio between their masses. $q_{{}_{MW}}$ and
$q_{{}_{M31}}$ 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 $g_{{}_{M}}$ found here depends on the angle $\theta$
between the MW$-$M31 line and the direction from the LG barycentre to the
perturber.

Using the long bracketed term in Equation 27 improves the accuracy of
$\ddot{d}_{{}_{p}}$ and allows us to refine their trajectory. We only do this
once as it leads to fairly small corrections. The final perturber trajectories
$d_{p}\left(t\right)$ are shown in Figure 2 for our best-fitting model
(Section 5).

Figure 2: Trajectories of the perturbers in Table 2, i.e. their distance from
the LG barycentre $\bm{r}_{{}_{LG}}$ over time. The MW is shown here for
comparison (M31 is much closer to $\bm{r}_{{}_{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 onto
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, figure 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
$\bm{r}_{{}_{LG}}=\bm{r}_{{}_{MW}}$ or $\bm{r}_{{}_{LG}}=\bm{r}_{{}_{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 $\bm{g}_{tan}$ between two objects
varies with their mass ratio (Appendix A).

In this way, we estimate the gravity $\bm{g}_{{}_{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
onto the line connecting them so as to estimate $g_{{}_{tide}}$, our symbol
for how much each perturber tidally alters the MW-M31 mutual attraction
$\ddot{d}_{{}_{M31}}$.

$\displaystyle g_{{}_{tide}}$ $\displaystyle=$
$\displaystyle\left(\left.\bm{g}_{{}_{LG}}\right|_{\bm{r}_{{}_{LG}}=\bm{r}_{{}_{M31}}}-\left.\bm{g}_{{}_{LG}}\right|_{\bm{r}_{{}_{LG}}=\bm{r}_{{}_{MW}}}\right)\cdot\widehat{\bm{d}}_{{}_{M31}}$
$\displaystyle\widehat{\bm{d}}_{{}_{M31}}$ $\displaystyle=$
$\displaystyle\frac{\bm{r}_{{}_{M31}}-\bm{r}_{{}_{MW}}}{\left|\bm{r}_{{}_{M31}}-\bm{r}_{{}_{MW}}\right|}$
(33)

We assume that the EF on the LG arises due to a mass concentration similar to
the Great Attractor (GA, Mieske et al., 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\left(t\right)$ 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, figure 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

$\displaystyle
g_{{}_{tide,GA}}~{}=~{}\frac{\left(2\cos^{2}\theta_{ext}-1\right)d_{{}_{M31}}g_{{}_{ext}}}{d_{{}_{GA}}}$
(34)

Here, $\theta_{ext}$ is the angle between the MW-M31 line and
$\bm{\widehat{g}}_{ext}$. The factor of $\left(2\cos^{2}\theta_{ext}-1\right)$
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 can’t be approximated as radial if they undergo a past
close encounter. For simplicity, we assume that their mutual orbital angular
momentum ${\bm{h}}_{{}_{LG}}$ has remained constant since their first
turnaround when the Universe was ${\,{\sim}\,3}$ Gyr old. This is because
${\bm{h}}_{{}_{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, figure 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 ${\bm{h}}_{{}_{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 Large Magellanic Cloud (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{\bm{h}}_{{}_{LG}}$ as ${\bm{h}}_{{}_{LG}}$ is
so small.444Mergers 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 halos in $\Lambda$CDM, mergers would be much rarer
in MOND, though its stronger long-range gravity would lead to more
gravitational focusing. In this case, our best-fitting model would have the
same ${\bm{h}}_{{}_{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 ${\bm{h}}_{{}_{LG}}$ has
changed since the MW-M31 flyby. Instead, we consider all possible orientations
for $\widehat{\bm{h}}_{{}_{LG}}$ and a wide range of values for
$\left|{\bm{h}}_{{}_{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}_{{}_{M31}}$ from cosmological acceleration
($\frac{\ddot{a}}{a}d_{{}_{M31}}$), angular momentum
(${{{h_{{}_{LG}}}^{2}}/{{d_{{}_{M31}}}^{3}}}$), self-gravity $g$ (Section 2.1)
and the tidal contribution $g_{{}_{tide}}$ due to each perturber (Equations
2.3.2 and 34), we get that

$\displaystyle\ddot{d}_{{}_{M31}}~{}=~{}\frac{\ddot{a}}{a}d_{{}_{M31}}~{}+~{}\frac{{h_{{}_{LG}}}^{2}}{{d_{{}_{M31}}}^{3}}~{}+~{}g~{}+~{}\sum_{Perturbers}g_{{}_{tide}}$
(35)

We use this to integrate the MW-M31 trajectory back from present conditions,
assuming a particular MW-M31 specific angular momentum $\bm{h}_{{}_{LG}}$ in
the period after their first turnaround but a purely radial orbit before then.
As the MW-M31 line $\widehat{\bm{d}}_{{}_{M31}}$ rotates, it changes
orientation relative to the EF and the directions towards nearby perturbers.
Thus, the results of our simulations depend on $\widehat{{\bm{h}}}_{{}_{LG}}$.

$\displaystyle\dot{\widehat{\bm{d}}}_{{}_{M31}}~{}=~{}\frac{\bm{h}_{{}_{LG}}\times\widehat{\bm{d}}_{{}_{M31}}}{{{d}_{{}_{M31}}}^{2}}$
(36)

$\widehat{\bm{d}}_{{}_{M31}}$ rotates substantially during our simulation, but
almost all of this rotation occurs close to the time of the MW-M31 flyby as
${d}_{{}_{M31}}$ is rather small then. Thus, it is not too important to
accurately know $\bm{h}_{{}_{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 backwards integration should leave the
MW-M31 separation satisfying

$\displaystyle\dot{d}_{{}_{M31}}~{}=~{}H_{{}_{i}}d_{{}_{M31}}~{}~{}\text{when
}t=t_{i}$ (37)

Figure 3: Three solutions for the MW-M31 separation
$d_{{}_{M31}}\left(t\right)$, 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 37). If $M$
is too low (red curve), $v_{pec}\propto\frac{1}{a}\to+\infty$ while if $M$ is
too high (blue curve), $v_{pec}\to-\infty$ as $a\to 0$. In neither case does
$d_{{}_{M31}}$ smoothly approach 0 at early times.

Our simulations cover the period from $t=t_{{}_{f}}$ (the present epoch) back
to $t=t_{{}_{i}}$, the limit of our backwards integration. We refer to this as
the start time of our simulations, even though in a numerical sense they start
at $t_{{}_{f}}$. We use Equation 24 to determine $t_{{}_{i}}$ by solving it
for when $a=0.05$ (redshift 19). This allows us to determine the Hubble
parameter $H_{{}_{i}}$ at that epoch.

In general, a particular set of model parameters will not satisfy Equation 37.
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 Figure 3, we show how $d_{{}_{M31}}\left(t\right)$ 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 a MW-M31
pericentre at the time of the Big Bang.

##  3 The MW and M31 disks

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 disks. 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 $\bm{g}$ governing our test
particles is axisymmetric, though the symmetry axis
$\widehat{\bm{d}}_{{}_{M31}}$ changes with time according to Equation 36. We
can thus obtain $\bm{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 utilising a ‘ring library’ which stores $\bm{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 $\bm{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 disk scale
length of either galaxy (Table 3, stellar disk used for MW). This should not
affect our results too much because we are mostly interested in their
satellites. These lie several disk scale lengths from their hosts, making a
point mass approximation reasonably accurate for a theory where galaxies lack
extended dark matter halos.

Disk scale length of.. | Value (kpc)
---|---
MW stars | 2.15
MW gas | 7
M31 | 5.3
Table 3: Parameters of the MW and M31 disks. We use kpc for the exponential
scale lengths of the MW stellar disk (Bovy & Rix, 2013), its gas disk
(McMillan, 2017) and the M31 disk (Courteau et al., 2011). We model the MW
disk surface density as a double exponential with 0.1764 of its mass in its
gas disk and the rest in its stellar disk (Banik & Zhao, 2018a, table 1). M31
is also treated as an exponential disk, though its mass is always
$\frac{7}{3}\times$ the total MW mass. Disk orientations are discussed in
Table 5.

Each test particle is initially set up at position $\bm{r}$ on a circular
orbit around galaxy $i$ of mass $M_{i}$ and position $\bm{r}_{{}_{i}}$. The
circular speed is based on a gravitational field directly towards galaxy $i$
of strength

$\displaystyle g$ $\displaystyle=$
$\displaystyle\frac{g_{{}_{N}}}{2}~{}+~{}\sqrt{\left(\frac{g_{{}_{N}}}{2}\right)^{2}+g_{{}_{N}}a_{{}_{0}}}$
(38) $\displaystyle g_{{}_{N}}$ $\displaystyle=$
$\displaystyle\frac{GM_{i}}{\left|\bm{r}-\bm{r}_{{}_{i}}\right|^{2}}$ (39)

We add the systemic velocity of galaxy $i$ to the circular velocity
$v_{{}_{c}}=\sqrt{\left|\bm{r}-\bm{r}_{{}_{i}}\right|g}$. Our estimates of
$v_{{}_{c}}$ 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 disks 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 disks 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 kpc 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 nonetheless 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 disk 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.

Parameter | Value for MW | Value for M31
---|---|---
Radial resolution (kpc) | 0.215 | 0.53
Azimuthal resolution | $2.8125^{\circ}$ | $2.8125^{\circ}$
Inner radial limit (kpc) | 2.15 | 5.3
Outer radial limit (kpc) | 107.5 | 106
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 disk
scale length of each galaxy (using the stellar disk for the MW).

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 timestep is quantised 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 timestep to be below
$\frac{1}{200}$ of the Universe’s age. If necessary, our algorithm can reduce
the timestep 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 timestep is generally used, though this never exceeds 6.8
Myr.

###  3.1 Disk orientations

Due to the MW-M31 flyby, their disks end up in a slightly different
orientation to what they start with. To account for this, we adjust the
initial disk spin vectors of the MW and M31 so as to match their presently
observed orientations. Because our simulations use test particle disks, it is
possible to adjust the disks independently. We use a 2D Newton-Raphson
procedure to match the presently observed orientation of each disk.555We
accelerate this process by assuming a change to the initial disk spin of
either galaxy causes an equal change to its final disk spin (e.g. tilting the
MW $5^{\circ}$ further north initially causes its present spin to also tilt
$5^{\circ}$ 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\times}$.

We determine the final MW disk 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 disk plane (using the present MW disk 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
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$\sigma$, corresponding to the 99% 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 disks, we need
their initial orientations to be as shown in Table 5. Fortunately, both disks
precess by only a few degrees. Thus, we are always sure which hemisphere
contains the total angular momentum of each disk. 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.

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 onto the MW and M31.

Initial disk | Initial | Final
---|---|---
spin vector of.. | direction | direction
Milky Way (MW) | $\begin{bmatrix}86.61^{\circ}\\\ -87.48^{\circ}\end{bmatrix}$ | $\begin{bmatrix}0^{\circ}\\\ -90^{\circ}\end{bmatrix}$
Andromeda (M31) | $\begin{bmatrix}243.35^{\circ}\\\ -23.61^{\circ}\end{bmatrix}$ | $\begin{bmatrix}238.65^{\circ}\\\ -26.89^{\circ}\end{bmatrix}$
Table 5: Initial spin vectors of the MW and M31 disks in Galactic co-
ordinates. These differ slightly from their presently observed orientations
(last column) because the disks precess during their close encounter. The
present M31 disk orientation is obtained from Banik & Zhao (2018b, table 1)
and changes by $5.37^{\circ}$ while that of the MW changes by only
$2.52^{\circ}$.

##  4 Comparison with observations

We compare each model with observations using a $\chi^{2}$ statistic that
includes contributions for the satellite plane orientations
$\chi^{2}_{{}_{SP,i}}$ (Section 4.1) and counter-rotating fractions
$\chi^{2}_{{}_{CR,i}}$ (Section 4.2) as well as $\chi^{2}_{{}_{TD}}$ for how
the MW disk is heated by the flyby (Section 4.3) and $\chi^{2}_{{}_{HVG}}$
(Section 4.4) for the alignment of the MW-M31 orbital plane with the plane of
HVGs identified by Banik & Zhao (2018b).

$\displaystyle\chi^{2}=\sum_{i=\text{MW,
M31}}\left(\chi^{2}_{{}_{SP,i}}+\chi^{2}_{{}_{CR,i}}\right)~{}+~{}\chi^{2}_{{}_{TD}}~{}+~{}\chi^{2}_{{}_{HVG}}$
(40)

Although the flyby must have affected both the MW and M31 disks, this is much
more difficult to test for M31 than for the MW because M31 is an external
galaxy viewed nearly face-on (Chemin et al., 2009). The accretion history of
M31 is also more complicated and includes recent interactions with other
galaxies, significantly complicating any analysis (Sadoun et al., 2014). For
example, the M31 disk 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 ${R_{{}_{max}}}$ of
the MW as this encompasses the region usually considered in the satellite
plane problem. In addition, we discard particles too close to its disk by
excluding the region with $r_{plane}<r_{{min}}$ and $\left|z\right|<z_{{min}}$
in cylindrical polar co-ordinates $\left(r_{plane},\phi,z\right)$ centred on
the MW with axis along the observed MW disk spin vector. Our adopted values
for these parameters are given in Table 6. We exclude everything within 15
disk scale lengths of its disk plane out to a large distance, completely
removing the MW disk 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 $z_{{}_{min}}$ and almost polar relative to the Galactic disk
(Pawlowski, 2018, figure 1). For M31, we use a similar procedure to the MW,
though with a slightly higher $z_{{}_{min}}$ to compensate for its larger
disk.

Quantity | Value for MW | Value for M31
---|---|---
$R_{{}_{max}}$ | 250 kpc | 250 kpc
$r_{{}_{min}}$ | 172 kpc | 250 kpc
$z_{{}_{min}}$ | 32.25 kpc | 53 kpc
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 $r_{{}_{min}}=R_{{}_{max}}$ for
both galaxies, in practice the shorter scale length of the MW meant a lower
$r_{{}_{min}}$ was sufficient.

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 (Figures 6 and 7).

After selecting ‘satellite’ particles in this way, we create a 2D histogram of
the orbital poles $\widehat{\bm{h}}\equiv\frac{\bm{h}}{\left|\bm{h}\right|}$
of these particles relative to their nearest galaxy, where

$\displaystyle\bm{h}~{}=~{}\left(\bm{r}-\bm{r}_{{}_{i}}\right)\times\left(\bm{v}-\bm{v}_{{}_{i}}\right)$
(41)

Here, the nearest galaxy to the particle is galaxy $i$ with position
$\bm{r}_{{}_{i}}$ and velocity $\bm{v}_{{}_{i}}$. We assign the statistical
weight of the particle to the bin in Galactic latitude and longitude defined
by its orbital pole $\widehat{\bm{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{\bm{h}}_{{}_{peak}}$ in each orbital pole histogram. We then
quantify the angular dispersion of this peak, considering only particles with
$\widehat{\bm{h}}$ in the same hemisphere as the peak i.e. with
$\widehat{\bm{h}}\cdot\widehat{\bm{h}}_{{}_{peak}}>0$. We also quantify the
angular dispersion of the remaining particles’ orbital poles about the
direction $-\widehat{\bm{h}}_{{}_{peak}}$. This allows for the possibility of
a counter-rotating peak i.e. a secondary peak in the orbital pole distribution
at $-\widehat{\bm{h}}_{{}_{peak}}$.

To make our results more robust, we apply a $\sigma$-clipping procedure to
reject outliers from the primary ($\equiv$ co-rotating) and counter-rotating
peaks. The threshold we use is 2.58$\sigma$, corresponding to the 99%
confidence interval of a Gaussian distribution.666If this threshold angle
exceeds ${90^{\circ}}$, we do not perform outlier rejection $-$ this often
happens for the MW. The angular dispersion of each peak is defined as the
weighted rms angle to $\pm\widehat{\bm{h}}_{{}_{peak}}$ of particles whose
orbital pole lies in the appropriate hemisphere. Considering the statistical
weight $W_{j}$ of each particle $j$, the angular dispersion of the counter-
rotating peak is thus

$\displaystyle\sigma_{{}_{-}}~{}=~{}\sqrt{\frac{\sum_{j}W_{j}\left[\overbrace{\cos^{-1}\left(-\widehat{\bm{h}}_{{}_{j}}\cdot\widehat{\bm{h}}_{{}_{peak}}\right)}^{<90^{\circ}}\right]^{2}}{\sum_{j}W_{j}}}$
(42)

After this outlier rejection stage, we obtain the total angular momentum
$\bm{h}_{{}_{tot,\pm}}$ of the remaining particles in each peak, taking care
of their statistical weights and using $+$ to denote the co-rotating side ($-$
for the counter-rotating side). We use this to refine our estimate of
$\widehat{\bm{h}}_{{}_{peak}}$, considering particles corresponding to both
orbital pole peaks.

$\displaystyle\widehat{\bm{h}}_{{}_{peak}}~{}\propto~{}\bm{h}_{{}_{tot,+}}-\bm{h}_{{}_{tot,-}}$
(43)

This leads to a revised measure of each peak’s angular dispersion. Using our
newly estimated ${\widehat{\bm{h}}_{{}_{peak}}}$ and ${\sigma_{{}_{\pm}}}$, we
repeat the $\sigma$-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 ${\,{\sim}\,10}$ steps. In this
way, we obtain the directions of the co-rotating and counter-rotating peaks as
well as their angular dispersions. We use ${\widehat{\bm{h}}_{{}_{peak}}}$ to
increment the $\chi^{2}$ of our model by

$\displaystyle\Delta\chi^{2}~{}$ $\displaystyle=$
$\displaystyle~{}\sum_{i={MW,M31}}\chi^{2}_{{}_{SP,i}}$ (44)
$\displaystyle\chi^{2}_{{}_{SP,i}}~{}$ $\displaystyle=$
$\displaystyle~{}\left[\frac{\cos^{-1}\left(\widehat{\bm{h}}_{{}_{peak}}\cdot\widehat{\bm{h}}_{{}_{SP,i}}\right)}{\sigma_{{}_{SP,i}}}\right]^{2}$
(45)

Satellite plane spin vector of.. | Direction | Error
---|---|---
Milky Way (MW) | $\left(176.4^{\circ},-15.0^{\circ}\right)$ | $25^{\circ}$
Andromeda (M31) | $\left(206.2^{\circ},7.8^{\circ}\right)$ | $15^{\circ}$
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 anti-centre. 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.

Here, $\widehat{\bm{h}}_{{}_{SP,i}}$ is the observed satellite plane spin
vector of galaxy $i$ (Table 7). This table also shows $\sigma_{{}_{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.
$\sigma_{{}_{SP,M31}}<\sigma_{{}_{SP,MW}}$).

###  4.2 Fraction of tidal debris that counter-rotate

The method described in Section 4.1 yields the amount of material in the co-
rotating and counter-rotating peaks, which we call $W_{+}$ and $W_{-}$,
respectively. This can be used to obtain a counter-rotating fraction

$\displaystyle f_{counter}~{}\equiv~{}\frac{W_{-}}{W_{+}~{}+~{}W_{-}}$ (46)

The values of $f_{counter}$ can be compared with the proportion of material in
each satellite plane on a counter-rotating orbit. To penalise models that have
a different $f_{counter}$ to what we expect, we add an appropriate
contribution to the total $\chi^{2}$ of the model.

$\displaystyle\Delta\chi^{2}~{}$ $\displaystyle=$
$\displaystyle~{}\sum_{i={MW,M31}}\chi^{2}_{{}_{CR,i}}$ (47)
$\displaystyle\chi^{2}_{{}_{CR,i}}~{}$ $\displaystyle=$
$\displaystyle~{}\left(\frac{f_{counter,i}-f_{counter,obs,i}}{\sigma_{f_{counter,i}}}\right)^{2}$
(48)

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 2 satellites could well be interlopers or genuine
counter-rotators. As there is no clear evidence for the latter, we assume that
$f_{counter,obs,M31}=0$ but allow an uncertainty of
$\sigma_{f_{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 ${\approx 10}$ Galactic
satellites with known proper motions, we require our models to yield
$f_{counter,obs,MW}=0.10\pm 0.05$.

###  4.3 Milky Way thick disk

####  4.3.1 Age

Our simulations of a past MW-M31 flyby must involve significant tidal torques
on their disks in order to explain their significant misalignments with their
satellite planes ($\,{\sim}\,50^{\circ}$ for M31 and $\,{\sim}\,75^{\circ}$
for the MW). As a result, the MW and M31 disks 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 disk was formed in this
event.

Originally, the Galactic thick disk was identified based on the density of
stars having a double-exponential law with distance out of the Galactic disk
plane (Gilmore & Reid, 1983). The lower mass thick disk component may well
have arisen from dynamical heating of the MW thin disk during its close
encounter with M31 (Zhao et al., 2013). Indeed, the thick disk seems to have
formed fairly rapidly from its thin disk (Hayden et al., 2015). The associated
burst of star formation (Snaith et al., 2014, figure 2) could have further
thickened the MW disk (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\pm 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 $\chi^{2}$ of

$\displaystyle\chi^{2}_{{}_{TD,age}}~{}=~{}\left(\frac{t_{f}-t_{flyby}-9Gyr}{1Gyr}\right)^{2}$
(49)

####  4.3.2 Thickness

In addition to its age, another important property of the thick disk is its
thickness or scale height. As our simulation lacks dissipation, it is most
analogous to that part of the MW stellar disk which existed before the flyby.
The gas disk would likely have settled down after the flyby, leading to a thin
component to the MW stellar disk that would not be accurately included in our
model. Thus, we expect to roughly match the observed scale height of the old
MW thick disk. This was recently measured by Ferguson et al. (2017), who fit
its vertical density profile using

$\displaystyle\rho\left(z\right)~{}\propto~{}\operatorname{sech}^{2}\left(\frac{z}{2H_{{}_{2}}}\right)$
(50)

Based on their figure 9, we take a value of $H_{{}_{2}}=0.7\pm 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 disks.

$\displaystyle z_{{}_{rms,obs}}~{}=~{}\frac{\pi}{\sqrt{3}}H_{{}_{2}}$ (51)

The observed thickness $z_{{}_{rms,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 disk heating arises almost entirely
due to the effect of M31. To account for this, we need to estimate how thick
the Galactic disk 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 $\sigma_{{}_{z}}$. This is
presumably because stars form out of dissipational gas which has settled
itself into a rather thin disk (Nakanishi & Sofue, 2016, figure 6). Older
stars have higher $\sigma_{{}_{z}}$, but this quickly saturates at
${\,{\sim}\,20}$ km/s. We assume that this is what $\sigma_{{}_{z}}$ would
have been without the M31 flyby. However, the thick disk has
$\sigma_{{}_{z}}\approx 40$ km/s (Minchev et al., 2014, figure 1). This
suggests that internal processes were responsible for
${\,{\sim}\,\frac{1}{4}}$ of the disk heating required to explain the Galactic
thick disk, 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 disk which is thinner than observed. To account for this,
we only require our models to reach

$\displaystyle z_{{}_{TD}}~{}=~{}\frac{3}{4}z_{{}_{rms,obs}}$ (52)

This can be compared with the simulated value of $z_{{}_{rms}}$. Following
Ferguson et al. (2017, figure 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 maximises the number
of particles used to determine the $z_{{}_{rms}}$ of our simulated thick disk.

Vertically, we allow particles up to 21.5 kpc either side of the MW disk,
ensuring this is properly included. We then find $z_{{}_{rms}}$ using a plane
fitting procedure similar to that described in section 5.1 of Banik & Zhao
(2018b) augmented by $\sigma$-clipping with a threshold of $2.58\sigma$,
corresponding to the 99% confidence interval of a Gaussian distribution. We
use as many $\sigma$-clipping stages as we need for the algorithm to converge,
without reintroducing particles excluded at a previous stage. Usually, this
process requires ${\,{\sim}\,10}$ iterations. Once it converges, we can
determine $\chi^{2}_{{}_{TD,height}}$ based on the appropriately scaled
uncertainty. This allows us to add the thick disk contribution to the total
$\chi^{2}$ of each model.

$\displaystyle\chi^{2}_{{}_{TD}}~{}=~{}\chi^{2}_{{}_{TD,age}}+\chi^{2}_{{}_{TD,height}}$
(53)

###  4.4 Alignment with the high-velocity galaxy 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, Figure
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 $\Lambda$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 a 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.

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 $\Lambda$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
$\sim\cos^{-1}0.9\approx 25^{\circ}$ (Banik & Zhao, 2018b, figure 7).777We use
$\sim$ for order of magnitude estimates and $\approx$ when making more precise
statements. An uncertainty of this order seems reasonable given the
$16^{\circ}$ 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 $\sigma_{{}_{HVG}}=25^{\circ}$
between the HVG plane normal $\widehat{\bm{n}}_{{}_{HVG}}$ and the MW-M31
orbital pole $\widehat{\bm{h}}_{{}_{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 $\chi^{2}$ of

$\displaystyle\chi^{2}_{{}_{HVG}}~{}=~{}\left(\frac{\cos^{-1}\left|\widehat{\bm{n}}_{{}_{HVG}}\cdot\widehat{\bm{h}}_{{}_{LG}}\right|}{\sigma_{{}_{HVG}}}\right)^{2}$
(54)

##  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 ${\chi^{2}}$
statistic for each model (Equation 40) 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 disks 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.

Quantity | Observed | $\sigma$ | Model | $\chi^{2}$
---|---|---|---|---
Time since flyby, Gyr | 9 | 1 | 7.65 | 1.82
Milky Way (MW) | $\begin{bmatrix}176.4^{\circ}\\\ -15.0^{\circ}\end{bmatrix}$ | $25^{\circ}$ | $\begin{bmatrix}180.6^{\circ}\\\ -42.8^{\circ}\end{bmatrix}$ | 1.26
satellite plane |  |
Andromeda (M31) | $\begin{bmatrix}206.2^{\circ}\\\ 7.8^{\circ}\end{bmatrix}$ | $15^{\circ}$ | $\begin{bmatrix}212.4^{\circ}\\\ 0.14^{\circ}\end{bmatrix}$ | 0.43
satellite plane |  |
MW-M31 orbital pole | $\begin{bmatrix}26.6^{\circ}\\\ 31.8^{\circ}\end{bmatrix}$ | $25^{\circ}$ | $\begin{bmatrix}24.6^{\circ}\\\ -16.3^{\circ}\end{bmatrix}$ | 3.70
(expect HVG plane) |  |
Counter-rotating | 0 | 0.05 | 0 | 0
fraction for M31 |  |  |  |
Counter-rotating | 0.1 | 0.05 | 0.22 | 5.84
fraction for MW |  |  |  |
MW thick disk $z_{{}_{rms}}$ | 0.95 | 0.20 | 0.75 | 1.03
at Solar Circle, kpc |  |  |  |
Table 8: Comparison of our best-fitting model with observational constraints
on the LG. The uncertainties $\sigma$ are shown in the central column e.g.
$\sigma_{{}_{HVG}}=25^{\circ}$. 0.0045 of the mass in the MW ends up in its
satellite plane, including both the co-rotating 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\times 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).

The total $\chi^{2}$ of our best model is 14.08. There are three model
parameters and 7 contributions to $\chi^{2}$. However, the orientation of a
plane is a 2D concept such that the satellite plane and HVG constraints should
really count for 6 constraints rather than 3. In this case, there are 10
constraints $-$ 3 model parameters $=$ 7 ‘degrees of freedom’, for which the
chance of an even higher $\chi^{2}$ is $P=4.98\%$ based on the chi2cdf
function of matlab${}^{\text{\textregistered}}$. 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\%$. 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 $\chi^{2}$ analysis neglecting one
contribution to $\chi^{2}$ at a time. The best-fitting model is unaffected by
the constraints from the thick disk age, the thickness of the MW disk, 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^{\circ}$ further south and
the model prefers a weaker present $g_{{}_{ext}}$ of 0.02$a_{{}_{0}}$. 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^{\circ}$ 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

Figure 4 shows the MW-M31 trajectory in our best-fitting model, with the main
model parameters summarised in Table 9. The flyby occurs 7.65 Gyr ago,
consistent with the thick disk age of ${9\pm 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\times 10^{11}M_{\odot}$. This is
slightly higher than our estimate of the total baryonic mass in their disks,
assuming the MW rotation curve flatlines at ${v_{{}_{f}}=180}$ km/s (Kafle et
al., 2012) while that of M31 flattens at 225 km/s (Carignan et al., 2006).
Taken at face value, these rotation speeds imply a total mass of $2.3\times
10^{11}M_{\odot}$ because the MOND dynamical mass of a galaxy is

$\displaystyle M_{dyn}~{}~{}=~{}~{}\frac{{v_{{}_{f}}}^{4}}{Ga_{{}_{0}}}$ (55)

Quantity | Value
---|---
Orbital pole of MW-M31 orbit | $\left(24.6^{\circ},-16.3^{\circ}\right)$
MW-M31 separation | 783 kpc
MW-M31 tangential velocity | 16.70 km/s
MW$+$M31 mass | $2.81\times 10^{11}M_{\odot}$
Flyby time (after Big Bang) | 6.17 Gyr
Closest approach distance | 19.6 kpc
Present external field on LG | $0.022~{}a_{{}_{0}}$
External field direction | $\left(276^{\circ},30^{\circ}\right)$
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.029a_{{}_{0}}}$ 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 figure 3). The grid of models we run has a step of
$0.001a_{{}_{0}}$ in the present EF on the LG, 1.15 km/s in the tangential
velocity of M31 (${\approx 1.7}$ kpc in closest approach distance) and
${5^{\circ}}$ in the MW-M31 orbital pole.

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 ${v_{{}_{f}}\approx 70}$ km/s (Alves & Nelson, 2000;
van der Marel & Kallivayalil, 2014) while ${v_{{}_{f}}\approx 120}$ km/s for
M33 (Corbelli, 2003; Corbelli & Salucci, 2007). Thus, the most luminous LG
satellites only contribute a mass of ${\approx 1.5\times 10^{10}M_{\odot}}$
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% 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).

Figure 4: _Top_ : The MW-M31 separation ${d\left(t\right)}$ 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. 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{\bm{h}}_{{}_{LG}}$ is not a
good guide to what $\widehat{\bm{h}}_{{}_{LG}}$ was several Gyr ago.

The high timing argument mass might also be a sign of hot gas halos
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\times 10^{10}M_{\odot}$ (see their
table 2 model A). This is consistent with the estimate of $2.7\pm 1.4\times
10^{10}M_{\odot}$ obtained by Salem et al. (2015) based on how the LMC gas
disc is truncated at a smaller radius than its stellar disk, 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, figure 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 et al., 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 ${\,{\sim}\,6}$ Gyr ago. This is much more recent than the observed
${9\pm 1}$ Gyr age of the Galactic thick disk 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 35.
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 figure 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{\bm{h}}_{{}_{LG}}$ is changed while not adjusting
$\bm{g}_{{}_{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.

Perturber | Estimated $g_{{}_{tide}}$ at | Estimated
---|---|---
| first apocentre, m/s2 |  $g_{{}_{tide}}$ now, m/s2
Cen A | $-1.85\times 10^{-13}$ | $+6.22\times 10^{-14}$
M81 | $-6.80\times 10^{-14}$ | $-4.11\times 10^{-15}$
IC 342 | $-5.33\times 10^{-13}$ | $+9.75\times 10^{-15}$
GA | $-8.34\times 10^{-15}$ | $+1.64\times 10^{-14}$
$\frac{\ddot{a}}{a}d_{{}_{M31}}$ | $-6.47\times 10^{-13}$ | $+6.06\times 10^{-14}$
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 Great Attractor
(GA), we use Equation 34 instead. The last row shows the contribution of the
cosmological acceleration term. The MW and M31 are expected to have a relative
acceleration of ${\approx 2\times 10^{-12}}$ m/s2 once the external field is
considered (Section 5.1.2).

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_{{}_{LG}}a_{{}_{0}}}}{d_{{}_{M31}}}=2.77\times 10^{-12}$ m/s2,
though this would be weakened ${\approx 30\%}$ due to the EF (Figure 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.

However, cosmic expansion can’t be the entire explanation because the tidal
forcing changes sign for Cen A and IC 342, something that expansion alone
can’t 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{\bm{d}}_{{}_{M31}}$ rotate
${233^{\circ}}$ 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 disk.

####  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 over-simplifying 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
${\,{\sim}\,120}$ Myr (Figure 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 55 suggests that
the LMC has 2.25% as much mass as the MW (Alves & Nelson, 2000; Kafle et al.,
2012) while M33 has 7.5% as much mass as M31 (Corbelli, 2003; Carignan et al.,
2006). Combining this information yields the relative velocity between the MW-
LMC barycentre and the M31-M33 barycentre. We fix the radial component of this
velocity at -93.4 km/s but let the tangential component vary so as to alter
the MW-M31 orbital pole $\widehat{\bm{h}}_{{}_{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 figure 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{\bm{h}}_{{}_{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 ${\approx 200}$ km/s 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, even if we only consider the more massive M33. If we
assume uncertainties of 10% on its rotation curve, this implies a 6 km/s
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 $\approx 200$ kpc (Patel
et al., 2017, table 1), the reflex motion of M31 must be ${\approx 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{\bm{h}}_{{}_{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 backwards 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 ${\,{\sim}\,75^{\circ}}$ of Galactic co-
ordinates $\left(227^{\circ},-34^{\circ}\right)$. 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, figure 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 (Figure 5). Both
orbital pole distributions show a lack of material on a similar orbit to the
disk 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: The distribution of orbital angular momentum directions (spin
vectors) for tidal debris around the MW and M31 disks 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 disks
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 disk 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 disk (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^{\circ}}$, similar to the observed spin vector of
the MW satellite plane. The Galactic latitude of the simulated satellite plane
is ${\,{\sim}\,25^{\circ}}$ 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 disk 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
${\approx 25^{\circ}}$ (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^{\circ}}$ for the MW satellite plane but only ${23^{\circ}}$ 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^{\circ}}$ (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^{\circ}}$ (Figure 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% 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 onto 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 two 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 co-rotating, 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 halos 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 co-rotating 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
${v_{{}_{f}}=70}$ km/s (Alves & Nelson, 2000) while the MW has
${v_{{}_{f}}=180}$ km/s (Kafle et al., 2012), then we expect the satellite
plane mass fraction to be $\,{\sim}\,\left(\frac{70}{180}\right)^{4}=0.02$
(Equation 55).

For M31, this quantity is only ${6.5\times 10^{-6}}$ in our model.
Observationally, M31 satellite plane galaxies are each ${\approx 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, figure 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 et al., 2014). Although that work considered
$\Lambda$CDM, this is also possible in MOND because M32 would experience some
dynamical friction against the baryonic disk 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 ${\approx 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 ${\approx
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

Figure 6: 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 disk (Section 3). The
particles are divided according to whether they are currently in the co-
rotating 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 et al., 2016), but most of
this structure is defined by less distant satellites (Pawlowski & Kroupa,
2013, figure 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, figure 1).

Our low predicted satellite plane masses may be related to our assumed
exponential disk 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 Figures 6 and 7, whose left 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 ${\approx 50}$ kpc from
the MW, even though its stellar disk has a scale length of only 2.15 kpc.
Thus, our simulated satellite plane masses are sensitive to assumptions
concerning the MW disk surface density at large radii. Partly to deal with
this, we also include a gas disk 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.

For M31, we use a similar initial radial range $-$ our simulations start with
M31 particles out to 106 kpc from M31 (20 disk scale lengths, Courteau et al.,
2011). This seems to be enough for our purposes (Figure 7). However, this is
such a large distance that our estimated M31 satellite plane mass is rather
sensitive to its disk surface density at large radii. Even if M31 can be
treated as an exponential disk, there is a ${10\%}$ 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 (Figure 7). This fraction would
be ${2.4\times}$ higher if we assume the M31 disk 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 Figure
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 (Figure 5). The
particles in the secondary peak start at similar radii to the co-rotating
particles (Figure 6). This result is different to that obtained by Pawlowski
et al. (2011), whose figure 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 7: Similar to Figure 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 (Figure 5).

####  5.2.5 Present radial distribution

The right sides of the panels in Figure 6 (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, figure
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 co-rotating particles. This is similar to the result obtained by
Pawlowski et al. (2011), as can be seen from their figure 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 co-rotating 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.

Figure 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 (bottom panel of Figure 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, figure 1).888At the 783 kpc distance of M31
(McConnachie, 2012), each square in this figure has a side length of 55 kpc.
Indeed, our simulated M31 satellite particle radial distribution peaks at
${\approx 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 59). As a
proxy, we show the cosine distribution of angles between the position and
velocity of each particle with respect to its host galaxy (Figure 8), with
weights assigned as in Section 3.

Figure 8: _Left_ : Histogram of $\widehat{\bm{r}}\cdot\widehat{\bm{v}}$, where
$\bm{r}$ and $\bm{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 co-rotating peak around M31 as there are
no counter-rotating particles (Figure 5).

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^{-1}0.9=26^{\circ}$ of the radial direction.
This is unlike the actual MW satellite plane as that is co-rotating 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{\bm{r}}\cdot\widehat{\bm{v}}$ suggests
that the MW satellite system has reached some sort of equilibrium in the
${\approx 8}$ Gyr since its interaction with M31.

We obtain rather different results for the satellite particles around M31
(right panel of Figure 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 Milky Way thick disk

Figure 9: Weighted histogram of the vertical distances of MW particles from
its disk. Here, we only show particles with cylindrical radius between 7.2 kpc
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).

An important constraint on our models is the extent to which they disrupt the
well-observed MW disk 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 $r_{plane}$ within 1 kpc of the
Solar value (Ferguson et al., 2017, figure 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 (Figure 9). Their distance from the
MW disk 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).

The lower scale height than observed may be due to the Galactic disk 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 52), the
MW disk 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 disk
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 disk. We thus use Figure 10 to show the rms height
of our simulated MW disk for several bins in $r_{plane}$, using a similar
method to that used in Section 4.3.2. This reveals that the thick disk 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.

However, the trend may well be opposite this if the MW disk 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 disk originated in a tidal heating event. If it did not, this would
cast doubt on our scenario.

Interestingly, the thick disk does appear to have a higher velocity dispersion
further from the MW (Banik, 2014). The apparent lack of accreted stars