pdftexdestiPropernation with the same hyperrefOption hyperrefToken
pdftex(dest)

# 3D hydrodynamic simulations for the formation of the Local Group satellite
planes

Indranil Banik1,2, Ingo Thies2, Roy Truelove1, Graeme Candlish3, Benoit
Famaey4, Marcel S. Pawlowski5, Rodrigo Ibata4 & Pavel Kroupa2,6
1Scottish Universities Physics Alliance, University of Saint Andrews, North
Haugh, Saint Andrews, Fife, KY16 9SS, UK
2Helmholtz-Institut für Strahlen und Kernphysik (HISKP), University of Bonn,
Nussallee 14$-$16, D-53115 Bonn, Germany
3Instituto de Física y Astronomía, Universidad de Valparaíso, Gran Bretaña
1111, Valparaíso, Chile
4Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de
Strasbourg, 11 rue de l’Université, 67000 Strasbourg, France
5Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, D-14482
Potsdam, Germany
6Astronomical Institute, Faculty of Mathematics and Physics, Charles
University, V Holešovičkách 2, CZ-180 00 Praha 8, Czech Republic Email:
<EMAIL_ADDRESS>(Indranil Banik)
<EMAIL_ADDRESS>(Ingo Thies)

###### Abstract

The existence of mutually correlated thin and rotating planes of satellite
galaxies around both the Milky Way (MW) and Andromeda (M31) calls for an
explanation. Previous work in Milgromian dynamics (MOND) indicated that a past
MW-M31 encounter might have led to the formation of these satellite planes. We
perform the first-ever hydrodynamical MOND simulation of the Local Group using
phantom of ramses. We show that an MW-M31 encounter at $z\approx 1$, with a
perigalactic distance of about 80 kpc, can yield two disc galaxies at $z=0$
oriented similarly to the observed galactic discs and separated similarly to
the observed M31 distance. Importantly, the tidal debris are distributed in
phase space similarly to the observed MW and M31 satellite planes, with the
correct preferred orbital pole for both. The MW-M31 orbital geometry is
consistent with the presently observed M31 proper motion despite this not
being considered as a constraint when exploring the parameter space. The mass
of the tidal debris around the MW and M31 at $z=0$ compare well with the mass
observed in their satellite systems. The remnant discs of the two galaxies
have realistic radial scale lengths and velocity dispersions, and the
simulation naturally produces a much hotter stellar disc in M31 than in the
MW. However, reconciling this scenario with the ages of stellar populations in
satellite galaxies would require that a higher fraction of stars previously
formed in the outskirts of the progenitors ended up within the tidal debris,
or that the MW-M31 interaction occurred at $z>1$.

###### keywords:

gravitation – Local Group – galaxies: formation – galaxies: interactions –
galaxies: kinematics and dynamics – hydrodynamics

††pubyear: 2022††pagerange: 3D hydrodynamic simulations for the formation of
the Local Group satellite planes–E

##  1 Introduction

It has been known since the early work of Kunkel & Demers (1976), Lynden-Bell
(1976), and Lynden-Bell (1982) that dwarf spheroidal galaxies around the Milky
Way (MW) have an anisotropic spatial distribution, as later confirmed by
Kroupa, Theis & Boily (2005) and Metz, Kroupa & Jerjen (2007). The orbital
poles that can be measured indicate that this Vast Polar Structure (VPOS) is
corotating (Metz et al., 2008; Pawlowski et al., 2013; Pawlowski et al., 2017;
Pawlowski & Kroupa, 2020; Li et al., 2021). This is in stark contrast with
expectations based on the standard $\Lambda$ cold dark matter ($\Lambda$CDM)
cosmological model (Efstathiou, Sutherland & Maddox, 1990; Ostriker &
Steinhardt, 1995) because the existence of such a corotating plane would a
priori imply a significant amount of dissipation. The issue is compounded by a
similar satellite plane (SP) discovered around Andromeda (M31; Ibata et al.,
2013; Ibata et al., 2014b) and Centaurus A (Cen A; Müller et al., 2018a;
Müller et al., 2021), with radial velocities (RVs) suggestive of corotation in
both cases.111Hints of the M31 SP were already evident in Metz et al. (2007)
and Metz, Kroupa & Jerjen (2009a). Hints of the Cen A SP were evident in Tully
et al. (2015), but the two planes they identified were later revealed to be
part of one thicker plane (Müller et al., 2016). Proper motions (PMs) have
recently become known for two M31 SP members, and indeed indicate that it too
is likely corotating within its plane (Sohn et al., 2020).

Such structures are difficult to explain using dissipationless haloes of CDM,
as first pointed out by Kroupa et al. (2005) and more recently by Pawlowski et
al. (2014), who considered and excluded a wide range of different proposed
explanations within the $\Lambda$CDM context. For instance, accreting most
satellites as a single group (Metz et al., 2009b) or along a single filament
would yield some anisotropy, but not enough to explain the very thin MW SP
(Shao et al., 2018). The two available PMs of M31 SP members are also in
tension with $\Lambda$CDM expectations because they indicate motion nearly
within the M31 SP (Pawlowski & Tony Sohn, 2021). Dissipationless collapse of
CDM halos is therefore insufficient to account for the observed SPs if they
are composed of primordial dwarfs. Including baryonic physics does not change
this picture very much (Ahmed et al., 2017; Pawlowski & Kroupa, 2020; Samuel
et al., 2021).

Since dissipation in the CDM component would cause a significant Galactic DM
disc that is in tension with observations (Buch et al., 2019; Widmark et al.,
2021), an obvious possibility is that the necessary dissipation occurred in
baryons. Given that the MW SP is almost polar with respect to its disc, this
would require a tidal interaction with another galaxy (Pawlowski, Kroupa & de
Boer, 2011, 2012). In this scenario, the MW SP would consist of tidal dwarf
galaxies (TDGs) that condensed out of gas-rich tidal debris expelled during a
past interaction. A similar scenario could have occurred around M31. The
formation of TDGs has been observed outside the Local Group (LG), for instance
around the Antennae (Mirabel, Dottori & Lutz, 1992) and the Seashell Galaxy
(Bournaud et al., 2004). The LG SPs may have formed analogously. A common
origin for both LG SPs is possible as a result of a major merger experienced
by M31 (Hammer et al., 2010) expelling tidal debris towards the MW (Hammer et
al., 2013), perhaps explaining the high Galactocentric velocities of MW
satellites (Hammer et al., 2021). Since the M31 SP is viewed close to edge-on
from our vantage point in the MW (Conn et al., 2013; Ibata et al., 2013;
Santos-Santos et al., 2020), it is possible that if M31 experienced a major
merger, then some of the tidal debris expelled from M31 formed its SP while
some reached a much larger distance and is now close to the MW.

In a $\Lambda$CDM context, TDGs would be free of DM due to its dissipationless
nature and its initial distribution in a dispersion-supported near-spherical
halo (Barnes & Hernquist, 1992; Wetzstein, Naab & Burkert, 2007). During a
tidal interaction, DM of this form is clearly incapable of forming into a thin
dense tidal tail which might undergo Jeans collapse into TDGs. As a result,
any TDGs would be purely baryonic and thus have only a very small escape
velocity, preventing them from subsequently accreting DM out of the Galactic
halo. For this reason, TDGs formed in the $\Lambda$CDM framework cannot
explain the high observed internal velocity dispersions
($\sigma_{\text{int}}$) of the MW satellites if they are in equilibrium. Non-
equilibrium solutions were considered by Kroupa (1997), Klessen & Kroupa
(1998), and Casas et al. (2012). While their proposed solutions match many
properties of the observed MW satellites, this scenario would not yield
elevated $\sigma_{\text{int}}$ values at Galactocentric distances $\ga 150$
kpc. Satellites in such a non-equilibrium phase would also be very fragile,
requiring us to be observing them at just the right epoch prior to total
disruption but after significant tidal perturbation. These issues are
undoubtedly a major reason why CDM-free galaxies are very rare in the latest
$\Lambda$CDM simulations (Haslbauer et al., 2019), making it very unlikely
that so many TDGs are in the LG right now around both the MW and M31.
Moreover, postulating that most of their satellites are TDGs would mean that
they have very few primordial satellites.

This conundrum presents an open invitation to consider a different theoretical
framework in which DM is dynamically irrelevant to holding galaxies together
against their high $\sigma_{\text{int}}$, removing the problem that TDGs are
expected to be free of DM. In this case, the effects conventionally attributed
to CDM must instead be due to a non-Newtonian gravity law on galactic scales.
The best developed such proposal is Milgromian dynamics (MOND; Milgrom, 1983).
In MOND, galaxies lack CDM but the gravitational field strength $g$ at
distance $d$ from an isolated point mass $M$ transitions from the Newtonian
$g_{{}_{N}}={GM/d^{2}}$ law at short range to

$\displaystyle g~{}=~{}\sqrt{g_{{}_{N}}a_{{}_{0}}}~{}~{}~{}\text{for
}~{}g_{{}_{N}}\ll a_{{}_{0}}\,.$ (1)

MOND introduces $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 (RCs; Begeman, Broeils & Sanders, 1991; Gentile, Famaey & de
Blok, 2011). With this value of $a_{{}_{0}}$, MOND predicts the detailed shape
of galaxy RCs very well using only their directly observed baryonic matter
(e.g. Kroupa et al., 2018a; Li et al., 2018; Sanders, 2019). In particular,
observations confirm the prior MOND prediction of very large departures from
Newtonian dynamics in low surface brightness galaxies (LSBs; McGaugh, 2021,
and references therein). More generally, there is a very tight empirical
‘radial acceleration relation’ (RAR) between the gravity inferred from RCs and
that expected from the baryons alone in Newtonian dynamics (McGaugh et al.,
2016; Lelli et al., 2017), with RCs asymptotically reaching a flatline level
$\propto\sqrt[{}^{4}]{M}$ as required by Equation 1 (McGaugh, 2012; Lelli et
al., 2019). The observed phenomenology of galactic RCs confirm all the central
predictions of Milgrom (1983), as reviewed in e.g. Famaey & McGaugh (2012).

MOND can also explain the X-ray temperature profile (Milgrom, 2012) and
internal dynamics of elliptical galaxies (figure 8 of Lelli et al., 2017),
which reveal a similar characteristic acceleration scale to spirals (Chae et
al., 2020a; Shelest & Lelli, 2020). At the low-mass end, MOND is consistent
with the $\sigma_{\text{int}}$ of pressure-supported systems like the
satellites of the MW (McGaugh & Wolf, 2010), M31 (McGaugh & Milgrom, 2013a,
b), non-satellite LG dwarfs (McGaugh et al., 2021), Dragonfly 2 (DF2; Famaey
et al., 2018; Kroupa et al., 2018b), and DF4 (Haghi et al., 2019a). These
predictions rely on correctly including the external field effect (EFE)
arising from the non-linearity of MOND (Milgrom, 1986), which we discuss
further in Section 2.4.3. For the more isolated galaxy DF44, the MOND
prediction without the EFE is consistent with the observed
$\sigma_{\text{int}}$ profile (Bílek et al., 2019; Haghi et al., 2019b). Note
that the situation with the EFE is however less clear for ultra-diffuse
galaxies located deep in the potential well of galaxy clusters (Freundlich et
al., 2022). The EFE also plays a role in accounting for the observed weak bar
of M33 (Banik et al., 2020), which is quite difficult to understand in the
presence of a live DM halo (Sellwood, Shen & Li, 2019) due to bar-halo angular
momentum exchange (e.g. Debattista & Sellwood, 2000; Athanassoula, 2002). This
problem is related to the bar pattern speeds in galaxies, which seem to be too
slow in $\Lambda$CDM cosmological simulations due to dynamical friction on the
bar exerted by the DM halo (Algorry et al. 2017; Peschken & Łokas 2019; Roshan
et al. 2021a, b; though see Fragkoudi et al. 2021). In addition, the
morphological properties of MW satellites are more easily understood in MOND
as a consequence of differing levels of tidal stability (McGaugh & Wolf,
2010). This is because in MOND, the lack of CDM haloes and the EFE render the
satellites much more susceptible to Galactic tides, which are not so relevant
in $\Lambda$CDM (see their figure 6). Neglecting tides leads to erroneous
conclusions regarding the viability of MOND since it predicts that much fewer
satellites $-$ especially at the ultra-faint end $-$ are amenable to
equilibrium virial analysis (Fattahi et al., 2018).

Further work will however be necessary to make rigorous MOND predictions in a
cosmological context, which needs a relativistic theory. A few such theories
exist, with one promising proof-of-concept being the relativistic MOND theory
of Skordis & Złośnik (2019) in which gravitational waves travel at the speed
of light. Such theories allow MOND calculations of weak gravitational lensing
by foreground galaxies in stacked analyses (e.g. Brimioulle et al., 2013),
which so far seems to agree with expectations (Milgrom, 2013; Brouwer et al.,
2017, 2021). At larger distances from the central galaxy, the EFE from
surrounding structures would cause the gravity law to become inverse square
and to depart from spherical symmetry (Banik & Zhao, 2018a), which may explain
some recent observations (Schrabback et al., 2021). More detailed calculations
would require knowledge of how large-scale structure forms in a Milgromian
framework and the resulting EFE on galaxies. This is a critical next step for
MOND, though care is required when comparing with observations as these could
have a rather different interpretation to what is usually assumed. One
possible smoking gun signature of MOND in weak lensing convergence maps would
be the discovery that the convergence parameter is negative in some regions
(Oria et al., 2021). This cannot arise in GR and is possible only in
gravitational theories that are non-linear in the weak-field regime.

At this stage, general statements can however already be made about a MONDian
cosmology, in which the key difference with $\Lambda$CDM would be faster
structure formation (Sanders, 1998). This might be relevant for the so-called
Hubble tension, i.e. the fact that, to fit the anisotropies in the cosmic
microwave background (CMB; Aiola et al., 2020; Planck Collaboration VI, 2020),
$\Lambda$CDM requires a local Hubble constant $H_{0}$ below the directly
measured value at high significance based on multiple independent techniques
(e.g. Riess, 2020; Di Valentino, 2021; Riess et al., 2021b). This tension
could be due to our position within a large local supervoid underdense by
$\approx 30\%$ out to a radius of $\approx 300$ Mpc. Such a large and deep
underdensity has indeed been observed in multiple surveys and is called the
KBC void after its discoverers (Keenan, Barger & Cowie, 2013). This is
incompatible with $\Lambda$CDM at $6.04\sigma$, one major reason for the high
significance of the tension being that the relevant observations cover 90% of
the sky (Haslbauer, Banik & Kroupa, 2020). However, a KBC-like void could
arise naturally in a MOND cosmology supplemented by light sterile neutrinos
playing the role of hot DM (HDM), as proposed by Angus (2009). The main
feature in the Haslbauer et al. (2020) void scenario is faster structure
formation than in $\Lambda$CDM. From an observational point of view, this is
also evident in the properties of the high-redshift interacting galaxy cluster
El Gordo (Asencio, Banik & Kroupa, 2021). The lack of analogues to El Gordo at
low redshift could well be due to our location within the KBC void, which
might also explain why the MOND simulations of Angus et al. (2013) seemingly
overproduced massive clusters when comparing their whole simulation volume
with low-redshift datasets. Therefore, MOND with HDM could potentially account
for astronomical observations ranging from the kpc scales of galaxies (where
HDM would play no role; Angus, 2010) all the way to the Gpc scale of the local
supervoid, without causing any obvious problems in the early Universe or in
galaxy clusters (Angus, Famaey & Diaferio, 2010). It is also possible to fit
the CMB in MOND without any sterile neutrino component (Skordis & Złośnik,
2021), though it is unclear whether this approach can explain the properties
of galaxy clusters. For a recent review of MOND that considers evidence from a
wide range of scales and discusses the cosmological aspects in some detail, we
refer the reader to Banik & Zhao (2022).

Concerning the history of the well-observed LG, MOND implies a very strong
mutual attraction between the MW and M31. 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), consistent with the timing of a few other events
putatively linked to the flyby (Section 4.2). An $N$-body simulation of this
interaction showed that it is likely to yield anisotropically distributed
tidal debris reminiscent of an SP (Bílek et al., 2018). Around the same time,
Banik, O’Ryan & Zhao (2018, hereafter BRZ18) considered a wider range of
orbital geometries using a less computationally intensive restricted $N$-body
approach in which the MW and M31 were treated as point masses surrounded by
test particle discs.222This leads to a numerically more tractable axisymmetric
potential. Their section 2 demonstrated that the MW-M31 trajectory is
consistent with negligible peculiar velocity in the early universe, a
constraint known as the timing argument (Kahn & Woltjer, 1959). Despite
lacking hydrodynamics or disc self-gravity, the initial setup of each galaxy
as a rotating disc was sufficient to cause significant clustering of the tidal
debri orbital poles. In some models, the preferred directions aligned with the
actually observed orbital poles of the LG SPs. Therefore, these structures
could well have formed as TDGs that condensed out of tidal debris orbiting in
the correct plane. Indeed, Bílek et al. (2021) conclude that all satellite
galaxies in the LG SPs are TDGs based on the earlier simulations of Bílek et
al. (2018).

TDGs are expected to be more resilient in MOND due to their enhanced self-
gravity, as explored with earlier high-resolution hydrodynamical MOND
simulations (Renaud, Famaey & Kroupa, 2016). Besides helping them to survive,
this would also explain the high observed $\sigma_{\text{int}}$ of satellite
galaxies around the MW (McGaugh & Wolf, 2010) and M31 (McGaugh & Milgrom,
2013a, b). However, the enhancement to the self-gravity would typically be
less than provided by a CDM halo, leading to a greater degree of tidal
susceptibility in MOND. This is more in line with the observed morphologies of
LG satellites (McGaugh & Wolf, 2010).

A key distinguishing characteristic of more recently formed TDGs is their high
metallicity for their mass (Duc et al., 2014). However, this relies on a long
process of metal enrichment in the disc of the progenitor galaxy. When
considering a very ancient interaction, there might not have been enough time
for such an enrichment (Recchi, Kroupa & Ploeckinger, 2015), especially as
TDGs are expected to form out of material that was initially several disc
scale lengths out (BRZ18). This could explain why the M31 satellites within
and outside its SP have similar properties, including in terms of metallicity
(Collins et al., 2015).

A past MW-M31 flyby also has interesting consequences for the rest of the LG.
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. As discussed further in Section 4.2.2, this could lead to the existence
of LG dwarfs with an unusually high RV in a $\Lambda$CDM context, such as the
dwarfs in the NGC 3109 association (Pawlowski & McGaugh, 2014; Peebles, 2017).
These could be backsplash from the MW-M31 flyby, a scenario that was
considered in detail by Banik & Zhao (2018c). Backsplash galaxies also exist
in $\Lambda$CDM, but very rarely have properties resembling NGC 3109 (Banik et
al., 2021).

In this contribution, we build on the earlier studies of Bílek et al. (2018)
and BRZ18 by conducting 3D hydrodynamical MOND simulations of the flyby using
phantom of ramses (por; Lüghausen et al., 2015; Nagesh et al., 2021). Our main
objective is to vary the MW-M31 orbital pole and pericentre distance to find
if there are models where the tidal debris around each galaxy aligns with its
observed SP. Achieving this simultaneously for both the MW and M31 is a highly
non-trivial test of the past flyby scenario. We also check if their discs are
preserved and end up with realistic properties.

In the following, we describe the initial conditions and setup of our
simulations (Section 2). We then present our results and analyses regarding
the MW-M31 trajectory (Section 3.1) and proper motion (PM; Section 3.2), the
tidal debris (Section 3.3), and the MW and M31 disc remnants (Section 3.4). We
discuss our results in Section 4 and conclude in Section 5. Videos of the LG
in our best-fitting simulation with frames every 10 Myr are publicly
available.333https://seafile.unistra.fr/d/6bb8e94212764324868e/

##  2 Methods

###  2.1 Poisson equation

The simulations presented in this paper are conducted with por, which solves
the governing equation of QUMOND:

$\displaystyle\overbrace{\nabla\cdot\bm{g}}^{\propto\rho_{\text{eff}}\,\equiv\,\rho_{{}_{\text{PDM}}}+\rho_{{}_{b}}}~{}=~{}\nabla\cdot\left[\bm{g}_{{}_{N}}\nu\left(g_{{}_{N}}\right)\right]\,,$
(2)

where $\bm{g}_{{}_{N}}$ is the Newtonian gravity determined from the baryonic
density $\rho_{{}_{b}}$ using standard techniques,
$v\equiv\left|\bm{v}\right|$ for any vector $\bm{v}$, and $\bm{g}$ is the true
gravity. It is often helpful to think of what density distribution would lead
to this $\bm{g}$ under Newtonian gravity. The required effective density
$\rho_{{}_{\text{eff}}}\equiv\rho_{{}_{\text{PDM}}}+\rho_{{}_{b}}$, where
$\rho_{{}_{\text{PDM}}}$ is the phantom dark matter (PDM) density which
captures the MOND corrections. Equation 2 is derived from a non-relativistic
Lagrangian (Milgrom, 2010), so QUMOND obeys the usual conservation laws
regarding energy and momentum. This is also true of an earlier version
(Bekenstein & Milgrom, 1984), though we do not consider it here as it is
computationally less efficient due to a non-linear grid relaxation stage.

To solve Equation 2, we must assume an interpolating function $\nu$ between
the Newtonian and Milgromian regimes. In spherical symmetry, this has the
effect that $\bm{g}=\nu\bm{g}_{{}_{N}}$, softening the transition between the
Newtonian inverse square law and Equation 1. In this work, we use the ‘simple’
form of the interpolating function (Famaey & Binney, 2005):

$\displaystyle\nu\left(g_{{}_{N}}\right)~{}=~{}\frac{1}{2}~{}+~{}\sqrt{\frac{1}{4}+\frac{a_{{}_{0}}}{g_{{}_{N}}}}\,.$
(3)

This is numerically rather similar to the function used by McGaugh et al.
(2016) and Lelli et al. (2017) to fit galaxy RCs, but can be inverted
analytically. Other reasons for using this function were discussed in section
7.1 of Banik & Zhao (2018d) in preference to functions with a sharper
transition. Equation 3 is quite accurate for
$g_{{}_{N}}\approx\left(0.1-10\right)a_{{}_{0}}$ as relevant for the MW-M31
flyby problem, but Solar system constraints imply that a more rapid
convergence to the Newtonian result is required for $g_{{}_{N}}\gg a_{{}_{0}}$
(Hees et al., 2014, 2016). However, the precise nature of this convergence is
not important if the relevant quantity is $\bm{g}$ rather than merely its
deviation from $\bm{g}_{{}_{N}}$.

For the computation of $\bm{g}_{{}_{N}}$, we use the standard boundary
condition that the Newtonian potential at large distance is

$\displaystyle\Phi_{N}~{}=~{}-\frac{GM}{r}\,,$ (4)

where $G$ is the Newtonian gravitational constant, $M$ is the total mass
within the simulation volume, and $\bm{r}$ is the position relative to its
barycentre. The boundary condition for the true potential $\Phi$ will be
discussed in Section 2.4 based on Equations 11 and 14.

###  2.2 Treatment of the MW-M31 orbit

An important aspect of our simulations is choosing an appropriate initial
position and velocity for the MW and M31 disc templates to be discussed in
Section 2.3. We do this using a semi-analytic backwards integration (hereafter
SAM) very similar to that used in section 2 of BRZ18, to which we refer the
reader for a detailed discussion of the time evolution of the MW-M31
separation $\bm{d}\left(t\right)$. Briefly, their separation $d$ in physical
coordinates (used throughout this paper) is governed by

$\displaystyle\ddot{d}~{}=~{}\frac{\ddot{a}}{a}d-g+\frac{h^{2}}{d^{3}}\,,$ (5)

where $a$ is the cosmic scale factor, $\bm{h}$ is the angular momentum with
magnitude $h$ and direction $\widehat{\bm{h}}$ (the orbital pole), $g$ is the
radially inward component of the mutual gravity between the MW and M31, and an
overdot denotes a time derivative. The term involving $a$ is present in a
homogeneously expanding universe, but $g$ and $h$ are zero in this case. As
discussed in section 3.1.1 of Haslbauer et al. (2020), we adopt a standard
Planck cosmology (Planck Collaboration XVI, 2014) at the background level
(Table 1). MOND can in principle explain the high locally measured $H_{0}$
(e.g., Riess et al., 2021a; Di Valentino, 2021) as arising from outwards
peculiar velocities induced by the observed KBC void (Keenan et al., 2013), so
our choice of cosmological parameters should be consistent with both early and
late Universe probes of the expansion rate.

Table 1: Our adopted cosmological parameters (Planck Collaboration XVI, 2014), with 0 subscripts denoting present values. We assume a standard flat background cosmology and neglect other components as we are not considering the very early universe. Parameter | Value
---|---
$H_{0}$ | 67.3 km/s/Mpc
$\Omega_{m,0}$ | 0.315
$\Omega_{\Lambda,0}$ | 0.685

The angular momentum barrier $h^{2}/d^{3}$ is necessary to prevent an
unrealistic direct collision between the MW and M31. However, SAM is a timing
argument analysis constrained to give zero peculiar velocity when $a=0.1$,
which implies $h=0$ then. We square this circle by assuming $h=0$ prior to
first apocentre, after which $\bm{h}$ instantaneously jumps to a particular
value that remains fixed until today. The discontinuous behaviour of $\bm{h}$
occurs at a time when the angular momentum barrier is least important to the
trajectory, minimizing numerical effects. Physically, it would be reasonable
if $\bm{h}$ was mostly gained from tidal torques around the time of apocentre,
but the more recent apocentre would be less relevant due to cosmic expansion
driving external perturbers much further from the LG (Section 3.2.1).
Importantly, our primary objective in this contribution is to conduct por
simulations of the flyby. These are initialized 1 Gyr before the flyby, which
is safely after the jump in $\bm{h}$.

The calculation of $g$ is rather complicated $-$ we summarize only the main
points here, and refer the reader to BRZ18 for a detailed discussion. In the
absence of any other bodies and assuming the deep-MOND limit (gravitational
fields $\ll a_{{}_{0}}$), we get a mutual gravity of

$\displaystyle g_{\text{iso}}~{}$ $\displaystyle=$
$\displaystyle~{}\frac{Q\sqrt{GMa_{{}_{0}}}}{d}\,,$ (6) $\displaystyle Q~{}$
$\displaystyle=$
$\displaystyle~{}\frac{2\left(1-{q_{{}_{\text{MW}}}}^{\frac{3}{2}}-{q_{{}_{\text{M31}}}}^{\frac{3}{2}}\right)}{3\,q_{{}_{\text{MW}}}q_{{}_{\text{M31}}}}\,,$
(7)

where $M$ is the total mass of the MW and M31, of which a fraction
$q_{{}_{i}}$ resides in galaxy $i$. The parameter $Q$ accounts for the finite
mass ratio between the MW and M31 (Zhao, Li & Bienaymé, 2010). Roughly
speaking, it is caused by the PDM halo of one galaxy being reduced in mass due
to gravity from the other galaxy, which is a manifestation of the EFE. We
assume that $q_{{}_{\text{MW}}}=0.3$ and $q_{{}_{\text{M31}}}=0.7$
(consistently with BRZ18, ), so $Q=0.7937$.

When the MW and M31 are near pericentre, it is less accurate to assume the
deep-MOND limit. However, the higher $g$ allows us to neglect the relatively
much weaker EFE on the whole LG, so we consider it as an isolated two-body
problem where we get $g$ numerically. Such Newtonian corrections are expected
to have only a minor impact on our results because only a small portion of the
SAM trajectory is subject to them.

When the MW and M31 are near apocentre, we expect the EFE from large-scale
structure to be important (Section 2.4.3). In the limit where
$g_{\text{ext}}\gg g_{\text{iso}}$, we can use the external field
(EF)-dominated analytic solution found by Banik & Zhao (2018a) to get that

$\displaystyle g_{{}_{\text{EFE}}}~{}$ $\displaystyle=$
$\displaystyle~{}\frac{GMa_{{}_{0}}}{d^{2}g_{\text{ext}}}\left(\frac{3+\cos^{2}\theta}{4}\right)\,,$
(8)

where $\theta$ is the angle between $\bm{d}$ and the EF $\bm{g}_{\text{ext}}$.
In reality, $d$ is never large enough for $g_{\text{ext}}$ to dominate,
requiring an interpolation between the isolated and EF-dominated regimes. We
achieve this using equation 14 of BRZ18, which is a fit to numerical results
for the case $\theta=0$.

In addition to considering a uniform $\bm{g}_{\text{ext}}$, we also follow the
BRZ18 approach to include the tidal effect of M81, IC 342, and Cen A. Since in
each case their gravity on the LG is $\ll g_{\text{ext}}$ from more distant
structures (see their section 2.3.1), we can superpose the gravitational field
of each perturber on the LG.444MOND becomes a linear gravity theory if the EFE
dominates (Banik & Zhao, 2018a).

For simplicity, we assume that $\bm{h}$ is constant in SAM except for the
above-mentioned discontinuity at first turnaround, which occurs well before
our por simulation starts. Our por calculations would miss changes in $\bm{h}$
due to tidal torques from perturbers beyond the LG, but we will show later
that such effects are quite small (Section 3.2.1). We therefore expect any
changes in $\bm{h}$ to arise mostly from the EFE (Section 3.2.2) and from
torques on the MW-M31 orbit during the flyby. Neither effect is included in
SAM, but both are included in por. Since the MW-M31 orbit is nearly radial,
changes in $\bm{h}$ should have little effect on the timing argument, which is
the main purpose of SAM.

In summary, we will hereafter use the present MW-M31 separation, direction,
and mutual RV as present-day constraints on the SAM. We vary their orbital
pole and mutual two-body angular momentum magnitude to match the observed SP
orbital poles. Constraints on the PM of M31 will not be taken into account in
this process.

###  2.3 Disc templates

SAM can only tell us the initial position and velocity of the MW and M31
discs. We therefore complement SAM with a system to generate a Milgromian disc
template, to which we then apply the appropriate rotation and Galilean
transformation. We generate two stable isolated Milgromian discs using the
procedures described in Banik et al. (2020), i.e. by using our adapted version
of the Newtonian code disk initial conditions environment (dice; Perret et
al., 2014). The MW and M31 are assumed to have exponential surface density
profiles, motivated by the fact that disc galaxies usually have an exponential
radial profile (Freeman, 1970), which arises naturally with about the right
mass-size relation when spherical gas clouds collapse under MOND gravity
(Wittenburg, Kroupa & Famaey, 2020). We implicitly assume that at the start
time of our simulations $\approx 8$ Gyr ago (redshift $z\approx 1$), the MW
and M31 discs had already formed with nearly their present masses. Thin
rotationally supported disc galaxies do exist at even higher redshift (Lelli
et al., 2018; Neeleman et al., 2020; Rizzo et al., 2020; Lelli et al., 2021).
The relatively isolated nature of the LG (e.g. Banik et al., 2021) suggests
that its major galaxies might well have attained nearly their present mass
rather early in cosmic history, especially in a framework with enhanced long-
range gravity (Peebles & Nusser, 2010) where mergers are less common due to
the lack of dynamical friction between extended CDM haloes (Kroupa, 2015;
Renaud et al., 2016). This is quite plausible given the difficulty faced by
the $\Lambda$CDM paradigm in explaining the high observed fraction of thin
disc galaxies, which could be due to mergers being too frequent in this
framework (Peebles, 2020; Haslbauer et al., 2022).

The orientation, barycentre position, and velocity of each disc are set to the
desired initial values by applying a rotation and Galilean transformation to
its particles using an adapted version of the ramses patch known as condinit.
This also assigns the density and velocity of each gas cell (Teyssier, Chapon
& Bournaud, 2010). To avoid severe thermal effects when the MW and M31
encounter each other, we use the same gas temperature of
$T2\\_ISM=4.65\times{10^{5}}$ K (465 kK) for both galaxies. We set a
temperature floor of $T2\\_star=0.8\,T2\\_ISM$ and disable star formation and
metallicity-dependent cooling, since at this exploratory stage we are mainly
trying to reproduce the observed orientations of the LG SPs. If a suitable
encounter geometry can be found, then it would be worthwhile to conduct a more
detailed simulation with realistic star formation and stellar feedback
prescriptions. However, this is beyond the scope of this work.

Unlike the por simulations of M33 in Banik et al. (2020), an important aspect
of the present por simulations is that we need to consider the outer parts of
the simulated discs in much greater detail because we expect these regions to
be the original source material for the SPs. Indeed, the restricted $N$-body
models of BRZ18 showed that the SPs mainly consist of material at an initial
galactocentric distance of ${\approx 50}$ kpc (see their figures 6 and 7).
Material at such a large distance would be very poorly resolved with a
computationally feasible number of equal mass particles. Therefore, we devise
a procedure to vary the stellar particle mass in our disc templates so as to
maintain a good resolution in the outer parts (see Appendix A). Near the disc
centre, the particle mass is approximately constant at $6\times
10^{5}M_{\odot}$ ($1.4\times 10^{6}M_{\odot}$) for the MW (M31). Each disc
template consists of $5\times 10^{5}$ particles. The spatial resolution of the
por gravity solver also needs to be sufficient $-$ this is discussed further
in Section 2.6. The resolution is 1.5 kpc in the best-resolved regions, though
we show that improving this to 0.75 kpc has little effect on our results
(Appendix E).

####  2.3.1 Initial disc parameters

The rotation curve of each disc template is calculated using MOND gravity with
the present value of $a_{{}_{0}}$, since throughout this work we assume that
$a_{{}_{0}}$ remains constant with time. Therefore, the initial MW and M31
discs would lie on the RAR defined by nearby galaxies. Possible consequences
of a time-varying $a_{{}_{0}}$ were discussed in Milgrom (2015) and section
5.2.3 of Haslbauer et al. (2020), though we do not consider this here for
simplicity. Constraining the evolution is challenging observationally because,
amongst other issues, it is not yet possible to use 21-cm observations of
neutral hydrogen to obtain rotation curves in the distant Universe, so the
H$\alpha$ line is typically used instead. The high-$z$ rotation curve data of
Genzel et al. (2017) are quite consistent with a time-independent
$a_{{}_{0}}$, which sets some limits on its evolution (Milgrom, 2017). Other
works also suggest that high-$z$ galaxies are consistent with MOND (Stott et
al., 2016; Harrison et al., 2017; Genzel et al., 2020; Sharma et al., 2021).
Tighter constraints on this issue would be valuable to better understand the
possible theoretical underpinnings of MOND, and more generally to test its
prediction that isolated galaxies in dynamical equilibrium at any fixed
redshift lie on a tight RAR.

We assume that the present MW and M31 disc scale lengths are as given in table
3 of BRZ18. When setting up their discs at the start of our simulation, we
scale their present lengths by 0.8 for M31 and 0.6 for the MW. Starting with
smaller discs is required to allow them to expand after the flyby to reach a
realistic present-day configuration, and is also in line with the observed
expansion of the stellar component of galaxies over cosmic time (Sharma et
al., 2021). We use the dice setting $Q\\_lim=1.25$ to ensure that all disc
components have an initial Toomre parameter $Q\geq 1.25$, with the MOND
generalization of the Toomre condition (Toomre, 1964) discussed further in
section 2.2 of Banik et al. (2020). The initial parameters of each disc are
summarized in Table 2.

Table 2: Parameters of the MW and M31 discs, each consisting of $5\times 10^{5}$ particles. Both galaxies have an initial gas temperature of 465 kK and a gas fraction of 0.5, which in the MW case is achieved by converting all of its more extended component and part of its less extended component into gas. Thus, the Galactic gas disc has a double exponential profile, while the stellar disc is a single exponential. We adopt an outer limit of 100 kpc for the M31 discs and for the more extended MW component, while 40 kpc is used for its less extended component. This is $\approx 25$ scale lengths in all cases. The initial spin vector of each disc is given in Galactic coordinates $\left(l,b\right)$. Galaxy and | MW | MW | M31 | M31
---|---|---|---|---
component | inner | outer | stars | gas
Total mass | $9.15\times 10^{10}\,M_{\odot}$ | $2.135\times 10^{11}\,M_{\odot}$
Fraction of mass | 0.8236 | 0.1764 | 0.5 | 0.5
Gas fraction | 0.3929 | 1 | 0 | 1
Scale length (kpc) | 1.29 | 4.2 | 4.24 | 4.24
Aspect ratio | 0.15 | 0.0461 | 0.15 | 0.15
Disc spin vector | $\left(55.11^{\circ},\,-84.89^{\circ}\right)$ | $\left(238.38^{\circ},\,-33.71^{\circ}\right)$

Gas dissipation in the tidal tails is likely important to obtaining thin SPs.
While we cannot include this process as rigorously as we would like due to
numerical limitations, it is certainly not appropriate to assume that the pre-
flyby MW and M31 had a similar gas fraction to what we observe. Both discs are
assigned an initial gas fraction of 0.5 because the flyby was $\approx 7-9$
Gyr ago (Zhao et al., 2013), when the gas fractions were likely much higher
than today (Stott et al., 2016). Chemical evolution modelling of the MW
indicates a gas fraction of $\approx 0.5$ at that time (figure 8 of Snaith et
al., 2015). For simplicity, we use the same gas fraction for M31.

The stellar and gas discs of M31 are assumed to have the same scale length of
4.24 kpc, so we request only one disc component in the M31 namelist for the
hydrodynamical version of dice. For the MW, a double exponential profile is
assumed, so two components are defined at that stage. The outer (more
extended) component corresponds to its gas disc today, so we assume this was
entirely gas at the start of our simulation. However, it only comprises 17.64%
of the total MW mass (Banik & Zhao, 2018b), so we also convert a substantial
fraction of the inner (less extended) MW component into gas to get a total gas
fraction of 0.5.

The aspect ratios of the M31 disc and the inner MW disc are set to 0.15, so
the inner MW component has a vertical density profile with a characteristic
$\operatorname{sech}^{2}$ scale height of 193.5 pc. The outer MW disc
component is set to an aspect ratio of 0.0461 so that both components have the
same scale height in pc. For both galaxies, the radial run of the gas disc
scale height is then found by dice to ensure it is as close to equilibrium as
practicable following section 2.3 of Banik et al. (2020). They also described
how particles in the dice template are written out with a reduced mass or not
at all, with the removed mass put back in as gas through the condinit routine
in por to ensure the correct gas fraction. Therefore, our hydrodynamical MOND
version of dice does not yield an equilibrium disc template by itself $-$ it
must be carefully combined with a modified version of por. All the algorithms
we use are publicly available
555https://bitbucket.org/SrikanthTN/bonnpor/src/master/ for reproducibility,
with an accompanying user manual (Nagesh et al., 2021).

Since it is not possible to guarantee that our disc templates are exactly in
equilibrium initially, we start our simulation 1 Gyr before the expected time
of the flyby that we found semi-analytically in Section 2.2. The initial
thickness profiles of the MW and M31 gas discs (shown in Figure 16) are
similar to that used by Banik et al. (2020) in their 100 kK model of M33.

We expect the MW and M31 discs to precess slightly from their initial
orientations (see section 3.1 of BRZ18, ). Therefore, the observed
orientations of the MW and M31 discs differ slightly from our adopted initial
orientations, which are given in final Galactic coordinates $\left(l,b\right)$
in Table 2. We use this system throughout this article when specifying
directions of vectors $-$ it is the system used in the simulation. To
iteratively correct for disc precession, we find the rotation matrix between
the final simulated and observed orientation of each disc, and then apply the
inverse of this rotation to the observed orientation to initialize the next
simulation. We will see later that the final simulated orientation of each
disc agrees quite closely with observations after just one such iterative
correction (Section 3.4).

###  2.4 Adding features to por

The SAM procedure discussed in Section 2.2 is very accurate for following the
overall behaviour of $d\left(t\right)$, but insufficient to model tidal debris
generated by the interaction. This is the main purpose of the por simulations
we will conduct in the present paper. There are some slight differences
between the physics considered in SAM and in por, which we try to rectify by
adding features to por and adjusting the initial conditions.

####  2.4.1 An allowance for tides

Tides from objects outside the LG are not directly included in the por
simulations, but are considered in SAM to estimate the flyby time as
accurately as possible. To approximately include tides in por, we estimate the
amount of energy gained by the MW-M31 system due to tidal compression in the 1
Gyr preceding the flyby, which we estimate directly from SAM using information
on the forces caused by each perturber. This energy is put into the radial
component of $\dot{\bm{d}}$ at the start of our por simulation, which has the
effect of slightly increasing how quickly the MW and M31 are approaching each
other at that time. Our por models neglect the impact of tides after the
flyby, which is justified as the perturbers are much further apart then due to
cosmic expansion.

####  2.4.2 Dark energy

The present MW-M31 separation of 783 kpc (McConnachie, 2012) is large enough
that our por model should include the cosmological term in Equation 5. This
partly consists of a decelerating term due to matter, which is included
automatically because the LG mass mainly resides in the MW and M31, which we
directly include. At late times, there is also an outwards repulsion from dark
energy. We include this while operating ramses in non-cosmological mode, since
this is required by the por patch. For some dark energy parameter
$\Omega_{\Lambda,0}$, the idea is to create an extra repulsive force

$\displaystyle\Delta\bm{g}~{}=~{}{H_{0}}^{2}\Omega_{\Lambda,0}\,\bm{r}\,,$ (9)

where $\bm{r}$ is the position relative to the barycentre.

We can reproduce Equation 9 with a standard Poisson solver if we adjust the
density and boundary condition. Since we want $\bm{g}_{{}_{N}}$ to be
calculated in a standard way so that it is correctly MONDified, we apply the
density increment only to the PDM.

$\displaystyle\Delta\rho_{{}_{\text{PDM}}}~{}=~{}-\frac{3{H_{0}}^{2}\Omega_{\Lambda,0}}{4\mathrm{\pi}G}\,.$
(10)

For consistency, we also adjust the boundary potential for only the MOND stage
by

$\displaystyle\Delta\Phi~{}=~{}-\frac{{H_{0}}^{2}\,\Omega_{\Lambda,0}}{2}\,.$
(11)

Equations 10 and 11 are implemented by appropriate adjustments to the por
algorithm, thereby yielding Equation 9 in the interior.

The impact of dark energy on the MW-M31 dynamics is quite small in MOND as
their mutual gravity is $\approx 30\times$ stronger than the cosmological
acceleration term in Equation 5 (table 10 of BRZ18, ). We nonetheless include
it for completeness. The fact that the cosmological acceleration is small
compared to the internal gravity means that our results are robust with
respect to uncertainty regarding how MOND should be applied in a cosmological
context (the ‘Hubble field effect’ discussed in section 5.2.3 of Haslbauer et
al., 2020). It is however possible to make some plausible assumptions and
simulate systems where the average enclosed density differs only slightly from
the cosmic mean, as done in that work and in several others (e.g. Katz et al.,
2013; Candlish, 2016).

####  2.4.3 The external field effect (EFE)

As discussed in section 2 of BRZ18, the EFE from large-scale structure has a
significant effect on the LG gravitational field when the MW and M31 are close
to apocentre, where they spend a significant amount of time. The EFE is a non-
standard phenomenon caused by the non-linearity of MOND (Equation 1). It leads
to the internal gravitational dynamics of a system being affected by
$\bm{g}_{\text{ext}}$ even in the absence of tidal effects (Milgrom, 1986).
Strong evidence for the EFE in field galaxies was recently reported by Chae et
al. (2020b, 2021) by comparing the RCs of galaxies in isolated and more
crowded environments (building on similar earlier work; Haghi et al., 2016;
Hees et al., 2016).

In the QUMOND approach where we must first get $\bm{g}_{{}_{N}}$, the main
change is to add the Newtonian-equivalent external field to the Newtonian
gravity sourced by the system under study.

$\displaystyle\bm{g}_{{}_{N}}~{}\to~{}\bm{g}_{{}_{N}}+\bm{g}_{{}_{N,\text{ext}}}\,.$
(12)

Assuming that $\bm{g}_{\text{ext}}$ is sourced by a distant point-like object,
we get that

$\displaystyle\bm{g}_{{}_{N,\text{ext}}}\overbrace{\nu\left(g_{{}_{N,\text{ext}}}\right)}^{\nu_{\text{ext}}}~{}=~{}\bm{g}_{\text{ext}}\,.$
(13)

If we know the EF $\bm{g}_{\text{ext}}$, this can be inverted to obtain
$\bm{g}_{{}_{N,\text{ext}}}$.

Once we have adjusted $\bm{g}_{{}_{N}}$ according to Equation 12, we use it to
find $\nabla\cdot\bm{g}$ using Equation 2 as before. The reason is that this
applies to the Newtonian gravity sourced by matter both within and beyond the
simulated domain, the extent of which is an arbitrary decision that should
have no bearing on the result.

Including the EFE also requires us to change the boundary condition, but only
for the MOND stage since the internal dynamics and external field are fully
separable in Newtonian gravity. For simplicity, the boundary should be in some
asymptotic regime where $\bm{g}$ has a well-understood analytic behaviour.
Normally, it is sufficient for the boundary to be distant enough that the
simulated system can be approximated as a point mass. If there is also a non-
negligible EFE, then the simplest option is to choose a boundary where
$\bm{g}_{\text{ext}}$ is much stronger than the internal gravity of the
system. Its internal potential then becomes EF-dominated (e.g., Banik & Zhao,
2018a):

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle-\frac{GM\nu_{\text{ext}}}{r}\left(1+\frac{K_{0}}{2}\sin^{2}\theta\right)\,,$
(14) $\displaystyle K_{0}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\frac{\partial\ln\nu_{\text{ext}}}{\partial\ln
g_{{}_{N,ext}}}\,,$ (15)

where $M$ is the total mass within the simulation volume, $\bm{r}$ is the
position relative to its barycentre, and $\theta$ is the angle between
$\bm{r}$ and $\bm{g}_{\text{ext}}$. Due to the $1/r$ dependence and the fact
that potentials from different sources can be superposed in this perturbative
framework, the result is reminiscent of standard Newtonian mechanics, so the
EF-dominated regime is also known as the quasi-Newtonian regime. Note that
Equation 14 alone is not our boundary condition because we also include a dark
energy adjustment (Equation 11).

For the Newtonian stage of solving the QUMOND Poisson equation, we continue to
use a boundary potential of $-GM/r$, ignoring the constant
$\bm{g}_{\text{ext}}$. This means that our simulations consider the internal
dynamics in a freely falling reference frame accelerating at
$\bm{g}_{\text{ext}}$, whose relevance to the internal dynamics is a violation
of the strong equivalence principle.

For simplicity, we keep $\bm{g}_{\text{ext}}$ fixed over the course of our por
simulation. However, to maximize the accuracy of the overall MW-M31
trajectory, SAM uses a time-dependent $g_{\text{ext}}$ as described in section
2.2 of BRZ18. In both cases, we follow the approach in that paper of assuming
that today, $\bm{g}_{\text{ext}}=0.022\,a_{{}_{0}}$ directed towards Galactic
coordinates $\left(276^{\circ},\,-30^{\circ}\right)$, the direction in which
the LG presently moves relative to the CMB (Kogut et al., 1993). This
direction is assumed fixed throughout cosmic history.

###  2.5 Iterative orbit adjustment

Despite our best attempts to ensure SAM and por handle the flyby problem as
similarly as possible, the two algorithms nonetheless use very different
techniques. Thus, advancing the SAM-generated initial conditions using por
does not quite yield the presently observed $\bm{d}$ or its direction
$\widehat{\bm{d}}$. In the presence of an EFE, the late-time
$\widehat{\bm{d}}$ affects the torque exerted by the EFE on the MW-M31 system
(Section 3.2.2), which in turn influences the present PM of M31. Moreover, an
incorrect final $\widehat{\bm{d}}$ suggests that the flyby took place in a
different orientation to how it occurred in the simulation, which would
influence the SPs.

To ensure the final $\widehat{\bm{d}}$ matches the observed sky position of
M31 as accurately as possible, we find the rotation matrix that takes the
observed $\widehat{\bm{d}}$ to the simulated final value. We then rerun SAM
with the inverse of this rotation applied to the presently observed
$\widehat{\bm{d}}$ and the relative velocity $\dot{\bm{d}}$ used in the
previous SAM simulation. The idea is that if $\widehat{\bm{d}}$ ends up
$10^{\circ}$ further north than observed, then we can get approximately the
correct final $\widehat{\bm{d}}$ by running SAM with a present
$\widehat{\bm{d}}$ that lies $10^{\circ}$ south of the actually observed
direction towards M31, because there is some additional physical effect in por
but not in SAM that pushes $\widehat{\bm{d}}$ further north by $10^{\circ}$.
The initial conditions generated in this way are used to rerun the por
simulation. We find that just one such iteration allows us to match the
observed $\widehat{\bm{d}}$ to within a few degrees, which we consider
sufficient.

###  2.6 Simulation setup and initial conditions

####  2.6.1 SAM

We begin by running SAM with the constraints given in Table 3. The main model
parameters that we vary are:

  1. 1.

the MW-M31 orbital pole $\widehat{\bm{h}}$, and

  2. 2.

the magnitude $h$ of their mutual two-body angular momentum.

As discussed in section 5.1.3 of BRZ18, we account for reflex motion induced
by the Large Magellanic Cloud (LMC) on the MW and by M33 on M31, leading to a
present RV of $\dot{d}=-93.4$ km/s.666Although the LMC should form out of
tidal debris expelled during the flyby, our model does not form individual
TDGs. The simulated Galactic disc thus does not experience recoil from an
orbiting massive satellite, as it might do in a more advanced model. Since the
timing argument is sensitive to the gravity between the MW and M31, any
massive satellites should be included in the mass and velocity of each galaxy.
We keep $\dot{d}$ fixed, but vary the present tangential velocity
$\bm{v}_{\text{tan}}$ used in SAM, which does a backwards integration.
Changing $h$ mainly affects the perigalacticon distance, but does not much
alter the relative speed then. Varying $\widehat{\bm{h}}$ allows us to
consider a wide range of possible orbital geometries. While our priority is to
match the SP orientations, we subsequently compare our best-fitting model with
the observed PM of M31 (Section 3.2). Interestingly, if we neglect both
considerations, the timing argument alone sets some constraints on
$\widehat{\bm{h}}$ because perturbers like Cen A have a different influence on
the internal dynamics of the LG depending on $\widehat{\bm{d}}$ prior to the
flyby. As a result, the timing argument mass of the LG is too high for a large
range of $\widehat{\bm{h}}$, which renders the models unlikely as the MOND
timing argument mass should correspond to the baryonic mass (section 5.1.3 of
BRZ18, ).

Table 3: Present-day constraints imposed on SAM, which does a backwards integration. The tangential velocity of M31 is varied to best match the observed SP orbital poles. Initial MW-M31 $\ldots$ | Value
---|---
distance | 783 kpc
direction | $\left(121.17^{\circ},-21.57^{\circ}\right)$
RV |  $-93.4$ km/s

The model parameters which we explore are the same as in BRZ18, except that we
do not vary the EF. It was indeed shown that varying the EF within the
plausible range has only a small impact on the final results, which is likely
due to the well-known LG velocity in the CMB frame setting some constraint on
the history of $\bm{g}_{\text{ext}}$. We therefore use the best-fitting model
of BRZ18 in which $\bm{g}_{\text{ext}}=0.022\,a_{{}_{0}}$ directed towards
Galactic coordinates $\left(276^{\circ},\,-30^{\circ}\right)$ at the present
time. To make the MW-M31 trajectory as realistic as possible, the
$g_{\text{ext}}$ assumed in SAM varies with time, as discussed in their
section 2.2. However, implementing a time-varying $\bm{g}_{\text{ext}}$ in por
would involve significant complications, so this uses the present
$\bm{g}_{\text{ext}}$ for the full duration of the simulation. The EFE is not
expected to influence details of the MW-M31 interaction, but causes them to
reach a larger apocentre at late times by weakening their mutual gravity. As a
result, the main role of the EFE is to alter the timing of the MW-M31 flyby,
which sets the time available for the tidal debris to settle down after the
flyby. We therefore prioritize making the flyby time as accurate as possible
in our por models, which requires a carefully prepared SAM. Due to these
slight differences and the more rigorous treatment of the flyby in por, the
MW-M31 trajectory is expected to differ somewhat compared to SAM. We do not
run the por simulations for a different length of time than that indicated by
SAM in order to obtain a better agreement with the observed distance to M31,
so this is expected to differ somewhat from the observed $783\pm 25$ kpc
(McConnachie, 2012).

It is important to realize that although the final M31 distance in por need
not match observations exactly, requiring even an approximate match places
non-trivial constraints on our model. This is because a very close encounter
would lead to significant dynamical friction between the baryonic discs,
causing a very low apocentre and a subsequent merger within a few Gyr. This
possibility was neglected in the models of BRZ18, where the $d\left(t\right)$
returned by SAM was assumed to be exactly correct. In general, a
hydrodynamical model of the interaction is obviously a significant advance on
the previous restricted $N$-body models, even if the much higher computational
cost reduces the scope for fine-tuning to match certain observables to high
precision. We tried $\approx 30$ por models to obtain a good fit to the SP
orientations, which we judged visually.777We also ran a few more models to
fine-tune the agreement with some observables like the present disc
orientations and the MW-M31 direction (Section 2.5).

Table 4: The initial position in kpc and velocity in km/s of the MW and M31 disc centres at the start of our best-fitting por simulation. The combined barycentre is at the origin because the MW:M31 mass ratio is 3:7 (Table 2), consistently with BRZ18. SAM indicates that these initial conditions are valid 8.1 Gyr ago. Galaxy | MW | M31
---|---|---
Direction | Position | Velocity | Position | Velocity
$x$ | -56.0 | 5.2 | 24.0 | -2.2
$y$ | 254.2 | -177.8 | -109.0 | 76.2
$z$ | 44.1 | -86.7 | -18.9 | 37.2

####  2.6.2 por

Our por simulations use very similar settings to those described in section
2.4 of Banik et al. (2020), so we briefly mention the main points here. We use
the non-cosmological particle-in-cell mode, activate MOND and the EFE, and use
the refinement conditions $m\\_refine=10^{3}M_{\odot}$ and
$n\\_subcycle=\left(1,1,2,2\right)$. Since we are simulating the whole LG, we
use a much larger box size of 6144 kpc, with $7-12$ levels of refinement. The
most poorly resolved regions thus have a resolution of $6144/2^{7}=48$ kpc,
which improves to $6144/2^{12}=1.5$ kpc for the best resolved regions.
Teyssier (2002) provides a more detailed description of ramses, including
default values of parameters that we do not alter. Perhaps the most important
of these is the $gravity\\_type$, whose default setting of 0 represents self-
gravity. We also use the default Poisson convergence parameter
$epsilon=10^{-4}$.

The results of our best-fitting model are discussed next. In this model, the
initial position and velocity of the MW and M31 are as given in Table 4 for
the centres of two galaxies each consisting of $5\times 10^{5}$ particles and
a gas fraction of 50%, with the disc orientations given in Table 2. The
initial gas temperature is $T2\\_ISM=465$ kK, with the temperature floor
$T2\\_star$ being 20% lower. According to SAM, the initial MW and M31
positions correspond to 1 Gyr before the flyby, so our por simulations can be
considered to start 5.72 Gyr after the Big Bang (8.1 Gyr ago). Starting the
simulations 1 Gyr before the flyby gives the discs some time to settle down
first, but avoids the need to simulate very early epochs at which the
cosmological term in Equation 5 would be significant.

##  3 Results

In this section, we present the results of our best-fitting model and compare
it with relevant observations. We use the simulation snapshot after 8.2 Gyr
because this corresponds as closely as possible to the lookback time estimated
by SAM for the initial conditions of Table 4 .

###  3.1 The MW-M31 orbit

The initial conditions are obtained by running SAM, which gives a timing
argument mass of $3.457\times 10^{11}\,M_{\odot}$ for the whole LG. This is
somewhat higher than the sum of the disc masses used in por (Table 2), which
we consider acceptable as there would also be some mass in, e.g., a halo of
gas around each galaxy, though we do not include a halo explicitly. This issue
was discussed further in section 5.1.1 of BRZ18, who argued that although the
observed MW and M31 RCs suggest a combined disc mass of $2.3\times
10^{11}\,M_{\odot}$, a modestly higher timing argument mass for the whole LG
is quite feasible due also to satellite galaxies like the LMC and M33.
Equation 1 implies that a 50% increase in the mass of a galaxy increases the
flatline level of its RC by just 11%. Such a small increase in the MW and M31
RCs at large ($\ga 100$ kpc) distances is difficult to rule out at present,
especially for M31 (Corbelli et al., 2010; Sofue, 2015).

Figure 1: The stellar particles in the central 500 kpc square of our por
simulation, viewed along the direction which would make both discs appear
edge-on in actual observations at $z=0$ (Section 3.4). Since the disc
orientations change somewhat during the flyby and we have not perfectly
adjusted for this, the discs are not perfectly edge-on in the last snapshots.
The projected surface density is shown in units of $1000\,M_{\odot}$/pc2. The
time since the start of the simulation is indicated at the top left of each
panel. The MW disc appears almost horizontal in this view. Notice that the MW
and M31 retain thin discs, outside of which there is very little material.
Figure 2: Similar to Figure 1, but now showing the gas. The smallest simulated
gas cell has sides of length 1.5 kpc.

Before conducting more detailed analyses, we show the central 500 kpc of the
por simulation as viewed from the direction which would make both discs appear
perfectly edge-on if their orientations are as observed $-$ this is very
nearly the case (Section 3.4). Figure 1 shows the stellar particles, while
Figure 2 shows the gas. The discs undergo closest approach $\approx 1$ Gyr
into the simulation. Dynamical friction during the flyby is small, allowing
the galaxies to reach a large post-encounter separation. Moreover, the flyby
does not disrupt the MW and M31 discs too severely, allowing them to retain a
thin disc by the end of our simulation $-$ we will investigate this in more
detail in Section 3.4. There is also a small numerical drift of the MW-M31
barycentre, which can be understood using much less computationally intensive
methods (see footnote 14 of Banik et al., 2020).

We extract the MW-M31 trajectory from our por simulation while it is running,
without extracting every simulation output. We briefly describe this technique
in Appendix C as it could also be useful for other projects and is part of the
publicly available version of por used here. In this way, we obtain the
trajectory shown in Figure 3, with the top and bottom panel used to show the
MW-M31 separation and relative velocity, respectively. The expected result
using SAM is shown as a dotted red line in each panel. Both methods give a
rather similar overall trajectory, indicating little dynamical friction during
the encounter. This is due to the fairly large pericentre distance of 81 kpc.
Table 5 summarizes information about the pericentre and apocentre in each
trajectory. The higher second apocentre in SAM was discussed in section 5.1.2
of BRZ18, who concluded that it is mainly driven by tides from Cen A due to
its relatively high mass and close alignment with the MW-M31 line after but
not before the flyby. Tides are not explicitly included in our por simulation,
which moreover starts shortly before the flyby and thus includes only the most
recent MW-M31 apocentre, when the perturbers are rather distant.

Figure 3: _Top_ : The MW-M31 separation as a function of time in our por (solid black) and SAM (dot-dashed red) models. Although our por simulation is not cosmological, its initial conditions are obtained from SAM 1 Gyr before the flyby using a total mass of $3.457\times 10^{11}\,M_{\odot}$, thus placing the por model in the cosmological context of a trajectory consistent with the timing argument (Section 2.2). The por model (with a slightly lower total mass as it considers only the discs) is advanced slightly into the future, but our detailed analyses focus on the snapshot 8.2 Gyr after the start of the simulation to most closely match the post-flyby time estimated by SAM. The vertical grey line in each panel shows the start of our por simulation. _Bottom_ : The MW-M31 relative velocity in SAM (dot-dashed red) and por (solid black). The solid blue line with more scatter shows their relative angular momentum as an equivalent tangential velocity at their present separation of 783 kpc (McConnachie, 2012), multiplied by 10 for clarity. The discontinuity in the SAM velocity around 3 Gyr arises from assuming no angular momentum prior to first apocentre, to allow the MW-M31 trajectory to reach zero peculiar velocity at very early times (Section 2.2). Notice the lack of a significant change in the angular momentum due to the MW-M31 flyby, indicating that it hardly alters their orbit. The gradual change over many Gyr is due to torque from the EFE (Section 3.2.2) and due to tidal debris falling back onto each galaxy, especially the MW (Section 3.3.3). Table 5: Information about the MW-M31 pericentre and apocentre in our SAM and por simulations, based on Figure 3. Times are relative to the Big Bang. Orbital phase | Perigalacticon | Apogalacticon
---|---|---
Algorithm | SAM | por | SAM | por
Time (Gyr) | 6.72 | 6.79 | 12.0 | 12.6
Distance (kpc) | 73.8 | 81.5 | 867.8 | 871.2
Velocity (km/s) | 530.9 | 452.1 | 45.1 | 18.6

The simulation snapshot that we analyse (8.2 Gyr after the start) has $d=846$
kpc, which slightly exceeds the observed $783\pm 25$ kpc distance to M31
(McConnachie, 2012). As explained previously, we do not rectify this issue by
running the simulation for longer. Doing so should have only a small effect on
the SPs, the main focus of this work. Importantly for the overall geometry,
the simulated $\widehat{\bm{d}}=\left(121.11^{\circ},\,-21.99^{\circ}\right)$,
which differs by only $0.42^{\circ}$ from the observed sky position of M31
$\left(121.17^{\circ},\,-21.57^{\circ}\right)$.888The simulated
$\widehat{\bm{d}}$ is actually found using the method discussed in Section 3.2
to identify the MW and M31 barycentres, but the results are almost identical.
Therefore, the overall MW-M31 trajectory in our por simulation is quite
reasonable, and should be consistent with cosmological initial conditions (the
timing argument) at much earlier times.

###  3.2 Final tangential velocity and PM

The PM of M31 provides an important constraint on the orbital geometry of our
best-fitting model. When scanning the parameter space (Section 2.6.1), we did
not consider this constraint, though we did not consider all possible orbital
geometries either as some are highly disfavoured by the timing argument alone
(Section 2.2). In the following, we describe how we obtain the predicted PM of
M31, and compare this to the latest observational constraints.

To find the separation and relative velocity of the MW and M31, we first have
to identify each galaxy’s centre of mass. We obtain an initial guess using the
iterative on-the-fly method described in Appendix C based on the particles
alone. Much more detailed analyses are possible using a simulation snapshot
because we also use our modified version of rdramses (section 3 of Banik et
al., 2020) to obtain a list of all gas cells, treating them as particles at
the cell centres. This allows the stars and gas to be analysed on an equal
footing. We therefore find the barycentre of all material whose position and
velocity lies within 250 kpc and 500 km/s, respectively, of our initial guess
for the barycentre. These thresholds are deliberately set quite wide to reduce
the risk of converging on the wrong density maximum. This process is repeated
iteratively until the barycentre position shifts by $<1$ kpc and its velocity
shifts by $<1$ km/s between successive iterations. The process converges very
rapidly because the initial guess for the barycentre based on particles alone
(Appendix C) is already highly accurate.

Table 6: The final position in kpc and velocity in km/s of the MW and M31 barycentres 8.2 Gyr into our best-fitting por simulation. According to SAM, this is the por simulation snapshot temporally closest to the present epoch. The combined barycentre is slightly offset from the origin due to numerical drift, as discussed further in footnote 14 of Banik et al. (2020). Galaxy | MW | M31
---|---|---
Direction | Position | Velocity | Position | Velocity
$x$ | 232.9 | -31.0 | -172.1 | -6.6
$y$ | -248.8 | 73.5 | 422.6 | 29.4
$z$ | 176.4 | -0.7 | -140.2 | -14.4

Table 6 shows the position and velocity of the MW and M31 in the simulation
reference frame, from which we get the MW-M31 relative separation and
velocity. We take the cross product of these vectors to get their orbital
angular momentum $\bm{h}$ and thus the orbital pole $\widehat{\bm{h}}$. We
also use $h$ to obtain the tangential velocity $v_{t,783}$ at a distance of
783 kpc, which slightly exceeds the simulated tangential velocity by a factor
of (846/783) because the final distance is slightly larger than observed. We
then find what PM components of M31 on the night sky would most closely mimic
$\bm{d}$ and $\bm{v}_{\text{tan}}$ in the analysed por simulation output. We
do this by adjusting the assumed M31 PM components
$\left(\mu_{\alpha,\star},\mu_{\delta}\right)$ in SAM to find the combination
which best matches the
$\widehat{\bm{h}}=\left(209.77^{\circ},\,3.29^{\circ}\right)$ and
$v_{t,783}=34.12$ km/s found from por, with $\mu_{\alpha,\star}$ and
$\mu_{\delta}$ being the angular velocity in the directions which most quickly
increase the right ascension and declination, respectively. In this
calculation, the present M31 direction and RV in SAM are fixed to the values
in Table 3. The Galactic circular velocity at the Solar circle is assumed to
be 239 km/s (McMillan, 2011), while the non-circular motion of the Sun is
taken from Francis & Anderson (2014) $-$ uncertainties in these parameters and
in the distance to M31 are much smaller than in its PM. We vary
$\left(\mu_{\alpha,\star},\mu_{\delta}\right)$ using a gradient descent
algorithm (Fletcher & Powell, 1963) to minimize the sum of squared errors in
$\widehat{\bm{h}}$ and $v_{t,783}$, with each error scaled to an uncertainty
of $1^{\circ}$ and 1 km/s, respectively. While it is possible to match
$v_{t,783}$ exactly, the best-fitting
$\left(\mu_{\alpha,\star},\mu_{\delta}\right)$ still gives a small error in
$\widehat{\bm{h}}$ because of a slight difference in the final
$\widehat{\bm{d}}$ between por and SAM. This is reduced by an iterative rerun
of the simulation. Since the mismatch is then only $0.42^{\circ}$, our
approach gives a good idea of what our simulation implies for the present M31
PM.

Figure 4: The heliocentric PM of M31 in our best-fitting por model
(intersection of blue and pink lines), found using the method described in
Section 3.2. The short blue (pink) line indicates the effect of varying the
LMC (M33) mass by $\pm 20\%$. In both cases, a heavier satellite increases the
predicted M31 $\mu_{\alpha,\star}$. The open blue circle shows the M31 PM for
a purely radial MW-M31 orbit, while the black cross shows the result of
reversing the Galactocentric tangential velocity of M31 given by our por
model. The large red error bars show the observed M31 PM using HST (van der
Marel et al., 2012), while the black error bars use instead Gaia eDR3 (Salomon
et al., 2021). The measurement of van der Marel et al. (2019) is not shown
here for clarity, but all three PMs are listed in Table 7. The horizontal
green line represents 10 km/s at the M31 distance of 783 kpc.

Our result is shown in Figure 4. As discussed in section 5.1.3 of BRZ18, this
includes a correction for the LMC and M33 altering the barycentric velocity of
the MW and M31, respectively. Though the LMC and M33 positions and velocities
are known fairly well, there is some uncertainty regarding their mass. This
affects the results slightly because our model does not form individual TDGs,
so it does not properly capture the induced recoil on the MW from the LMC,
which ideally should form in a simulation of the flyby. The uncertainties are
rather small in a MOND context because all galaxies are purely baryonic and we
know the position and velocity of M33 and the LMC rather well. We therefore
vary the mass of each satellite by $\pm 20\%$ and repeat the above-mentioned
PM calculation, which assumes the MW (M31) position and velocity in SAM refers
to the barycentre of the MW-LMC (M31-M33) system. As discussed in section 4.4
of Banik & Zhao (2016), the PM correction due to the LMC is not as significant
because its Galactocentric velocity is mostly away from M31. This is evident
in the very short length of the dark blue line through the model PM in Figure
4, which indicates the impact of a 20% uncertainty in the LMC mass. The longer
pink line shows the same for M33, whose effect is somewhat larger. Even so,
the uncertainty on the M31 PM is $<1\,\mu$as/yr in both cases, which is much
smaller than the difference between its predicted PM and that for a purely
radial MW-M31 orbit (open blue circle).

Table 7: The PM of M31 in $\mu$as/yr according to different investigations, which in order are: HST alone (van der Marel et al., 2012), HST and Gaia DR2 (van der Marel et al., 2019), and Gaia eDR3 (Salomon et al., 2021) using the $B_{\text{pm}}$ sample in their table 4. The latter result is in a heliocentric frame, including instrumental corrections for PMs of background quasars and using an outlier rejection system. Astrophysical motions beyond the Solar System are not accounted for in these results, since we do the corrections ourselves. The final column gives the PM in our best-fitting model. |  | HST $+$ |  |
---|---|---|---|---
Component | HST | Gaia DR2 | Gaia eDR3 | Model
$\mu_{\alpha,\star}$ | $~{}45\pm 13$ | $~{}49\pm 11$ | $~{}48.0\pm 10.2$ | 43.4
$\mu_{\delta}$ | $-32\pm 12$ | $-38\pm 11$ | $-38.4\pm 7.8$ | $-33.2$

We compare the predicted M31 PM with the measurements in Table 7. For clarity,
we only plot the observed values from the Hubble Space Telescope (HST) prior
to Gaia (van der Marel et al., 2012) and the recent result from Gaia early
data release 3 (eDR3; Salomon et al., 2021). In the latter case, we use their
so-called $B_{\text{pm}}$ sample in their table 4 to include corrections for
instrumental effects as estimated from PMs of background quasars, but without
including corrections for astrophysical motions beyond the Solar System $-$
these are handled in our analysis (see, e.g., section 2.3 of Banik & Zhao,
2016). The combination of HST results with Gaia data release 2 (DR2; van der
Marel et al., 2019) is not shown in Figure 4 as it gives a similar but
somewhat less precise result to the $B_{\text{pm}}$ sample in Gaia eDR3. All
three PM estimates are mutually consistent within uncertainties. They also
agree very well with the PM of our por model (Figure 4). The PM uncertainties
are now small enough that this is not guaranteed. We demonstrate this by
showing the expected PM of M31 if the por-determined Galactocentric tangential
velocity of M31 is reversed, which reverses $\bm{h}$. The result is the black
cross towards the upper left of Figure 4, which is now inconsistent with the
latest determination of the M31 PM (Salomon et al., 2021).

Figure 5: The probability of a higher $\chi^{2}$ between the observed PM of
M31 and in our por model if the predicted Galactocentric tangential velocity
of M31 is rotated around the simulated direction towards M31 by the angle
shown on the $x$-axis. There is a clear preference for rotation angles close
to zero, i.e. for an MW-M31 orbital pole similar to that in the model.

To explore this further, we rotate the por-determined $\widehat{\bm{h}}$
around the M31 direction $\widehat{\bm{d}}$ by all possible angles between
$-180^{\circ}$ and $+180^{\circ}$ in steps of $1^{\circ}$. We then determine
the $\chi^{2}$ statistic relative to the results of van der Marel et al.
(2012) and Salomon et al. (2021). Since this $\chi^{2}$ is based on two
parameters (the PM components), the probability of a higher $\chi^{2}$ arising
by chance if the model were correct is $\exp\left(-\chi^{2}/2\right)$. Figure
5 shows this probability for all possible rotation angles, with $0^{\circ}$
corresponding to the actual $\widehat{\bm{h}}$ of our por model. It is clear
that for fixed $h$, the direction $\widehat{\bm{h}}$ preferred by the observed
M31 PM (Table 7) corresponds quite closely to that of our best-fitting por
model.

At this point, it is worth emphasizing that we did not consider the PM of M31
when selecting the best-fitting por simulation $-$ this was based purely on
the phase space distribution of the tidal debris. The agreement between the
observed M31 PM and in this model constitutes a non-trivial success thereof.

####  3.2.1 Tidal torques

Our por model neglects the late-time effect of tidal torques from perturbers
outside the LG, which would also affect the PM of M31. To estimate the tidal
torque from each perturber, we assume the deep-MOND limit and that the EFE
dominates. The former assumption is clear because of the large distances to
the perturbers, while the latter assumption was justified in section 2.3.1 of
BRZ18. As we are only interested in a rough estimate, we approximate that the
gravity towards any object of mass $M$ is
$GMa_{{}_{0}}/\left(r^{2}g_{\text{ext}}\right)$ at distance $r$ (Equation 8),
neglecting the small angular dependence. Using also the distant tide
approximation and assuming that $g_{\text{ext}}=0.022\,a_{{}_{0}}$, the
relative MW-M31 acceleration in the direction orthogonal to their separation
has a magnitude

$\displaystyle
g_{\text{tan}}~{}=~{}\frac{3\,GMa_{{}_{0}}d\cos\theta\sin\theta}{r^{3}g_{\text{ext}}}\,,$
(16)

where $r$ is the distance from the MW-M31 mid-point to the perturber,
$\cos\theta$ is the angle between this direction and $\widehat{\bm{d}}$, and
$M$ is the mass of the perturber. $r$ is measured from the MW-M31 mid-point
because the MW and M31 are treated as independent test particles freely
falling in the gravitational field of the perturber, so the individual MW and
M31 masses are irrelevant. Thus, the accuracy of our calculations would be
maximized if using the MW-M31 mid-point. For a perturber with $r\gg d$, it
does not matter exactly which point we use as the ‘centre’ of the LG.

Table 8: Our estimated tidal torque due to each perturber, defined as the change it causes in the MW-M31 relative tangential velocity (Equation 16). The assumed perturber properties are as listed in table 2 of BRZ18. Despite the high mass of Cen A, it has only a small effect because it is almost on the MW-M31 line. The torquing effect of a perturber is maximized if $\left|\cos\theta\right|\approx 1/\sqrt{2}$, which is almost the case for IC 342. However, even its effect is substantially smaller than the PM uncertainty of $\approx 50$ km/s (Figure 4). Perturber | $\left|\cos\theta\right|$ | Change in tangential
---|---|---
velocity over 5 Gyr (km/s)
Cen A | 0.99 | 2.43
M81 | 0.32 | 3.79
IC 342 | 0.76 | 8.22

Assuming the same perturber properties as listed in table 2 of BRZ18 and that
$d=783$ kpc, we obtain the tidal torque estimates given in Table 8 over a
period of 5 Gyr, which is roughly the amount of time for which $d$ is similar
to its maximum value (Figure 3). Tidal torques would be much less significant
around the time of the flyby. As shown by the values in Table 8 and the green
line in Figure 4 representing 10 km/s, it is clear that tidal torques from
perturbers hardly affect the present PM of M31. Its measurement accuracy would
need to improve another order of magnitude for such subtle effects to be
discernible, which would then necessitate more detailed modelling. The small
effect of tidal torques arises mainly because of the isolated nature of the
LG, but also because one of the most massive external galaxies which could
tidally affect it (Cen A) lies almost on the MW-M31 line (Ma et al., 1998),
limiting its tidal torque on the LG.

####  3.2.2 The external field torque (EFT)

In MOND, the mutual gravity between two masses does not necessarily align with
their separation if there is an external field (e.g., Banik & Zhao, 2018a).
This creates a non-tidal torque on the MW-M31 system due to objects outside
it. Since this is induced by the EFE, we call this the EFT. As described in
appendix A of BRZ18, the tangential component of the relative gravity between
two masses due to the EFT can be estimated as

$\displaystyle g_{\text{tan}}~{}=~{}$
$\displaystyle\sin\theta\left[\cos\theta+\frac{4}{5}\sin^{2}\theta\cos\left(\mathrm{\pi}q_{{}_{\text{MW}}}\right)\frac{\tilde{r}}{1+\tilde{r}^{2}}\right]$
(17)
$\displaystyle\times\frac{GMa_{{}_{0}}r}{2g_{\text{ext}}\left(r^{3}+{r_{t}}^{3}\right)}\,,$

with the “$-$” sign changed to “$+$” in front of the $4/5\ldots$ (see section
4.4 of Banik et al., 2020). $\widetilde{r}$ is the separation relative to that
at which the problem becomes EF-dominated:

$\displaystyle\tilde{r}~{}\equiv~{}\frac{r}{r_{t}}\,,\quad
r_{t}~{}=~{}\frac{\sqrt{GMa_{{}_{0}}}}{Q\,g_{\text{ext}}}\,.$ (18)

$\theta$ is the angle between $\bm{d}$ and $\bm{g}_{\text{ext}}$, while
$Q=0.7937$ for an MW mass fraction of $q_{{}_{\text{MW}}}=0.3$ (Equation 7).
Note that in appendix A of BRZ18, the two masses considered were an external
perturber (e.g. Cen A) and the whole LG, with the latter placed at the
position of the MW or M31 depending on the circumstances. Their derivation
leads directly to Equation 17 if instead we take the two masses to be the MW
and M31.

The present work improves significantly upon BRZ18 in that the EFT is directly
included in por (Section 2.4.3). However, this still relies on knowing the
appropriate value of $\bm{g}_{\text{ext}}$ and its even more uncertain time
dependence. We can use Equation 17 to estimate the significance of these
factors. Since the torque is small for an isolated system ($\widetilde{r}\ll
1$), it is mainly important when the MW and M31 are close to apocentre. If we
use $r=900$ kpc and integrate Equation 17 over a 5 Gyr period (as in Section
3.2.1) assuming constant $g_{\text{tan}}$, we get a tangential impulse of
$\Delta v_{t,783}=33.6$ km/s for $g_{\text{ext}}=0.03\,a_{{}_{0}}$, including
also a factor of 900/783 to account for changes in angular momentum due to the
EFT having a greater impact on the present tangential velocity of M31 if it is
currently closer to us. If instead we assume that
$g_{\text{ext}}=0.02\,a_{{}_{0}}$, then $\Delta v_{t,783}$ is only 25.9 km/s.
This difference of $\approx 8$ km/s is also apparent at other separations,
e.g. if we use $r=600$ kpc, we get that $\Delta v_{t,783}$ differs by 11 km/s
depending on whether $g_{\text{ext}}$ is 2% or 3% of $a_{{}_{0}}$.

In principle, the functional form of the dependency on $q_{{}_{\text{MW}}}$ is
not known as BRZ18 only had numerical results for the case
$q_{{}_{\text{MW}}}=0$. This means we can replace the factor of
$\cos\left(\mathrm{\pi}q_{{}_{\text{MW}}}\right)$ by any function with value 1
when $q_{{}_{\text{MW}}}=0$ and which is anti-symmetric with respect to
$q_{{}_{\text{MW}}}\to 1-q_{{}_{\text{MW}}}$. However, due to the low value of
$\sin^{2}\theta=0.167$, these results hardly change if instead we use a mass
ratio dependency of $\left(1-2q_{{}_{\text{MW}}}\right)$ or
$\cos^{3}\left(\mathrm{\pi}q_{{}_{\text{MW}}}\right)$ instead of
$\cos\left(\mathrm{\pi}q_{{}_{\text{MW}}}\right)$.

Figure 6: The MW-M31 orbital pole in our por simulation, shown as Galactic
latitude (blue) and Galactic longitude less $180^{\circ}$ (red). The deviation
from the initial orbital pole is shown using a solid black line. The vertical
solid grey lines around 7 Gyr and 13 Gyr show the times of perigalacticon and
the most recent apogalacticon, respectively, while the dotted grey line shows
the present time of 13.82 Gyr. Notice the lack of precession around
pericentre. Combined with the magnitude of the MW-M31 angular momentum hardly
changing then (bottom panel of Figure 3), it is clear that the MW-M31
interaction has little torquing effect on their orbit. Their orbital pole
precesses more significantly around apocentre due to the EFT (Section 3.2.2).
The present orbital pole direction
$\left(209.77^{\circ},\,3.29^{\circ}\right)$ is very similar to that of the
M31 SP.

Since $r$ around the time of apocentre typically has some intermediate value
between 600 kpc and 900 kpc (Figure 3), we estimate that uncertainty in the
behaviour of $g_{\text{ext}}$ causes an $\approx 10$ km/s uncertainty in the
PM of M31. This would not much alter the present MW-M31 orbital pole because
$v_{t,783}=34.12$ km/s in our best-fitting por model. Indeed, Figure 6 shows
that $\widehat{\bm{h}}$ does not change much due to the EFT already included
in por, so uncertainties in the EFT should have only a small effect.
Interestingly, $\widehat{\bm{h}}$ also hardly precesses due to the flyby,
which is indicated with a vertical solid grey line at pericentre (the later
solid grey line represents apocentre). As a result, the final
$\widehat{\bm{h}}$ is fairly similar to its initial orientation. The more
rapid precession of $\widehat{\bm{h}}$ after the most recent apocentre is
caused by the MW-M31 orbit becoming almost radial (notice the low value of $h$
at that time in Figure 3). It could also be related to tidal debris falling
back onto the MW (Section 3.3.3).

Uncertainty in the EFT is mitigated by the geometric factors being well known:
the MW-M31 line at apocentre must be quite close to its presently observed
direction as they are on a nearly radial orbit (van der Marel et al., 2012;
van der Marel et al., 2019; Salomon et al., 2021). $\bm{g}_{\text{ext}}$ is
also constrained by the observed motion of the LG with respect to the
CMB/surface of last scattering (section 2.2 of BRZ18, ). Since the EF on the
LG mostly arises from rather distant sources, the direction of
$\bm{g}_{\text{ext}}$ would not have changed much in the last 5 Gyr. As a
result, the EFT at late times can only affect the PM of M31 along one
particular direction. However, it is not too useful to speculate further about
this because the PM of M31 still has an uncertainty $\gg 10$ km/s (Figure 4).

###  3.3 Tidal debris

We now turn to the distribution of tidal debris around the MW and M31 disc
remnants. One of the main goals of our por simulations is to check whether the
tidal debris around each galaxy prefers a particular orbital pole, and if so,
to compare this with the observed orbital pole of its SP. This requires us to
define a ‘satellite region’ around each disc, which we take to be within 250
kpc of the barycentre found in Section 3.2 and $>z_{\text{max}}=50$ kpc from
the disc plane. For simplicity, we assume that the orientation of each disc
matches observations $-$ we show later that this is a fairly good assumption
(Section 3.4) due to iterative adjustments to the initial orientation of each
disc (Section 2.3.1). Even if the actual disc orientations differ somewhat,
our choice of $z_{\text{max}}$ should completely exclude the disc. However,
only a small portion of the SPs would be lost as the satellite galaxies of the
MW and M31 go out much further (see, e.g., figure 1 of Pawlowski, 2018). We
note that using a spherical excluded region for the disc does not work very
well since e.g. excluding the central 30 kpc still leaves a considerable
amount of material close to the disc plane at larger radii. This makes it very
difficult to clearly disentangle the disc and SP, at least unless a much
larger inner radius is considered $-$ which then loses much of the satellite
region. We therefore consider only a disc-shaped excluded region in what
follows, or alternatively focus on precisely this region when analysing the
disc remnants (Section 3.4). Appendix D shows the effect of varying
$z_{\text{max}}$.

####  3.3.1 Orbital pole distribution and mass

Figure 7: _Left_ : The orbital pole distribution of material in the MW
satellite region in the 8.2 Gyr snapshot of our por simulation, shown in
Galactic coordinates. The title shows the total mass in the satellite region
and the fraction of this in stars. The colour of each pixel indicates the mass
in a square of $\left(l,b\right)$ with sides of length $6^{\circ}$. The pink
star shows the observed orbital pole of the Galactic SP (Table 9). The solid
pink curve shows a cone around this direction with opening angle of
$30^{\circ}$. The pink dot shows the preferred orbital pole in the simulation
(Equation 19), with the dot-dashed pink curve showing the estimated orbital
pole dispersion (Equation 20). The lack of material at very low and high
Galactic latitudes is caused by excluding material within $z_{\text{max}}=50$
kpc of the disc plane. Results appear similar for a different $z_{\text{max}}$
(Appendix D). _Right_ : The red bars show the cosine distribution of the angle
between the orbital pole $\widehat{\bm{h}}$ and
$\widehat{\bm{h}}_{\text{tot}}$, the preferred orbital pole of the simulated
SP (Equation 19). The open black bars show the orbital pole distribution of
classical MW satellites at Galactocentric distances of $50-250$ kpc based on
the positions and orbital poles listed in table 2 of Pawlowski & Kroupa
(2020). The angles shown here are relative to the observationally preferred
direction listed in Table 9.

Figure 8: Similar to Figure 7, but showing results for M31. Its disc spin vector is shown as a black star. The lack of material with an orbital pole close to this direction or its opposite is caused by the exclusion of material close to the M31 disc plane. Notice the smaller simulated disc-SP misalignment compared to the MW, which is likely due to M31 having a higher mass and a shorter scale length for the outer disc (Table 2). The simulated and observed orbital pole dispersions are also much smaller than for the MW. Results appear similar for a different $z_{\text{max}}$ (Appendix D). Table 9: The observed orbital pole of the MW SP (section 3 of Pawlowski & Kroupa, 2013) and of the M31 SP (section 4 of Pawlowski et al., 2013) in Galactic coordinates. We also show the angle between the spin vector of each SP and that of its parent disc, and the angle between the SP spin vectors. The MW-M31 orbital pole (shown in Figure 6) aligns fairly closely with that of the simulated M31 SP (shown in Figure 8). Galaxy | MW | M31
---|---|---
SP spin vector | $\left(176.4^{\circ},\,-15.0^{\circ}\right)$ | $\left(206.2^{\circ},\,7.8^{\circ}\right)$
Disc-SP misalignment | $75^{\circ}$ | $47^{\circ}$
Angle between SP spins | $37^{\circ}$

We find the angular momentum of all particles and gas cells in the satellite
region of each galaxy relative to the galaxy’s barycentre. The important
quantity for us is the direction of this angular momentum, which we use to
build up a distribution in Galactic coordinates. The procedure is the same as
that described in section 4.1 of BRZ18. The resulting orbital pole
distribution is shown in Figure 7 for the MW and in Figure 8 for M31. The
observed orbital poles of the MW and M31 SPs (summarized in Table 9) are
obtained from section 3 of Pawlowski & Kroupa (2013) and section 4 of
Pawlowski et al. (2013), respectively. These are shown as pink stars on the
left-hand panels, which our por model matches fairly well. In Appendix D, we
show that the appearance remains quite similar if we alter $z_{\text{max}}$ to
40 kpc or 60 kpc.

An important aspect of our analysis is an estimate for the dispersion in
orbital pole directions. On the observational side, we use $30^{\circ}$ for
the MW based on section 4 of Pawlowski & Kroupa (2013). For M31, we note that
its SP has an aspect ratio $3\times$ smaller than for the MW (table 3 of
Pawlowski et al., 2013). We therefore adopt an orbital pole dispersion of
$10^{\circ}$ for M31. These dispersions are illustrated by drawing a cone with
an opening angle equal to the estimated dispersion and an axis aligned with
the observed SP orbital pole direction listed in Table 9. These cones are
shown on the left-hand panels of Figures 7 and 8 using solid pink curves.

The dashed pink curves on these figures show analogous results for the
simulated tidal debris, whose preferred orbital pole
$\widehat{\bm{h}}_{\text{tot}}$ is shown with a pink dot in each case. We
estimate this using an iterative procedure where the initial guess is the
centre of the pixel in $\left(l,b\right)$ with the most mass. We then find

$\displaystyle\widehat{\bm{h}}_{\text{tot}}~{}\propto~{}\sum_{i}m_{i}\widehat{\bm{h}}_{i}\,,$
(19)

where each particle or gas cell in the satellite region has mass $m_{i}$ and
orbital pole direction $\widehat{\bm{h}}_{i}$ relative to its host galaxy. The
sum is taken over only those particles whose $\widehat{\bm{h}}_{i}$ aligns
with $\widehat{\bm{h}}_{\text{tot}}$ to better than $30^{\circ}$ for M31 or
$90^{\circ}$ for the MW, i.e. $3\times$ the observed orbital pole dispersion
in both cases. This restriction causes $\widehat{\bm{h}}_{\text{tot}}$ to
influence which particles and gas cells contribute to the sum, so we need to
repeat the process a few times until convergence is reached. We find that only
a handful of iterations are required to reach convergence in
$\widehat{\bm{h}}_{\text{tot}}$ to within machine precision. We then calculate
the orbital pole dispersion $\theta_{\text{rms}}$ using

$\displaystyle\theta_{\text{rms}}~{}=~{}\sqrt{\frac{\sum_{i}m_{i}{\theta_{i}}^{2}}{\sum_{i}m_{i}}}\,,$
(20)

where
$\theta_{i}\equiv\cos^{-1}\left(\widehat{\bm{h}}_{i}\cdot\widehat{\bm{h}}_{\text{tot}}\right)$
is the angle between $\widehat{\bm{h}}_{\text{tot}}$ and the orbital pole of
particle $i$. The sum is again taken over only those particles or gas cells in
the satellite region whose $\widehat{\bm{h}}_{i}$ is within the above-
mentioned cone around $\widehat{\bm{h}}_{\text{tot}}$. The so-obtained orbital
pole dispersion is $\theta_{\text{rms}}=49^{\circ}$ for the MW and
$19^{\circ}$ for M31, so our model naturally yields a lower
$\theta_{\text{rms}}$ around M31. The simulated dispersions are larger than
the observed ones, which we ascribe to the somewhat high 372 kK temperature
floor of our por simulations due to resolution limitations. It is also likely
that individual TDGs would form only in the densest regions, leading to a
narrower spread of orbital poles than for the tidal debris considered as a
whole.

The left-hand panels of Figures 7 and 8 reveal the expected gap in the orbital
pole distribution around the disc spin vector and the opposite direction
arising from our definition of the satellite region. We clarify this in the
M31 case by displaying its observed disc spin vector
$\left(238.65^{\circ},\,-26.89^{\circ}\right)$ as a black star in Figure 8
(section 2.1 of Banik & Zhao, 2018c). This is omitted for the MW because by
definition its disc spin points towards the south Galactic pole, leading to a
lack of material at very low and high Galactic latitudes in Figure 7.

The right-hand panels of these figures show the mass-weighted distribution of
$\cos\theta_{i}$, which should be uniform for a completely isotropic
distribution. In the MW case, we use open black bars to show the observed
distribution for classical satellites at Galactocentric distances of $50-250$
kpc (table 2 of Pawlowski & Kroupa, 2020). The result is similar to what we
obtain in our por simulation: both the model and observations show a
significant clustering of orbital poles, though observationally the clustering
is somewhat tighter. This could be related to the resolution and temperature
floor of our model. The orbital poles are much more clustered for the tidal
debris around M31 than for the MW, which as explained earlier is in line with
observations as the M31 SP is much thinner than that of the MW.

The titles of Figures 7 and 8 indicate that the satellite region of each
galaxy has a mass of $\approx 10^{9}M_{\odot}$. Since we expect the SP
material to mostly have been quite far out initially (section 5.2.4 of BRZ18,
), our estimated SP masses are quite sensitive to the reliability of
extrapolating the assumed exponential disc law to large radii. The relatively
small amount of material here means that altering the initial mass
distribution at large radii would hardly affect the gravitational field,
leaving the simulated SP orientations unchanged. Therefore, the SP masses are
not a strong test of our model. Bearing this in mind, we note that an SP mass
of $10^{9}M_{\odot}$ is reasonable for the MW because the LMC dominates the
baryonic mass in the Galactic SP. The RC of the LMC has a flatline level of
$\approx 70$ km/s (Alves & Nelson, 2000; van der Marel & Kallivayalil, 2014;
Vasiliev, 2018), which in a MOND context implies a mass of $1.5\times
10^{9}M_{\odot}$. While this is in reasonable agreement with our model,
another important aspect of the LMC is its rather high specific angular
momentum $h$ (Kallivayalil et al., 2013). This is quite difficult to reproduce
in our model, with only a small fraction of the tidal debris having a higher
$h$ than that of the LMC (Section 3.3.4). It is difficult to solve this
problem by postulating a much larger amount of tidal debris as then the
Galactic disc would need to be damaged much more significantly. Instead, the
problem may lie with resolution and the temperature floor: colder gas would
need more support from rotational motion, which should lead to higher $h$
initially. It is not presently clear whether this will solve the problem of
the LMC, so for the time being its properties remain somewhat problematic for
our flyby scenario.

In the M31 case, the simulated SP mass suggests that M32 is likely part of
this structure (see section 5.2.3 of BRZ18, ), and/or that all of the mass in
the SP has not condensed into individual satellites. This is quite possible
given that the orbital pole distribution is a non-uniform ring rather than a
single point (Figure 8). It could well be that only in the densest part of
this ring was the density high enough for the gas to undergo Jeans collapse
into TDGs. Moreover, gas accreted onto a newly formed TDG could be
subsequently expelled by feedback, with perhaps only a small fraction ending
up in stars.

To address these issues more thoroughly, a higher resolution simulation would
be required in which bound satellite galaxies form out of the tidal debris,
which is beyond the scope of this project. A previous high-resolution
hydrodynamical MOND simulation of interacting disc galaxies suggests that TDGs
should form out of the tidal debris (Renaud et al., 2016). Their work focused
on the Antennae (Mirabel et al., 1992), but our work indicates that a similar
process could well have played out in the much better observed LG. We hope
that the initial conditions of our best-fitting model (Table 4) serve as a
starting point for further work on the LG in MOND.

####  3.3.2 Radial distribution

Figure 9: The radial distribution of material in the MW (_top_) and M31
(_bottom_) satellite regions, showing contributions from stars (red) and gas
(blue). The observed satellite distribution is shown using a stepped solid
black line representing the number of satellites in four equally wide bins
covering $50-250$ kpc. These were obtained from the bold entries in table 4 of
Pawlowski et al. (2013) and scaled for clarity as indicated in the legend. Our
por simulation suggests that the SPs would initially have been dominated by
gas, though this need not be the case today as we do not simulate star
formation and feedback.

We use the top and bottom panel of Figure 9 to show the radial distribution of
material in the satellite region of the simulated MW and M31, respectively.
Each bar is divided into a red part indicating stars and a blue part
indicating gas. It is clear that the satellite regions are completely
dominated by gas. For the MW, this might be caused by its initial distribution
of gas being more extended than that of its stars (Table 2). However, the
stars and gas in M31 have the same initial surface density profile. The
dominance of gas in the satellite region could indicate that this is more
easily removed from the disc than stars due to ram pressure effects.

Our results can tentatively be compared with the observed distances of LG SP
members from their host galaxy. We obtain these from the bold entries in table
4 of Pawlowski et al. (2013), restricting further to only those satellites
within $50-250$ kpc of their parent galaxy. These observational results are
shown in Figure 9 using a solid black line. There is good overall agreement
with our model, but observationally there is a lack of observed satellites
$200-250$ kpc from their host. This could be due to the tidal debris at such
large distances being too sparse to form TDGs. Selection effects could also
play a role, especially around the MW where more distant satellites would be
fainter. While this effect is less significant around M31, the limited size of
the Pan-Andromeda Archaeological Survey (Ibata et al., 2014a; McConnachie et
al., 2018) means that observations of M31 do not extend to galactocentric
distances much beyond 150 kpc.

Our simulations do not allow star formation, so the high gas fraction in the
satellite regions merely indicates that the flyby created an initially gas-
rich distribution of tidal debris. In a more advanced simulation, the gas may
well condense into a small number of TDGs which then form stars (Renaud et
al., 2016). While the formation of TDGs and their internal evolution is beyond
the scope of this project, their overall distribution should be similar in a
higher resolution simulation that allows star formation. In this respect, it
is interesting that the radial distribution of simulated material in the
satellite region of each galaxy broadly agrees with that of its actual SP
members.

####  3.3.3 Orbital eccentricity

An important aspect of the LG satellite planes is their small dispersion of
orbital pole directions (Pawlowski, 2021). Another aspect to consider is the
motions of satellite galaxies within their preferred plane. To quantify this,
we determine the position $\bm{r}$ and velocity $\bm{v}$ of each particle or
gas cell relative to its parent galaxy (Section 3.2). We then find the unit
vectors $\widehat{\bm{r}}$ and $\widehat{\bm{v}}$, where
$\widehat{\bm{n}}\equiv\bm{n}/n$ for any vector $\bm{n}$. The distribution of
the angle between $\widehat{\bm{r}}$ and $\widehat{\bm{v}}$ provides a measure
of whether the velocity dispersion tensor is isotropic or tangentially biased,
and also whether there is any net radial motion that shows up as an asymmetry
between $\pm\widehat{\bm{r}}\cdot\widehat{\bm{v}}$.

We use Figure 10 to plot the distribution of
$\widehat{\bm{r}}\cdot\widehat{\bm{v}}$. The top panel shows the result for
the MW. If the $\widehat{\bm{v}}$ for each particle or gas cell were
distributed randomly relative to the outwards radial direction
$\widehat{\bm{r}}$, then $\widehat{\bm{r}}\cdot\widehat{\bm{v}}$ would be
distributed uniformly over the range $\left(-1,1\right)$. This is a reasonable
description of the figure, so there is not much tendency for the orbits to be
tangentially biased, somewhat at odds with observations (Cautun & Frenk, 2017;
Riley et al., 2019; Hammer et al., 2021). This is probably due to the high gas
temperature floor of 372 kK, which is required given numerical limitations. In
this sense, our result for the MW is similar to figure 8 of BRZ18. One
difference is that unlike in their work, our results show an asymmetry between
$\pm\widehat{\bm{r}}\cdot\widehat{\bm{v}}$ such that radial inflow is
preferred over radial outflow. This indicates that according to our
simulation, some tidal debris should should still be falling back towards the
MW. This could explain the properties of the Fornax dwarf spheroidal satellite
(Yang et al., 2022).

Figure 10: The cosine distribution of the angle between the position and
velocity of each particle relative to its host galaxy, shown for the MW
(_top_) and M31 (_bottom_) satellite regions. Contributions are shown from
stars (red) and gas (blue). If the velocities were oriented randomly relative
to the radial direction, then the distribution would be uniform. The asymmetry
between $\pm\widehat{\bm{r}}\cdot\widehat{\bm{v}}$ for the MW indicates that
material is still falling towards it. Notice the strong preference for
tangential motion in the M31 case.

The bottom panel of Figure 10 shows the
$\widehat{\bm{r}}\cdot\widehat{\bm{v}}$ distribution for M31. There is a sharp
peak near 0, indicating a strong preference for tangential motion. This is
quite different to the result obtained by BRZ18, whose figure 8 shows only a
mild preference for tangential motion for the tidal debris around M31. The
hydrodynamical nature of our simulations is probably the main reason for this
difference. Since a higher resolution simulation with a lower temperature
floor should have even more dissipation, our por simulation indicates that the
M31 SP members should have rather low orbital eccentricities. This is in line
with the observed PMs of the only two M31 SP members whose PMs are currently
known (Sohn et al., 2020; Pawlowski & Tony Sohn, 2021). There is not much
tendency for radial inflow as opposed to outflow: at low values of
$\left|\widehat{\bm{r}}\cdot\widehat{\bm{v}}\right|\la 0.5$, there is more
material at positive values, indicating a tendency for outflow. However, this
trend is reversed at high values of
$\left|\widehat{\bm{r}}\cdot\widehat{\bm{v}}\right|\ga 0.5$. Thus, the M31 SP
should be more nearly virialized than the MW SP. The strong tendency for
tangential motion in the M31 SP also suggests that it should be dynamically
colder than the MW SP, thereby having a lower aspect ratio and orbital pole
dispersion. This is indeed the case in our model (Equation 20) and
observationally (e.g., Pawlowski et al., 2013). The full 6D phase space
structure of the M31 SP is still unclear because this would require accurate
PMs for its relatively faint member satellites, which have only recently
become available in two cases (Sohn et al., 2020). While these are indicative
of corotation, a more detailed analysis will need to await further data
(Pawlowski & Tony Sohn, 2021).

####  3.3.4 Analogues to the LMC

The phase space structure of the Galactic SP is much better known because of
its proximity. This is clearly dominated by the LMC, so a realistic model
should have material on an orbit similar to the LMC. To investigate this, we
conduct a very similar analysis to section 5.4 of BRZ18. The LMC is fairly
close to pericentre at the moment (Kallivayalil et al., 2013), which could
well be something that even a quite realistic model does not reproduce simply
because the considered snapshot time is not when the LMC analogue is at
pericentre. Indeed, analogues to the Magellanic Clouds are not very common
observationally, so it may be quite difficult to reproduce certain aspects of
them even with the correct theory of nature (Robotham et al., 2012).
Therefore, we focus on quantities which should be rather more robust around
the whole orbit.

The high tangential velocity of the LMC for its position indicates a rather
high $h$ in the Galactocentric frame. To explain this, we should have material
in the satellite region with a similar $h$. The distribution of $h$ is shown
in Figure 11. The LMC value is based on table 5 of Kallivayalil et al. (2013),
which implies a Galactocentric tangential velocity of $314\pm 24$ km/s at a
distance of $r_{{}_{\text{LMC}}}=49.39$ kpc (Pietrzyński et al., 2013). More
recent observations indicate a similar LMC PM (Gaia Collaboration, 2018). It
is clear that the high $h$ of the LMC is somewhat challenging for our model
given its mass, though the simulated $h$ distribution for material in the MW
satellite region does extend beyond the LMC value.

Figure 11: The specific angular momentum distribution of material in the MW
satellite region, shown relative to that of the LMC. The rightmost bar
includes all material with higher values.

An important aspect of the LMC is its high orbital eccentricity. This is not
necessarily captured merely by looking at $h$, because $h$ would be rather
high even for material on a circular orbit if it is sufficiently far out. We
address this by following equation 51 of BRZ18 in defining an adjusted
specific energy

$\displaystyle\widetilde{E}~{}\equiv~{}\Phi\left(r\right)+\frac{1}{2}v^{2}-\frac{h^{2}}{2{r_{{}_{\text{LMC}}}}^{2}}\,,$
(21)

where the distance $r$ and velocity $v$ are Galactocentric. $\widetilde{E}$ is
the specific radial kinetic energy of a particle at $r=r_{{}_{\text{LMC}}}$,
with $\widetilde{E}<0$ indicating that the particle is not capable of reaching
this distance. It would have little difficulty doing so if $r\approx
r_{{}_{\text{LMC}}}$ in the simulation. But if $r\gg r_{{}_{\text{LMC}}}$,
then the particle would need to undergo substantial radial oscillations to
ever resemble the LMC at any point along its orbit. If it does so, then
ideally it should match the $\widetilde{E}$ value of the LMC, i.e. have an RV
of 64 km/s at $r=r_{{}_{\text{LMC}}}$. To calculate if this is the case, we
need to assume some form for the Galactic potential $\Phi\left(r\right)$. We
find this by treating the MW as an isolated point mass, neglecting M31 and the
EFE. This is justified because the relevant distances are much larger than the
extent of the MW disc (Bovy & Rix, 2013), but much smaller than the distance
to M31 (McConnachie, 2012) or the distance beyond which the EFE dominates
(Banik & Zhao, 2018b). The form of $\Phi$ follows from integrating the simple
interpolating function (Famaey & Binney, 2005) for a single point mass using a
hyperbolic substitution.

$\displaystyle\Phi\left(r\right)$ $\displaystyle=$
$\displaystyle\sqrt{GMa_{{}_{0}}}\left[\ln\left(1+\sqrt{1+\tilde{r}^{2}}\right)-\frac{1}{\tilde{r}}-\sqrt{\frac{1}{\tilde{r}^{2}}+1}\right]\,,$
$\displaystyle\tilde{r}$ $\displaystyle\equiv$
$\displaystyle\frac{2r}{r_{{}_{M}}}\,,\quad
r_{{}_{M}}\equiv\sqrt{\frac{GM}{a_{{}_{0}}}}\,.$ (22)

Here, $M=9.15\times 10^{10}\,M_{\odot}$ is the initial mass of the MW, while
$r_{{}_{M}}$ is its MOND radius, the distance beyond which MOND effects become
significant. Since ${r_{{}_{M}}\approx 10}$ kpc, it is not much smaller than
the LMC distance, so it is not accurate to assume that the problem is in the
deep-MOND limit.

Figure 12: Distribution of the adjusted energy (Equation 21) and specific
angular momentum relative to that of the LMC, shown for material in the MW
satellite region assuming a point mass Galactic potential (Equation 22). The
LMC value is shown using a white diamond, with the horizontal dashed line
through it showing a 24 km/s uncertainty in its Galactocentric tangential
velocity (see the text). We neglect uncertainty in its RV because uncertainty
in its PM is much larger. The solid blue curve shows values for particles on
circular orbits, with the orbital radius indicated in kpc as a text label next
to the corresponding filled blue dot.

The joint distribution of $\widetilde{E}$ and $h$ is shown in Figure 12, which
also has a blue curve showing values for particles on circular orbits at
different $r$ (indicated in kpc). Our results show that if considering both
parameters jointly, the observed LMC value is near the edge of the simulated
distribution. Its $\widetilde{E}$ is somewhat on the low side, while its $h$
is on the high side. Nonetheless, the observed combination of
$\left(\widetilde{E},h\right)$ for the LMC is still within the simulated
distribution for material in the MW satellite region. Therefore, the LMC does
not obviously present an insurmountable difficulty for our model. Given also
that material in the satellite region shows a sharp concentration of orbital
poles at about the right direction (Figure 7), our model provides a plausible
explanation for how one could get a large satellite with a rather high angular
momentum on a quite eccentric nearly polar orbit aligned with the SP.

Our model strongly suggests that the LMC formed as a TDG out of the MW-M31
flyby. Its high mass could indicate that it formed at the tip of the Galactic
tidal tail, where a massive TDG is more likely to form due to the weaker tidal
stress and EFE combined with the tendency for material to pile up at
apocentre. In this scenario, the LMC must have been orbiting the MW for
${\approx 7}$ Gyr. This is quite feasible in a MOND context as there is almost
no dynamical friction 40 kpc from a purely baryonic MW. Moreover, the time
dependence of the gravity from M31 could cause perigalacticon distances to
vary somewhat between orbits (see, e.g., figure 20 of Banik & Zhao, 2018d).
Assuming that the SMC also formed as a TDG on a similar orbit to the LMC as
part of the Galactic SP, one can imagine that the LMC and SMC subsequently
underwent one or more interactions, helping to explain the observed properties
of the Magellanic Stream (Hammer et al., 2015; Lucchini et al., 2020; Lucchini
et al., 2021). Indeed, it seems difficult to altogether avoid close
interactions between a large number of TDGs confined within a 2D structure by
virtue of how they formed.

The top panel of Figure 9 indicates that the Galactic SP would initially have
been mostly gas, though a small contribution from stars is also expected. This
might explain how the LMC came to have stars older than the flyby (Harris &
Zaritsky, 2009; Nidever et al., 2021). Our scenario would struggle to explain
a significant fraction of stars older than the flyby, though this depends on
the efficiency with which gas in the tidal arm condenses into stars in the
LMC. A lower efficiency would increase the relative importance of stars formed
before the flyby. The star formation histories of MW satellites are discussed
further in Section 4.1.

A common origin during a past MW-M31 flyby provides a natural explanation for
the alignment between the Galactic orbit of the LMC, the LMC-SMC orbit, and
the MW SP (figure 1 of Pawlowski, McGaugh & Jerjen, 2015). These planes also
align well with that defined by the Magellanic Stream (figure 1 of Pawlowski &
Kroupa, 2020). If instead the LMC and SMC were recently accreted by the MW due
to dynamical friction with its CDM halo (Besla et al., 2007), then these
geometric alignments would be fortuitous. This is because the orientation of
the MW SP would a priori be unrelated to the LMC-SMC binary orbit, which could
itself be oriented differently to their barycentre’s motion around the MW
(section 5.3.1 of Kroupa, 2015). One exception is if the dynamical influence
of the LMC is responsible for most Galactic satellites having orbital poles in
a very narrow range of directions, as proposed recently in the $\Lambda$CDM
context where the LMC should be more massive relative to the MW than in MOND
(Garavito-Camargo et al., 2021).999Conroy et al. (2021) recently obtained
tentative evidence for the perturbation induced by a massive LMC on the
Galactic stellar halo. However, their idealized simulation only works because
the test particles they consider have a very low specific angular momentum, so
even a very small tidal torque can significantly change the orbital pole. This
scenario is untenable for the classical MW satellites as they have quite high
specific angular momenta, which is related to their strong preference to be
moving tangentially rather than radially (Hammer et al., 2021). Consequently,
including the influence of the LMC in a backwards integration (Correa Magnus &
Vasiliev, 2022) still leaves the MW satellites with a clustered orbital pole
distribution prior to the LMC infall, implying that this scenario does not
explain the clustered orbital pole distribution of MW satellites in the
$\Lambda$CDM context (Pawlowski et al., 2021).

The LMC could also be related to the MW SP rather more directly if the LMC
brought in its own retinue of satellites (Samuel et al., 2021). However, their
claim to explain the MW SP in $\Lambda$CDM contradicts the fact that group
infall is already included in cosmological $\Lambda$CDM simulations.101010An
exception is if MW analogues very rarely have an LMC analogue around them, but
these rare cases are also those with a satellite system similar to that of the
MW. In this case, $\Lambda$CDM would not be consistent with the LMC. Based on
these, it was previously argued that “having most satellites accreted as a
single group or along a single filament is unlikely to explain the MW disc of
satellites” (Shao et al., 2018), an important reason being that dwarf galaxy
groups are typically much larger than the thickness of the LG SPs (figure 1 of
Metz et al., 2009b). It is difficult for an LMC-mass galaxy to bring in too
many satellites due to its lower mass than the MW, so the Galactocentric orbit
of the LMC should be aligned with a pre-existing chance alignment of Galactic
satellites into a planar structure. A recently fallen in group would also span
only a narrow range of angles on the sky (Santos-Santos et al., 2021), whereas
the MW SP members go around most of the sky (figure 2 of Pawlowski & Kroupa,
2020). Group infall is already included in cosmological $\Lambda$CDM
simulations like Illustris, so postulating it as a reasonably likely
explanation for the MW SP is not meaningful when the low likelihood of this
structure in a $\Lambda$CDM universe has already been demonstrated in a
cosmological context (Pawlowski & Kroupa, 2020), especially since high
resolution hydrodynamical simulations indicate that the impact of baryonic
physics is small (Samuel et al., 2021). Moreover, the velocities of MW
satellites are biased tangentially (Hammer et al., 2021), even though a radial
bias is expected from infall (Angus, Diaferio & Kroupa, 2011). Our model is
quite promising in this respect because it yields a significant tangential
bias to the velocities for the tidal debris around M31, though it does not
achieve this for the MW (Figure 10). However, a dissipative origin for the
Galactic SP remains a promising explanation for why its members have
tangentially biased orbits, especially when considering the numerical
limitations of our simulation.

A successful theory should explain not only the MW SP but also the SPs around
M31 (Ibata et al., 2013; Ibata et al., 2014b) and Cen A (Müller et al., 2018a;
Müller et al., 2021) without postulating a highly unusual chance alignment in
all cases. Rather, it should provide a mechanism for generating the
significant observed anisotropy in the only three systems where the 3D
distribution of satellites is well known. Ideally, this would be done with
minimal assumptions, e.g. based on the known existence of dissipative baryons
and the known formation of phase space-correlated TDGs (Mirabel et al., 1992;
Bournaud et al., 2004). In keeping with Occam’s Razor, it is also preferable
to use only one hypothetical galaxy interaction to explain both SPs in the LG.

###  3.4 The simulated MW and M31 discs today

Our results in Section 3.3 show that only a small fraction of the initial
$3.05\times 10^{11}\,M_{\odot}$ in the MW and M31 discs ends up in their
satellite regions. This indicates that the interaction is not very strong and
should still leave recognizable disc remnants, as also suggested by Figures 1
and 2. The reason is that the perigalacticon distance of 81.47 kpc is much
more than the disc scale length of either galaxy (Table 2), which also causes
there to be only a small amount of dynamical friction (Figure 3). This is
encouraging for our scenario because a strong encounter that transforms the MW
and M31 discs into ellipticals would not be realistic. In this section, we
take a closer look at the MW and M31 disc remnants (Sections 3.4.1 and 3.4.2,
respectively). In each case, our analyses focus on the region within 250 kpc
of the galaxy’s barycentre and within 50 kpc of its disc plane. Combining our
analyses of the disc and satellite regions therefore covers the entire 250 kpc
sphere around each galaxy.

We begin by finding the total angular momentum of all material within the disc
region of each galaxy relative to the galaxy’s barycentre. This gives the disc
spin vector of each galaxy in the analysed por simulation snapshot, which we
list in Table 10 in Galactic coordinates. This also shows the initial and
observed disc spin vector in each case, and the amount by which the disc has
precessed during the simulation. As discussed in section 4.2.5 of Banik et al.
(2020), this precession is partly due to the EFE, though we expect the flyby
to also induce disc precession. The total precession is $\approx 10^{\circ}$
in both cases, so differential precession of different parts of the MW/M31
disc should not warp it too much. Importantly, the simulated discs end up
oriented similarly to observations thanks to our iterative adjustments
(Section 2.3.1), especially in the M31 case where the disc-SP misalignment is
smaller (Figure 8). In what follows, we analyse each disc in a reference frame
rotated so the $z$-axis corresponds to the observed disc orientation as listed
in Table 10. We have checked explicitly that the results of our analyses
differ very little if we instead align the $z$-axis to the spin vector of each
simulated disc as given in Table 10. We avoid doing this for simplicity as it
would require iterative adjustments to the simulated spin vector of each disc.

Table 10: Spin vectors of the MW and M31 discs at the start of our por simulation and in its 8.2 Gyr snapshot, comparison of which yields the angle by which each disc has precessed. We also show the observed spin vector of each disc, which for M31 follows from its observed inclination and kinematic position angle (section 2.1 of Banik & Zhao, 2018c). Galaxy | MW | M31
---|---|---
Initial | $\left(55.11^{\circ},\,-84.89^{\circ}\right)$ | $\left(238.38^{\circ},\,-33.71^{\circ}\right)$
Final | $\left(231.05^{\circ},\,-85.63^{\circ}\right)$ | $\left(237.72^{\circ},\,-26.15^{\circ}\right)$
Observed | $\equiv\left(0^{\circ},\,-90^{\circ}\right)$ | $\left(238.65^{\circ},\,-26.89^{\circ}\right)$
Precession | $9.47^{\circ}$ | $7.59^{\circ}$
Mismatch with | $4.37^{\circ}$ | $1.11^{\circ}$
observations

####  3.4.1 MW stellar disc remnant

Figure 13: _Top left_ : The stellar surface density map of the MW as viewed
from along its observed disc spin vector, which is very similar to that in our
por simulation (Table 10). The positions are relative to the MW barycentre,
found as described in Section 3.2. The dot-dashed circles show radii of 2, 4,
and 8 kpc. _Top right_ : The MW in the cylindrical $rz$ projection (Banik et
al., 2020). _Bottom left_ : The stellar surface density profile of the MW (red
curve), found by binning the stellar particles in cylindrical polar radius.
The blue line shows an exponential profile whose scale length is an
observationally motivated 3.6 kpc (Jurić et al., 2008). _Bottom right_ : The
initial (dot-dashed blue) and final (solid red) stellar $\sigma_{z}$ of the MW
disc as a function of radius. The data point at the Solar circle radius of
8.178 kpc (Gravity Collaboration, 2019) shows $\sigma_{z}=40\pm 10$ km/s based
on results for different stellar samples in Yu & Liu (2018). Appendix E shows
a version of this panel for a higher resolution simulation.

The top left panel of Figure 13 shows the face-on view of the stars in the MW
8.2 Gyr after the start of our por simulation. This shows a regular
appearance. We also show the MW in the cylindrical $rz$ projection system
described in section 3.2 of Banik et al. (2020). The so-obtained $rz$ view of
the stars in the MW is shown in the top right panel. This confirms that a disc
remains despite the flyby and other simulated effects.

Our results allow us to obtain the surface density profile of the MW’s stellar
component. Since our simulation does not allow star formation, this should
correspond to something akin to the old Galactic thick disc (Gilmore & Reid,
1983), whose stars pre-date the flyby. By binning the results in cylindrical
polar radius, we obtain the stellar surface density profile
$\Sigma_{\star}\left(r\right)$ shown in the bottom left panel of Figure 13.
The red line shows the simulated result, while the blue line shows an
exponential law with a scale length of 3.6 kpc (Jurić et al., 2008). This is
consistent with Jayaraman et al. (2013), who found that the thick disc scale
length is $\approx 4$ kpc (see their section 4.4). It is also in line with Li
& Zhao (2017) and Mateu & Vivas (2018).

Another test of our model is the vertical velocity dispersion $\sigma_{z}$ at
different $r$, especially at the Solar circle. We find this by binning the
stars in $r$ and finding their mass-weighted $\sigma_{z}$ using an iterative
outlier rejection procedure with threshold of $3.29\sigma$, corresponding to
the 0.1% tail of a Gaussian distribution (the procedure is described further
in section 3.2.1 of Banik et al., 2020). Our result is shown in the bottom
right panel of Figure 13. At the Solar circle of 8.2 kpc (McMillan, 2017;
Gravity Collaboration, 2019), the simulated $\sigma_{z}$ rises from an initial
value of $\approx 10$ km/s to a present value of $\approx 50$ km/s. Note that
since our model does not allow star formation, this should be compared with
only those stars which existed at the time of the MW-M31 flyby. In this
regard, we mention that the old metal-poor sample at large height in e.g.
figure 4 of Yu & Liu (2018) has $\sigma_{z}\approx 40$ km/s, though results
differ by $\approx 10$ km/s depending on the exact sample used. We therefore
show a data point at $40\pm 10$ km/s at the Solar circle. In a higher
resolution simulation with all other parameters fixed, $\sigma_{z}\approx 40$
km/s at the Solar circle $-$ we present this in Appendix E. While results at
this level of detail may differ in a more advanced simulation, the broad
agreement in the disc scale length and $\sigma_{z}$ is encouraging.

####  3.4.2 M31 stellar disc remnant

Figure 14: Similar to Figure 13, but for M31. _Top left_ : Face-on view of
M31, showing a bar (unlike the simulated MW). _Top right_ : Cylindrical $rz$
projection of M31. _Bottom left_ : The stellar surface density profile of M31
(red curve). The blue line shows an exponential law with an observationally
motivated scale length of 5.3 kpc (Courteau et al., 2011). _Bottom right_ :
The initial (dot-dashed blue) and final (solid red) stellar $\sigma_{z}$ of
M31 as a function of radius, showing that it is expected to be dynamically
hotter than the MW (c.f. Appendix E).

The face-on view of M31 is shown in the top left panel of Figure 14, revealing
a bar. Observationally, M31 is indeed a barred spiral galaxy (Beaton et al.,
2007), with the bar having a length of $\approx 4$ kpc and a ratio of
corotation radius/bar length of $\mathcal{R}=1.6\pm 0.2$ (Blaña Díaz et al.,
2018). Detailed analysis of the bar goes beyond the scope of this
contribution, though bars have been considered before in MOND (Tiret & Combes,
2007, 2008; Combes, 2014). We refer the reader to Banik et al. (2020) for a
recent study of the M33 bar in MOND, and to Sellwood et al. (2019) for a
similar study in $\Lambda$CDM. A general study of galactic bars in CDM and
different theoretical frameworks is available in Roshan et al. (2021a), which
focused on galaxies with a central surface density similar to the MW and
M31.111111In an acceleration-dependent theory like MOND, there is a critical
surface density. Thus, galaxies with a similar central surface density should
behave similarly, albeit with dimensional quantities like lengths
appropriately scaled (see section 5.2 of Roshan et al., 2021a). Their work
showed that in isolated MOND simulations, we would typically expect such a
galaxy to have a bar with $\mathcal{R}\approx 1-1.2$ (see their figure 21), so
the bar of M31 is somewhat slower but consistent within uncertainties.

The $rz$ projection of M31 (top right panel of Figure 14) indicates that our
simulation retains a rather thin M31 disc. Its surface density profile is
shown in the bottom left panel, which also shows a blue line representing an
exponential law with an observationally motivated scale length of 5.3 kpc
(Courteau et al., 2011). This provides a fairly good description of our
simulated M31 disc remnant, especially when bearing in mind that their
estimated scale length has an uncertainty of 0.5 kpc.

Finally, the most interesting result about the M31 remnant is shown in the