even higher. As a result, the initially dynamically hotter disc in Roshan et
al. (2021a) remains dynamically hotter by the end, which is also apparent from
comparing the edge-on projections. The higher $\Sigma_{0}$ in their model
requires it to be dynamically hotter to remain Toomre stable, as can be
understood analytically (Banik et al., 2018b). We would thus expect that in
MOND, higher surface density galaxies would be less likely to have a
detectable bar due to being dynamically hotter. This is borne out by
comparison of the bar strength evolution between figure 13 of Roshan et al.
(2021a) and the lower $\Sigma_{0}$ model shown in our Figure 17.

At this stage, it is far from certain that MOND can solve the very significant
failure of $\Lambda$CDM with respect to the fraction of disc galaxies that
have a bar (Figure 18). Existing results suggest that MOND should give the
right declining trend with $M_{\star}$, even though the simulation of Banik et
al. (2020) was motivated by a desire to obtain a weak bar resembling M33,
while the work of Roshan et al. (2021a) considered only one model and was
targeted at the pattern speed of the bar rather than its strength. Though
these unlooked-for results are encouraging, further MOND simulations are
needed covering a wider range of galaxy properties. Eventually, cosmological
simulations will be necessary to handle issues like the EFE and interactions
between galaxies, and more generally to put galaxies into their proper
cosmological context (Section 9).

###  4.5 Bar pattern speed

Another crucial property of a galactic bar is its pattern speed $\Omega_{p}$,
the angular frequency at which the bar rotates almost as a solid body. Since
galaxies cover a wide range of mass, size, and rotation speed, they have a
huge range of characteristic angular frequency. It is therefore necessary to
normalize $\Omega_{p}$ in some way to arrive at a dimensionless quantity. For
this, we follow Elmegreen (1996) in using the parameter

$\displaystyle\mathcal{R}~{}\equiv~{}\frac{R_{c}}{R_{b}}\,,$ (68)

where $R_{b}$ is the semi-major axis of the bar and $R_{c}$ is its corotation
radius, the radius where $\Omega_{p}$ is also the angular velocity $v_{c}/r$
of a star on a circular orbit in the galactic potential. Bars are said to be
fast if $\mathcal{R}=1.0-1.4$ (Binney & Tremaine, 2008), while bars with
$\mathcal{R}>1.4$ are said to be slow as $\Omega_{p}<v_{c}/r$ at the end of
the bar. Theoretically, an ultrafast bar ($\mathcal{R}<1$) should be unstable,
so bars should not extend beyond their corotation radius (Contopoulos &
Papayannopoulos, 1980; Contopoulos, 1981; Contopoulos & Grosbol, 1989).

A bar rotating through a CDM halo is expected to experience dynamical friction
for much the same reason as any other massive moving object (Sellwood, 1980).
Their figure 11 shows that the disc rapidly starts losing angular momentum to
the halo as soon as the bar forms, though substantial bar slowdown was
prevented by this pioneering simulation having a limited duration of 1.2 Gyr.
Rapid slowdown of the bar was demonstrated in the work of Weinberg (1985) and
confirmed in later studies (Debattista & Sellwood, 1998, 2000; Holley-
Bockelmann et al., 2005; Klypin et al., 2009). However, not all such
simulations give slow bars (e.g. Athanassoula, 2003; Athanassoula et al.,
2013). For a review of isolated CDM-based models of bar evolution in disc
galaxies, we refer the reader to section 5.4.1 of Roshan et al. (2021a), which
also considers possible reasons for the different results. One of the most
crucial factors they identified was whether the halo parameters are “truly
what one expects in the $\Lambda$CDM paradigm”, or if they are unrealistic due
to inaccuracies like a small truncation radius for the CDM halo.

Figure 19: Time evolution of the bar pattern speed $\Omega_{p}$ in the
modified gravity theories MOND (black), MOG (green), and non-local gravity
(NLG; blue). Two Newtonian models are shown in which the CDM halo is either
included self-consistently (red) or as a fixed extra contribution to the
potential (orange), which can be thought of as the halo consisting of rigid
particles that provide but do not respond to gravity. This unphysical model
yields a nearly constant $\Omega_{p}$, as do the modified gravity theories.
The bar decelerates significantly in the realistic CDM model, as shown also in
the inset which evolves it for longer. Reproduced from figure 17 of Roshan et
al. (2021a).

These uncertainties can be reduced by considering a large galaxy population
drawn from a high-resolution cosmological hydrodynamical simulation of the
$\Lambda$CDM paradigm. Before turning to this, we discuss the results of
isolated galaxy simulations conducted by Roshan et al. (2021a) which give a
deeper insight into the problem. If $\Lambda$CDM is correct, the halo
parameters must reproduce the RAR (Figure 3), leaving little wiggle room.
Importantly, those authors also considered several other modified gravity
theories (summarized in their section 2). The evolution of $\Omega_{p}$ in
these models is shown in our Figure 19, which reproduces their figure 17. Two
models were considered in Newtonian gravity: a self-consistent live Plummer
halo (LPH) model and a model with a rigid Hernquist halo (RHH) that provides
no dynamical friction, achieved by treating the halo as providing a fixed
extra contribution to the potential. In both cases, the halo parameters were
chosen to reproduce the RAR in an exponential disc galaxy with
$\Sigma_{0}=10\,\Sigma_{M}$, intermediate between the MW and M31 values (see
their section 3.3). The unphysical RHH model illustrates how suppressing
dynamical friction from the halo allows the bar to maintain a nearly constant
$\Omega_{p}$, whereas the bar in the LPH model decelerates significantly, as
is particularly apparent in the extra 2 Gyr covered by the inset to Figure 19.
All considered modified gravity theories maintain a nearly constant
$\Omega_{p}$.

Figure 20: Posterior inference on the mean value of the $\mathcal{R}$
parameter (Equation 68) and its intrinsic dispersion in log-space based on
considering different snapshots of the simulations shown in Figure 19, with
the same colour used for the same model. The $1\sigma$ confidence region is
shown in all cases. The dotted red contour shows the extended evolution of the
CDM model in a live halo, whose shift to higher $\mathcal{R}$ is directly
linked with the bar deceleration evident in the inset to Figure 19. The dotted
green contour shows the observational estimate for the barred disc galaxy
population. While this cannot be directly compared to the above-mentioned
models, it can be compared to the galaxy population in the EAGLE100
cosmological hydrodynamical simulation of $\Lambda$CDM, results of which are
shown at redshift 0.5 (purple), 0.27 (yellow), and 0 (thick cyan) based on
$\mathcal{R}$ determinations by an independent team (Algorry et al., 2017).
Here also the mean $\mathcal{R}$ rises with time, and at $z=0$ is in
$7.96\sigma$ tension with local observations. Reproduced from figure 20 of
Roshan et al. (2021a).

These differences in the evolution of $\Omega_{p}$ translate directly into
differences in the distribution of $\mathcal{R}$. This is because the RCs
shown in figure 4 of Roshan et al. (2021a) rise almost linearly with radius in
the central region, so a small decrease in $\Omega_{p}$ translates into a
significant rise in $R_{c}$ and thus also in $\mathcal{R}$ (Equation 68). To
build up its statistics, those authors treated the same galaxy viewed at
different times as representative of a population of galaxies viewed at the
same time in the theory under consideration. The distribution of $\mathcal{R}$
was then analysed in log-space to yield a mean value and intrinsic dispersion.
For comparison, the study also considered different observational studies and
$\mathcal{R}$ parameter determinations for barred galaxies in the EAGLE100
simulation, as summarized in their table 2. We show the results in our Figure
20, which reproduces figure 20 of Roshan et al. (2021a). The modified gravity
theories yield $\mathcal{R}\approx 1$ with little intrinsic dispersion, quite
similar to the observations. In contrast, the decelerating bar in the LPH
model has a rather high $\mathcal{R}$, which rises further as the simulation
is evolved. This is because of angular momentum exchange between the bar and
halo, as demonstrated directly in figure 12 of Roshan et al. (2021a) by
considering how the angular momentum of each evolves over time. An important
aspect of these results is the distribution of $\mathcal{R}$ in EAGLE100
galaxies. This also shows a rising trend with time, and at $z=0$ is strongly
inconsistent with observations at $7.96\sigma$ confidence.

Given the high significance of this falsification, we briefly describe how
$\mathcal{R}$ is determined observationally. Available observational
determinations of the $\mathcal{R}$ parameter were compiled in the study of
Cuomo et al. (2020), which also contains a good review of how observers
determine the bar length and corotation radius. We therefore summarize only
the main points below.

Since the RC is generally well known in comparison to the other parameters,
obtaining $R_{c}$ is mainly a question of reliably determining $\Omega_{p}$.
This is easy from a movie of a galaxy, but remarkably is also possible using a
single snapshot with the model-independent Tremaine-Weinberg (TW) method
(Tremaine & Weinberg, 1984). The basic idea is to determine the so-called
kinematic and photometric integrals, defined as the luminosity-weighted
average LOS velocity $V_{\text{LOS}}$ and position $X$ parallel to the major
axis of the particles/stars, respectively.

$\displaystyle\langle V\rangle~{}$ $\displaystyle\equiv$
$\displaystyle~{}\frac{\int V_{\text{LOS}}\Sigma\,dX}{\int\Sigma\,dX}\,,$ (69)
$\displaystyle\langle X\rangle~{}$ $\displaystyle\equiv$
$\displaystyle~{}\frac{\int X\Sigma\,dX}{\int\Sigma\,dX}\,,$ (70)

with $\Sigma$ being the surface density of the tracer in arbitrary units. The
integrals have to be measured along apertures (or pseudo-slits) parallel to
the line of nodes between the disc and sky planes, between which the
inclination angle is $i$. The units of $\Sigma$ are not relevant, reducing
sensitivity to the assumed $M_{\star}/L$. The pattern speed is calculated from
the kinematic and photometric integrals as

$\displaystyle\Omega_{p}\sin i~{}=~{}\frac{\langle V\rangle}{\langle
X\rangle}\,.$ (71)

The bar length $R_{b}$ can be obtained in a few different ways from a
photograph. The main methods are isophotal ellipse fitting and Fourier
analysis. In the former, isopohotes are drawn joining points of constant
surface brightness, with the bar length defined as the semi-major axis of the
isophote having the highest ellipticity (Aguerri et al., 2009). In the latter,
the azimuthal surface brightness profile is decomposed into Fourier modes,
which are then used to estimate the ratio $I$ of surface brightness in the bar
direction and at right angles. $I$ usually rises rapidly to a maximum before
declining to a low value in the outskirts. Typically the idea is to find the
mean between the minimum and maximum $I$ and then find the corresponding
radius within the disc at which $I$ has this value, with the highest
determination used in case of multiple solutions (Aguerri et al., 2000).

To check if the EAGLE100-based falsification of $\Lambda$CDM persists with
more recent higher resolution simulations, Roshan et al. (2021b) analysed all
$z=0$ disc galaxies in EAGLE50, EAGLE100, TNG50, and TNG100 using the TW
method, thereby allowing a direct comparison with measurements of
$\mathcal{R}$ in real galaxies. Their figure 2 illustrates the steps involved
in determining $R_{b}$ and $\Omega_{p}$ for a single simulated galaxy with
identifier number 229935 in TNG50. The reliability of the TW method is
confirmed in their figure 3, which considers a simulated galaxy with known
$\Omega_{p}$ from time-resolved data. With $\Omega_{p}$ in hand, it is not
difficult to find $R_{c}$, as illustrated in their figure 4. On the
observational side, 42 galaxies were used which passed all the quality cuts to
be included in the final statistical analysis. The highest resolution
simulation considered in the study was TNG50, from which 209 galaxies were
usable. EAGLE100 had fewer usable galaxies, while TNG100 had more (see their
table 2). EAGLE50 had only 5 usable galaxies, so we do not discuss it further
given also that it has the same resolution as EAGLE100.

The final result of Roshan et al. (2021b) is their figure 8, reproduced here
as our Figure 21. The different $\Lambda$CDM simulations all give consistent
results, but the region they converge on is strongly excluded observationally.
For TNG50, the tension is $12.62\sigma$ based on 209 galaxies. Due to the
larger sample size of 745, TNG100 yields a slightly higher tension of
$13.56\sigma$, while the 70 EAGLE100 galaxies used in the final analysis yield
a slightly smaller tension of $9.69\sigma$. The authors also used a different
method to obtain the RC that should capture the maximum possible uncertainty
in it (see their section 3.10), but this only reduces the significance levels
by $\approx 1.5\sigma$ (see their table 5). Given also the similarity in bar
fraction between TNG100 and the higher resolution TNG50 (Figure 18), there is
compelling evidence that the latest highest-resolution cosmological
$\Lambda$CDM simulations have converged with respect to the $\mathcal{R}$
parameter distribution in barred disc galaxies, thereby falsifying
$\Lambda$CDM at overwhelming significance when confronted with the
observational compilation in Cuomo et al. (2020). For $\Lambda$CDM to be
consistent with the observed approximate equality between bar length and
corotation radius, the bar length would have to rise over time to keep up with
the rising corotation radius. However, observed bar lengths evolve little over
time (Kim et al., 2021a; Lee et al., 2022).

Figure 21: Posterior inference on the mean $\mathcal{R}$ and its intrinsic
dispersion for the low redshift galaxy population according to observations
(black; Cuomo et al., 2020) and the $\Lambda$CDM cosmological simulations
EAGLE100 (red), TNG100 (green), and TNG50 (blue). In all cases, the $1\sigma$,
$3\sigma$, and $5\sigma$ confidence regions are shown. Notice the consistency
between the different $\Lambda$CDM simulations, whose differently sized error
ellipses are correlated with differences in the sample size. Reproduced from
figure 8 of Roshan et al. (2021b).

This argument is statistical in nature, leaving plenty of room for _some_
galaxies to have fast bars in $\Lambda$CDM. Indeed, many such cases are
evident in figures 5 and 6 of Roshan et al. (2021b). It is therefore not
completely surprising that some studies found fast bars in the $\Lambda$CDM
paradigm (e.g. Fragkoudi et al., 2021). Those authors used 16 barred galaxies
from the 30 Auriga zoom-in galaxy simulations (Grand et al., 2017), finding a
mean $\mathcal{R}$ close to 1. Possible reasons for this result were discussed
in section 6 of Roshan et al. (2021b), whose main points we review here. Since
a much larger sample of galaxies is available from TNG50, the main advantage
of a zoom-in simulation would be if it reaches a higher resolution. In this
regard, Fragkoudi et al. (2021) checked if resolution effects were responsible
for their surprising result, but concluded in their section 3 that “resolution
is likely not the main culprit for the high $\mathcal{R}$ values found in
previous simulations, such as EAGLE and Illustris.” TNG50 is not mentioned
because their work pre-dates the publication of Roshan et al. (2021b), but the
resolution in TNG50 is almost as good as in the lower resolution simulations
of Fragkoudi et al. (2021), which reveal no major difference with the higher
resolution runs. Instead of resolution, the typically fast bars reported by
Fragkoudi et al. (2021) are almost certainly related to the sample selection.
Indeed, those authors mention that the galaxies they study are not on the
abundance matching relation expected in $\Lambda$CDM, even though the full
simulation volume usually is because the adjustable parameters are calibrated
partly to ensure this is so.

More generally, it is unclear how Fragkoudi et al. (2021) selected their
sample $-$ and whether this introduced biases. For example, requiring a strong
bar may mean that the galaxy necessarily had an interaction in the not too
distant past. If so, there would be less time for the bar to slow down via
dynamical friction against the halo, thus closing the gap with observations in
a fake way (Figure 20). Moreover, it is possible that strongly barred galaxies
make some aspects of the TW analysis less robust, with the analysis perhaps
giving an incorrectly long bar that leads to an underestimated $\mathcal{R}$
(Cuomo et al., 2021). Those authors showed how observed galaxies with
apparently ultrafast bars can largely be attributed to complications with the
TW analysis for galaxies with fast bars if they also have features like an
inner ring and/or strong spiral arms, artificially raising the estimated bar
length. It is easy to imagine that this also occurs in some simulated
galaxies, where it is known that bar-spiral arm alignment can lead to an
erroneously low measurement of $\mathcal{R}$ (Hilmi et al., 2020). Selecting
strongly barred galaxies might have inadvertently selected galaxies where just
such issues arise, perhaps because a bar aligned with the spiral arm would
also lead to an overestimated bar strength, raising the likelihood that the
galaxy passes the criteria in appendix A of Fragkoudi et al. (2021). It is
also possible that the original Auriga sample was itself not representative of
the full EAGLE100 disc galaxy population, whose $\mathcal{R}$ parameter
distribution is in $9.69\sigma$ tension with observations (Figure 21). Other
zoom-in simulations of $\Lambda$CDM galaxies also yield a significant fraction
of slow bars (e.g. Zana et al., 2018, 2019), which seems inevitable given the
predicted dynamical friction against the halo. A resimulation of a typical
TNG50 disc galaxy with a slow bar might help to clarify this, though the
resolution is already sufficient according to Fragkoudi et al. (2021).
Removing DM from the central region through enhanced feedback might reduce
dynamical friction on the bar, but this would make it difficult to satisfy the
RAR, especially for lower mass galaxies. Indeed, their study does not go down
to as low a mass as the observational compilation of Cuomo et al. (2020).

The failure of $\Lambda$CDM with respect to bar pattern speeds could be
rectified by assuming that the DM consists of low-mass particles like axions
with a long de Broglie wavelength, suppressing dynamical friction on kpc
scales (for a review, see Lee, 2018). However, a significant reduction in the
density fluctuations on this scale would make it difficult to form dwarf
galaxies in sufficient numbers (see section 6 of Roshan et al., 2021b). It is
also necessary to solve other serious problems with $\Lambda$CDM discussed
elsewhere in this review, many of which are on much larger scales where the
above models would lead to much the same behaviour as $\Lambda$CDM (see in
particular Sections 7 and 8). We therefore argue that it is untenable to solve
the many problems faced by $\Lambda$CDM merely by making a small adjustment to
the properties of DM on kpc scales. Indeed, this is unlikely to be a viable
solution even if we only consider evidence on this scale (Rogers & Peiris,
2021).

In summary, $\Lambda$CDM appears to get incorrect the fraction of galaxies
that are thin discs (Figure 14), the fraction of these discs that have bars
(Figure 18), and the pattern speeds of the bars which do form (Figure 21). The
number of spiral arms is also less than predicted, indicative of an incorrect
dispersion relation for perturbations to the disc (Section 4.2). Since the RC
and velocity dispersion are measurable, the issue is almost certainly with the
self-gravity term. Therefore, $\Lambda$CDM does not correctly address
gravitational dynamics in galaxies. Meanwhile, MOND provides a good physical
explanation for the general trends apparent in observations, with our
intuitive understanding backed up by several pioneering $N$-body and
hydrodynamical simulations. Further work is required to address the behaviour
of Milgromian disc galaxies in a cosmological context.

##  5 Interacting galaxies and satellite planes

We can learn a great deal about the properties of particles by colliding them
at high energy. Similarly deep insights might be gained by studying
interactions between galaxies. In this section, we consider what such studies
reveal about the likely composition of galaxies and the gravity law that
governs them. As observations improve, an important test using interacting
galaxies may become possible (Section 11.2).

###  5.1 Tidal stability

One of the most common types of galaxy interaction is when a primary galaxy
with mass $M_{p}$ tidally disturbs an orbiting dwarf companion with mass
$M_{c}\ll M_{p}$. Tidal stability considerations imply a maximum companion
size or tidal radius $r_{t}$ at fixed $M_{c}$ and distance $R$ from the host.
The basic principle is that the dwarf self-gravity should be comparable to the
difference in host gravity across the spatial extent of the dwarf. If the host
exerts gravity $g_{{}_{e}}$ on the dwarf, then the approximate Newtonian
condition for the maximum allowed dwarf size $r_{t}$ is

$\displaystyle\overbrace{\frac{GM_{c}}{{r_{t}}^{2}}}^{g_{{}_{N,i}}}~{}\approx~{}g_{{}_{e}}^{\prime}r_{t}\,,$
(72)

where ′ is used to denote a radial derivative. Notice that the host properties
and the host-centric distance are relevant only insofar as they affect the
tidal stress $g_{{}_{e}}^{\prime}$ on the dwarf.

The basic idea is the same in MOND. Since the difference in host gravity
across the dwarf is comparable to its self-gravity $g_{i}$ at the tidal
radius, it is clear that the host gravity itself must be much stronger than
$g_{i}$. This means that at the tidal radius of a dwarf, its internal dynamics
are necessarily EFE-dominated ($g_{{}_{e}}\gg g_{{}_{i}}$). Consequently, we
can use Equation 41 to estimate $g_{i}$. Neglecting factors of order unity
such as the angular dependence, we have that

$\displaystyle
g_{{}_{i}}~{}\approx~{}\frac{GM_{c}a_{{}_{0}}\nu_{{}_{e}}}{{r_{t}}^{2}}\,,$
(73)

where $\nu_{{}_{e}}$ is the value of $\nu$ at the location of the dwarf due to
$\bm{g}_{{}_{e}}$ alone. For a point-like host and assuming the deep-MOND
limit,

$\displaystyle\nu_{{}_{e}}~{}=~{}\sqrt{\frac{a_{{}_{0}}}{g_{{}_{N,e}}}}~{}=~{}R\sqrt{\frac{a_{{}_{0}}}{GM_{p}}}\,.$
(74)

In this configuration, the tidal stress is
$g_{{}_{e}}^{\prime}=\sqrt{GM_{p}a_{{}_{0}}}/R^{2}$, so the tidal radius of
the dwarf is:

$\displaystyle r_{t}~{}\approx~{}R\left(\frac{M_{c}}{M_{p}}\right)^{1/3}\,.$
(75)

Including the factors of order unity that we have neglected so far, Zhao
(2005) derived the following more accurate estimate (see its equation 12):

$\displaystyle\frac{r_{t}}{R}~{}=~{}\overbrace{0.374}^{2^{1/6}/3}\left(\frac{M_{c}}{M_{p}}\right)^{1/3}\,.$
(76)

Zhao & Tian (2006) derived a generalization of this in their equation 27 under
the approximation that $g_{{}_{e}}$ follows a power-law in $R$. This leads to
additional corrections of order unity to account for the fact that in general,
the host gravity does not follow an inverse distance law because it is an
extended source not fully in the deep-MOND regime. It is possible to write
their result in terms of $g_{{}_{e}}$ and $g_{{}_{e}}^{\prime}$, the only
aspects of the host potential which affect the internal dynamics of the dwarf.
Notice that both $g_{{}_{e}}$ and $g_{{}_{e}}^{\prime}$ are relevant in MOND
due to the EFE, while only $g_{{}_{e}}^{\prime}$ is relevant in Newtonian
gravity because it satisfies the strong equivalence principle (see Equation
72).

In MOND, the tidal radius $r_{t}$ of a dwarf can be calculated without
reference to its internal dynamics or its actual radius $r_{h}$. For a virial
analysis (e.g. Equation 56) to be valid in any theory, it is necessary that
$r_{h}\ll r_{t}$. McGaugh & Wolf (2010) showed in their figure 6 that many
dwarf spheroidal satellites of the MW would not satisfy this condition in
MOND, which may well explain why they deviate from the BTFR (Equation 53), or
more generally from MOND expectations for objects in dynamical equilibrium.
Indeed, there is a fairly tight correlation between the extent to which the
dwarfs deviate from MOND equilibrium expectations and their sky-projected
ellipticity (see their figure 3). Such a correlation would have to be a
coincidence in the $\Lambda$CDM framework because the analysed dwarfs should
all have $r_{h}\la 0.1\,r_{t}$, so tidal effects would be negligible and the
elliptical images should have a different cause. It would be helpful to
distinguish these scenarios by identifying or ruling out clearer signatures of
ongoing tidal disruption.

An important consistency check on the MOND scenario is that none of the dwarfs
analysed by McGaugh & Wolf (2010) have $r_{h}\gg r_{t}$, as shown in their
figure 6. Interestingly, this shows that the distribution of $r_{h}/r_{t}$
marginally reaches 1 in MOND. This could be an indication that dwarf galaxies
formed around the MW with a range of sizes and Galactocentric radii, but
dwarfs which were tidally unstable are no longer detectable. The fact that the
distribution of $r_{h}/r_{t}$ has an upper limit of 1 in the MOND context
would need to be a complete coincidence in the $\Lambda$CDM context because
the upper limit lies at 0.1 instead, which would have to be understood in some
other way.

Figure 22: Data on dwarf galaxies in the Fornax cluster (Venhola et al.,
2018). The $x$-axis shows the sky-projected separation from the cluster
centre. The $y$-axis shows the ratio between the optical effective radius
$R_{e}$ and the maximum radius $r_{{}_{\text{max}}}$ that the dwarf could have
given its stellar mass to still remain detectable above the survey
sensitivity. The colour of each point indicates whether the dwarf appears
disturbed (red) or undisturbed (blue) in a manual visual classification
(Venhola et al., 2022). There is a clear lack of dwarfs above the diagonal
line marked as a tidal edge. Notice that dwarfs just below this edge are much
more likely to appear disturbed. LSB dwarfs are missing close to the cluster
centre, but these are clearly detectable in the survey as many such dwarfs are
evident further out. The predicted location of the tidal edge is expected to
differ between theories, so comparison with the observed location should be a
promising way to distinguish them. Credit: Elena Asencio.

Tides are expected to be important in cases other than the MW, for instance in
the Fornax galaxy cluster 20 Mpc away (Jordán et al., 2007). The dwarfs in
this cluster have been surveyed quite deeply as part of the Fornax Deep Survey
(FDS; Venhola et al., 2018). When the dwarfs are plotted in the space of
surface brightness and projected separation from the cluster centre, a clear
tidal edge is apparent (Figure 22). The location of this tidal edge could be a
good way to understand how much self-gravity the dwarfs have, which is
expected to differ substantially depending on the gravity law and whether
dwarf galaxies have their own CDM halos. The FDS could thus provide a
promising way to test $\Lambda$CDM and MOND based on their ability to
reproduce this feature and other aspects of the Fornax dwarf population,
especially given that a recent study found strong evidence tor tidal
interactions playing a role (see section 7.4 of Venhola et al., 2022).
Preliminary results suggest a significant problem for $\Lambda$CDM because the
dwarfs should be immune to tides, but this issue seems to be resolved in MOND
where they would be much more easily disrupted (see
https://darkmattercrisis.wordpress.com/2021/08/02/, accessed on 28/05/2022). A
similar analysis might also be possible in the Coma cluster (Chilingarian et
al., 2019; Freundlich et al., 2022). In both cases, the results would be
sensitive to the additional physics required to reconcile MOND with galaxy
cluster observations, be this sterile neutrinos or something else (Section
7.1). While the discrepancy between the observed gravity and the MOND gravity
of the baryons alone is larger in Coma, it is also apparent in Fornax
(Samurović, 2016b).

The higher tidal susceptibility of dwarfs in MOND is due to the EFE and the
lack of protective CDM halos. As a result, dwarfs which should be suitable for
an equilibrium analysis in $\Lambda$CDM might not be amenable to this in MOND.
This fact was not fully appreciated in the MOND virial analysis of Fattahi et
al. (2018), leading to the erroneous conclusion that MOND faces a “possibly
insurmountable challenge” explaining the high observed $\sigma_{{}_{i}}$ in
some cases. While they considered the EFE, they did not consider whether tides
invalidated the assumption of equilibrium, even though tides were extensively
discussed in the $\Lambda$CDM part of their work. In principle, it is much
simpler to estimate the tidal susceptibility of a galaxy in MOND because the
purely baryonic nature of galaxies makes it much easier to estimate their
mass. Only dwarfs with a low tidal susceptibility in MOND are suitable for
analysis with Equation 56 or its generalization in Figure 5. More generally,
this can be inaccurate if the internal crossing time within a dwarf becomes
comparable to its orbital period around its host (Brada & Milgrom, 2000b;
McGaugh & Wolf, 2010). This is especially true for an eccentric orbit because
if the EFE changes sufficiently quickly, it is difficult for the dwarf to
adjust, leading to a memory effect lasting a few internal crossing times (Wu &
Kroupa, 2013). However, our results indicate that this is only a major issue
if the dwarf $r_{h}$ exceeds its MOND tidal radius (Zhao, 2005; Zhao & Tian,
2006). This is supported by the very small memory effect experienced by DF2 as
it orbits NGC 1052 in the simulation results shown in Figure 6.

###  5.2 Tidal streams and EFE-induced asymmetry

If a dwarf satellite galaxy has $r_{h}>r_{t}$, we expect it to undergo tidal
disruption. By this, we mean that stars would escape through the inner or
outer Lagrange point and enter an independent host-centric orbit little
affected by the disrupting dwarf. This would lead to a tidal stream, which
might be detectable for a long time $-$ possibly even outliving the satellite
itself. The shape and kinematics of the tidal stream can be used to set
constraints on the host potential, especially if a gravitationally bound
remnant is still recognisable and its kinematics can be measured. Due to the
faintness of tidal features and the importance of observing them in detail to
obtain possibly fundamental insights about gravity, we will focus on tidal
streams around the MW. However, it is sometimes possible to conduct meaningful
analyses of tidal streams around other galaxies (e.g. around M31; Fardal et
al., 2013).

One of the best-studied and most prominent Galactic tidal streams is that of
the Sagittarius dwarf satellite, which goes out far beyond the Galactic disc
(Ibata et al., 2001). Using $N$-body models in a $\Lambda$CDM context to
constrain the shape of the Galactic potential, those authors found that the MW
halo needs to be quite round, with an aspect ratio $>0.7$. This was a priori
not expected, as mentioned in their abstract. Another argument for a near-
spherical halo is the observed bifurcation of the Sagittarius tidal stream
(Fellhauer et al., 2006), which suggests that its orbital pole has precessed
very little. As orbital pole precession is generically expected in a non-
spherical halo, they argued for a nearly spherical halo.

Since there is no direct evidence for the DM halo of the MW, it is possible to
choose a near-spherical halo consistent with the Sagittarius tidal stream.
This would make for a similar halo to the phantom halo in MOND, which would be
almost spherical at distances exceeding a few disc scale lengths (Lüghausen et
al., 2015), though the EFE from other galaxies becomes relevant at distances
of several hundred kpc and causes the potential to become flattened in a well-
defined manner (as discussed in Section 2.4; see also Oria et al., 2021).
Consequently, it is difficult to distinguish between $\Lambda$CDM and MOND
using the Sagittarius tidal stream (Read & Moore, 2005). Even so, tidal
streams provide an important test of the much less flexible MOND scenario.
Moreover, other constraints can be used to check whether the DM halo
properties required to fit the Sagittarius tidal stream are reasonable in a
$\Lambda$CDM context.

The benchmark $\Lambda$CDM model of the Sagittarius tidal stream involves a
triaxial DM halo whose intermediate axis aligns with the disc (Law & Majewski,
2010). The failure to reproduce the observed stream bifurcation might indicate
that the Sagittarius dwarf was originally rotating (Peñarrubia et al., 2010).
Rotation has indeed been detected in its present-day remnant (del Pino et al.,
2021). Even so, the intermediate axis configuration proposed by Law & Majewski
(2010) appears to be unstable (Debattista et al., 2013) and is not consistent
with the tidal stream of the globular cluster Palomar 5 (Pearson et al.,
2015). They identified ‘stream fanning’ as a characteristic signature of a
triaxial halo, but the observed lack of stream fanning in Palomar 5 argues
against this. One reason is that the lower mass of Palomar 5 makes for a
colder tidal stream, allowing its morphology to be studied in more detail.
This bodes well for future surveys that should uncover additional faint tidal
streams. A reasonably good match can be obtained to the Sagittarius tidal
stream in CDM-based models by including a massive LMC (Vasiliev et al., 2021)
and allowing more flexibility in the halo density distribution. Their high
required LMC mass of $\left(1.3\pm 0.3\right)\times 10^{11}M_{\odot}$ is
consistent with the mass inferred from a $\Lambda$CDM timing argument analysis
of the LG (Peñarrubia et al., 2014, 2016; Banik & Zhao, 2017). These results
are consistent with a more thorough 3D timing argument analysis that gave the
LMC a mass of $2.03\times 10^{11}M_{\odot}$ (Banik & Zhao, 2018c), as
explained in footnote 1 to Banik et al. (2021a).

The first study to consider the Sagittarius tidal stream in MOND was able to
fit most of its observed properties quite well despite very little flexibility
(Thomas et al., 2017). For instance, the total luminosity of the stream
tightly constrains the initial progenitor mass, with the stellar mass-size
relation for field galaxies then constraining the initial size (Dabringhausen
& Kroupa, 2013). The only major deficiency identified by Thomas et al. (2017)
was with the RVs in part of the leading arm, an issue also faced by Newtonian
models (see their conclusions). In addition to the possibility that stars in
this sky region assigned to the stream are not truly part of it, the leading
arm properties are better reproduced if the MW is assumed to have a hot gas
corona flattened similarly to the MW satellite plane (Section 5.6). Though the
models used QUMOND, results should be similar in AQUAL (Banik & Zhao, 2018a).

The Sagittarius tidal stream is sufficiently distant that dynamical friction
on the progenitor would be negligible in the MOND scenario. This is not
necessarily true in $\Lambda$CDM, where a sufficiently massive progenitor
moving through a sufficiently dense halo would experience strong dynamical
friction, altering the path of the tidal stream and its kinematics. As a
result, observations place an upper limit on how much dynamical friction the
Sagittarius progenitor might have experienced, which in turn constrains a
combination of its mass and that of the MW halo (Dierickx & Loeb, 2017). Their
work highlights a path to a potentially clear falsification of MOND:
unambiguous evidence for a satellite experiencing dynamical friction despite
being outside the baryonic distribution of the host galaxy. Null detection
would set an upper limit to the dynamical friction, which may be more
stringent than how much dynamical friction is expected in $\Lambda$CDM. This
would favour the MOND scenario.

MOND has also been applied to the tidal stream of Palomar 5 (Thomas et al.,
2018), though without attempting to obtain a detailed fit. Those authors
explored the EFE in greater detail, finding that it might explain the observed
asymmetry between the leading and trailing arms. While the EFE-dominated point
mass potential (Equation 41) is symmetric with respect to
$\pm\bm{g}_{{}_{e}}$, the potential becomes asymmetric if the internal and
external gravity are comparable, as discussed in more detail in sections 4.4
and 4.5 of Banik et al. (2020). The asymmetry of the Palomar 5 tidal stream
might have other explanations like a perturbation by the Galactic bar (Bonaca
et al., 2020) or by a subhalo consisting almost entirely of CDM (Banik & Bovy,
2019). However, the latter would be a random process because the subhalo could
pass through the leading or trailing arm, so the asymmetry could well be in
the opposite sense to that expected in MOND. More such examples would be
required to identify if the agreement with MOND in the case of Palomar 5 is
purely coincidental or part of a more general trend.

###  5.3 Polar ring galaxies

Polar ring galaxies (Bertola & Galletta, 1978) consist of a disc galaxy with a
substantially misaligned ring of usually gaseous material (for a review, see
Combes, 2014a). These offer an interesting probe of gravity because RCs can be
determined for two separate planes in the same system. Dynamical models of
polar rings in Newtonian and Milgromian gravity reveal that MOND can naturally
explain the higher rotation velocity in the polar ring compared to the host
(Lüghausen et al., 2013). While this can be explained to some extent in the
Newtonian model, the MOND interpretation was favoured by Lüghausen et al.
(2014) due to its smaller degree of flexibility. As explained in their
conclusions, the higher rotation velocities in the polar ring arise from the
orbits here being close to circular, while orbits in the host galaxy are more
eccentric due to perturbations caused by the polar ring. The difference arises
mainly because the polar ring tends to be more extended than the host galaxy,
presumably to avoid destructive disc crossings. As a result, the host galaxy
can be considered a point mass at the centre of the polar ring, thus yielding
a nearly axisymmetric potential around the ring. However, if the polar ring
has a non-negligible mass compared to the host galaxy and is not too distant,
then the potential in the host galaxy cannot be axisymmetric with respect to
its rotation axis (Lüghausen et al., 2013). The higher induced eccentricity of
orbits in the host reduces its observationally inferred rotation velocity,
which theoretically can be calculated by considering a test particle on a
closed orbit moving through the numerically determined potential (see their
section 4).

###  5.4 Shell galaxies

The photographs of some galaxies reveal shells (Schweizer, 1980), which likely
arose when a dwarf galaxy encountered a much more massive host (Quinn, 1984).
The basic idea is that a dwarf galaxy on a nearly radial orbit gets disrupted
at pericentre. When the liberated stars reach apocentre, they spend a
relatively large amount of time in a narrow range of radii, leading to a much
higher surface brightness in a photograph (Hernquist & Quinn, 1987).

Shell galaxies test several aspects of galaxy physics. The dwarf is expected
to undergo dynamical friction, affecting the pericentre times. Some dynamical
friction is expected in MOND as the baryonic parts of the galaxies would
typically overlap at pericentre, but in the outer parts of the orbit,
dynamical friction would be significant only in the CDM scenario. In addition,
the galactic potential also affects the shell positions.

Hernquist & Quinn (1987) argued that the observed shells in NGC 3923 rule out
MOND, but their work suffered from several very serious problems that
completely invalidate their conclusion (Milgrom, 1988a). In addition to
observational difficulties and the possibility of missing shells, the model
may not be realistic at very small radii because the dwarf would in general
not be on a purely radial orbit and would generally start losing stars before
reaching pericentre. While these issues can be mitigated with deeper imaging
and excluding very small radii from the fits, a more complicated issue is that
the dwarf may not be completely destroyed at pericentre. As a result, each
shell can be assigned to one of many pericentric passages of the dwarf. The
particular passage chosen is known as the generation of the shell. A plausible
dynamical model should fit all the shell radii quite well using only a small
number of generations, assume only a small number of missing shells, and
perform poorly only at very small radii.

The shells of NGC 3923 have been analysed more recently in the MOND context
(Bílek et al., 2013). Using a semi-analytic method, those authors obtained a
very good fit to the positions of 25 shells out of the observed 27 using a
model with three generations. The maximum deviation in shell radius was only
5.4%, with the mismatch correlating well with the shell’s positional
uncertainty. The two shells that could not be matched well are at very low
radii, but these might be part of the fourth generation $-$ a possibility not
considered by the authors in order to limit the flexibility of their model.

The use of shells to test MOND was discussed further in Bílek et al. (2015a)
and Vakili et al. (2017), with the latter authors also considering MOG
(subsequently falsified in Figure 12). In general, the strong dynamical
friction expected in the $\Lambda$CDM scenario suggests that large faint
shells should be accompanied by small bright shells because the launch
velocity would drop substantially between pericentres (Vakili et al., 2017). A
more uniform distribution of shell radii would be expected in the absence of
CDM, which is more in line with observations of NGC 3923 (Bílek et al., 2016).
More information can be gained by taking spectra and determining the LOS
velocity distribution of the shells, which are expected to have a quadruple-
peaked profile (Bílek et al., 2015b).

Some of these techniques were recently applied to the shell galaxy NGC 474
(Bílek et al., 2022). Their $N$-body MOND simulation of it agrees well with
the shell positions and even with the velocity dispersion of the shell where a
comparison was made. In general, the authors “were able to find a scenario
that fits all available constraints well.” The flexibility of the model is
fairly small due to the lack of CDM and the use of only two shell generations
arising from two pericentre passages, a consequence of the first interaction
having been only 1.3 Gyr ago. The more flexible Newtonian model also yielded
an acceptable match, though it proved difficult to produce a long tidal stream
similar to the northern stream, possibly due to the stronger dynamical
friction at large separation in this model. Further exploration of the
parameter space may yield an improved match. Even if that were the case, we
consider the less flexible MOND approach to be favoured by its fairly good
agreement with the tests conducted so far in shell galaxies (Equation 79).

###  5.5 Tidal dwarf galaxies (TDGs)

When gas-rich galaxies interact, tidal forces can pull out a thin stream of
stars and gas known as a tidal tail (Toomre & Toomre, 1972). If this tail is
dense enough, part of it can undergo self-gravitating Jeans collapse into a
bound structure, which for some range of properties would be considered a TDG.
The formation of TDGs has long been known from an observational perspective,
e.g. in the interacting galaxies known as the Antennae (Mirabel et al., 1992)
and the Seashell (Bournaud et al., 2004). For a review of TDGs, we refer the
reader to Lelli et al. (2015).

If only a few TDGs are produced in each interaction and survive to the
present, then most dwarf galaxies might actually be TDGs (Okazaki & Taniguchi,
2000). Later studies suggested that TDGs probably comprise only a few percent
of dwarf galaxies at low redshift (Kaviraj et al., 2012). However, this is
only a lower limit because the tidal features that allow dwarfs to be
confirmed as TDGs are expected to fade away as the tidally expelled material
becomes more diffuse. Whether the TDG also gets tidally disrupted is less
clear, but the existence of confirmed TDGs up to 4 Gyr old suggests that they
can be long-lived (Duc et al., 2014). It can therefore be very difficult to
distinguish ancient TDGs from primordial dwarfs, especially if their
properties are similar because both are purely baryonic. As a result, the
actual frequency of TDGs is not reliably known. The strong correlation between
the interaction histories of nearby galaxies and how many satellites they have
could be an indication that TDGs are actually commonplace, with the
interaction history quantified from the bulge fraction (Kroupa et al., 2010;
López-Corredoira & Kroupa, 2016; Javanmardi et al., 2019) or from more
detailed stellar population studies (Smercina et al., 2022).

| Expected $\sigma_{{}_{i}}$ if gravity $\ldots$
---|---
Origin of dwarf | Newtonian | Milgromian
Primordial | High | High
Tidal | Low | High
Table 2: Expectations for the internal velocity dispersion of a dwarf galaxy
depending on how it formed and the gravity law. If Newtonian gravity applies,
the cosmological context is assumed to be the $\Lambda$CDM paradigm in which
primordial galaxies form inside CDM halos. TDGs would be more common in MOND
due to their stronger self-gravity making them harder to disrupt.

One possible distinguishing feature of TDGs regardless of the gravity law is
that they are expected to be relatively metal-rich for their mass (Duc et al.,
2014; Recchi et al., 2015; Haslbauer et al., 2019). Unfortunately, even this
signature would be undetectable for a really ancient TDG because the
progenitor disc galaxies would not have been very enriched at sufficiently
early times (Recchi et al., 2015). This raises the possibility that many
nearby dwarf galaxies are actually ancient TDGs (Okazaki & Taniguchi, 2000),
which in turn means that we might be able to measure the $\sigma_{{}_{i}}$ of
a nearby TDG. In the $\Lambda$CDM context, $\sigma_{{}_{i}}$ would be rather
low (Barnes & Hernquist, 1992; Elmegreen et al., 1993; Wetzstein et al., 2007;
Ploeckinger et al., 2015). This would cause a significant discrepancy with
MOND expectations as dwarfs generally have a low surface brightness that puts
them in the MOND regime (Equation 52). It is therefore possible to falsify
MOND by observing a virialized TDG if the $\Lambda$CDM scenario is correct,
even if the tidal origin of the dwarf is unclear. But if MOND is correct, the
$\sigma_{{}_{i}}$ of a TDG would be typical of a primordial dwarf, so a TDG
governed by MOND would generally be misinterpreted as a primordial dwarf in
the $\Lambda$CDM context. Table 2 clarifies the $\sigma_{{}_{i}}$ expectations
of both paradigms for dwarfs that formed in different ways.

Figure 23: The observed mass-radius relation of globular clusters and ultra-
compact dwarfs (green circles), ordinary and dwarf elliptical galaxies (red
diamonds), and TDG candidates (blue squares). Ram pressure-stripped TDGs (blue
triangles) were discussed further in Haslbauer et al. (2019). The black
crosses show TDGs in numerical CDM simulations (Wetzstein et al., 2007). The
solid black lines show tidal radii 100 kpc from a galaxy whose mass in
$M_{\odot}$ is (from top to bottom) $10^{10}$, $10^{11}$, and $10^{12}$.
Notice that observed TDGs do not define a distinct sequence in the mass-radius
plane, even though this is expected in $\Lambda$CDM (Haslbauer et al., 2019)
due to their lack of CDM, which makes their self-gravity much weaker for the
same baryonic content (see Table 2). Reproduced from figure 2 of Dabringhausen
& Kroupa (2013).

Unlike in MOND, the self-gravity of a $\Lambda$CDM dwarf galaxy differs
substantially depending on how it formed. We discuss later the possibility of
directly testing this through $\sigma_{{}_{i}}$ or $v_{c}$ measurements of a
virialized TDG, which represents a very challenging measurement (Section
11.2). The difference between the matter content of similarly luminous
primordial and tidal dwarfs in the $\Lambda$CDM context causes them to follow
distinct tracks in the mass-radius plane (Haslbauer et al., 2019). However,
the splitting evident in their figure 9 is not observed when we compare field
dwarfs with confirmed TDGs (Dabringhausen & Kroupa, 2013). This problem is
evident in their figure 2, reproduced here as our Figure 23. These results
strongly suggest that there is no fundamental distinction between dwarf
galaxies that formed primordially or out of tidal debris. Since TDGs should be
purely baryonic, the most logical conclusion is that all dwarfs are purely
baryonic, thus favouring the modified gravity solution to the missing gravity
problem on galaxy scales.

###  5.6 The Local Group satellite planes

Due to the lack of any fundamental difference between primordial and tidal
dwarfs in MOND, the critical issue becomes how to demonstrate that a dwarf
galaxy is a TDG. One silver lining is that as with a primordial dwarf, the
self-gravity of a TDG would be enhanced by the modification to gravity, making
TDGs much more robust to perturbations like tides and gas removal by stellar
feedback. TDGs are thus expected to be more common in MOND (as shown with
numerical simulations; Renaud et al., 2016). This raises the possibility that
a rather nearby dwarf galaxy is actually a TDG, allowing it to be observed in
sufficient detail to both measure $\sigma_{{}_{i}}$ and confirm its tidal
origin.

The nearest dwarfs are the satellites of the MW and M31, which have elevated
$\sigma_{{}_{i}}$ values that reveal large amounts of missing gravity. With
the MW satellites, it has been demonstrated that this cannot be explained by
tides in a Newtonian context (McGaugh & Wolf, 2010). It would moreover be very
unusual if we were observing so many dwarfs right at the moment of
dissolution, even if this might be a viable explanation in some cases (Kroupa,
1997; Klessen & Kroupa, 1998). The dwarf spheroidal satellites of M31 also
have rather high $\sigma_{{}_{i}}$ values consistent with MOND expectations
(McGaugh & Milgrom, 2013a, b). This is not obviously problematic for
$\Lambda$CDM if the LG satellites are mostly primordial dwarfs (Table 2). If
on the other hand they are mostly TDGs, this would rule out $\Lambda$CDM. We
therefore review whether the LG satellites are better understood as primordial
or as tidal dwarfs.

To try and distinguish these possibilities, we exploit the fact that TDGs
would typically form a phase space-correlated structure like a satellite plane
(Pawlowski et al., 2011; Haslbauer et al., 2019). Anisotropy of the satellite
system might therefore be one of the longest-lasting lines of evidence for the
TDG nature of a dwarf. Such evidence would necessarily be statistical in
nature as many individual dwarfs are needed to detect the anisotropy.
Nonetheless, the large number of LG dwarfs (e.g. McConnachie, 2012) suggests
that this could be a promising approach.

The classical satellite galaxies of the MW do indeed form a highly flattened
system orthogonal to the Galactic disc, as was already evident prior to the
formulation of MOND (Lynden-Bell, 1976, 1982). This was later confirmed from
the positions of subsequently discovered satellites (Kroupa et al., 2005).
Their larger sample size of 16 finally allowed the authors to conclude that
the anisotropic distribution of Galactic satellites strongly contradicts
$\Lambda$CDM expectations. Today, it is known that the ultra-faint satellites,
globular clusters, and tidal streams also align with the plane of classical
satellites (Pawlowski et al., 2012). This realization has led to the notion
that the orbital motion of material in this plane came about due to a past
encounter between the MW and another galaxy (Pawlowski et al., 2011).
Importantly, proper motion data confirm that most classical MW satellites
share a common orbital plane (Pawlowski & Kroupa, 2013, 2020) aligned with the
plane normal defined by the satellite positions alone (Santos-Santos et al.,
2020). Their velocities show a significant bias towards the tangential
direction, as occurs for a rotating disc (Cautun & Frenk, 2017; Riley et al.,
2019; Hammer et al., 2021). This is highly suggestive of a dissipational
origin, and was certainly not predicted a priori in the $\Lambda$CDM context
(Libeskind et al., 2009; Deason et al., 2011). The velocity data are quite
crucial to the argument because if the Galactic satellites are distributed
isotropically but observations cover only one great plane on the sky, the
distribution of velocities would still be random relative to this plane.

Such incomplete sky coverage effects can be minimized by considering instead
M31, whose satellite region subtends a much smaller region of the sky. M31 has
been homogeneously observed by the Pan-Andromeda Archaeological Survey
(McConnachie et al., 2009, 2018). This allowed Ibata et al. (2013) to confirm
earlier hints of a satellite plane around M31 (Metz et al., 2007).
Approximately half of the M31 satellites are part of its dominant satellite
plane, a much smaller fraction than for the MW. Since we observe the M31
satellite plane almost edge-on, a coherently rotating satellite plane should
have a gradient in RV across it, which is observed (Ibata et al., 2013).
Proper motions have recently become available for the plane members NGC 147
and NGC 185, demonstrating that their 3D velocities align with the plane
within uncertainties (Sohn et al., 2020). This significantly worsens the
problem for $\Lambda$CDM (Pawlowski & Tony Sohn, 2021). For a review of the LG
satellite planes and the recently discovered one around Centaurus A (Section
5.7), we refer the reader to Pawlowski (2018), who also considered whether
various proposed scenarios might or might not provide viable explanations for
these structures. Pawlowski (2021b) explores correlations between satellites
more broadly, in particular by also considering satellite pairs and asymmetry
of the satellite distribution. Some of the review is spent on discussing
possible explanations to the observations which do not sit well with
$\Lambda$CDM, while also highlighting the areas in good agreement. The
satellite planes sit alongside other problems for $\Lambda$CDM on galaxy
scales, many of which seem to be resolved in MOND (de Martino et al., 2020).
Their review also considers other modified gravity theories like MOG and
standard gravity solutions where the dark halos of galaxies have constituents
with different properties to CDM.

After carefully considering several proposed scenarios for how primordial CDM-
rich satellites might come to lie in a thin plane, Pawlowski et al. (2014)
concluded that none of them agree with observations for either the MW or M31.
Structures as anisotropic as their satellite planes are extremely rare in
cosmological simulations (Ibata et al., 2014b; Pawlowski & McGaugh, 2014b),
including hydrodynamical simulations (Ahmed et al., 2017; Shao et al., 2019;
Pawlowski & Kroupa, 2020) and simulations which approximately include the
effects of a central disc galaxy (Pawlowski et al., 2019). Previous arguments
against the group infall and filamentary accretion scenarios (Metz et al.,
2009; Pawlowski et al., 2014) were later independently confirmed (Shao et al.,
2018) using the EAGLE simulation (Crain et al., 2015; Schaye et al., 2015),
thereby including hydrodynamics in a $\Lambda$CDM cosmological context. If
group infall is somehow the solution, then the most plausible scenario is that
the LMC brought in a significant number of the classical MW satellites, in
which case they cannot be considered independently. Group infall is already
included in cosmological $\Lambda$CDM simulations, but it was nonetheless
suggested that the LMC could have brought in enough MW satellites to explain
the MW satellite plane in $\Lambda$CDM (Samuel et al., 2021). However, the LMC
could not plausibly have brought in enough satellites to explain the whole
Galactic satellite plane (Nichols et al., 2011). Indeed, a more recent study
showed that the LMC should have brought in “about 2 satellites with
$M_{\star}>10^{5}\,M_{\odot}$” (Santos-Santos et al., 2021). Given also that
the direct gravitational effect of the LMC is insufficient to induce a
sufficiently strong clustering of orbital poles (Garavito-Camargo et al.,
2021; Correa Magnus & Vasiliev, 2022; Pawlowski et al., 2022), some other
explanation is required for the 8/11 classical MW satellites that orbit within
a common plane (Pawlowski & Kroupa, 2020).

In summary, each LG satellite plane is in $3.55\sigma$ tension with
$\Lambda$CDM (table 3 of Banik et al., 2021a). The problems faced by
$\Lambda$CDM with these structures motivate us to reconsider their origin,
especially given past misunderstandings that if corrected lead to a much
bleaker assessment for standard cosmology (Pawlowski, 2021a). The anisotropy
of the LG satellite planes strongly suggests a tidal origin in both cases, but
the TDG scenario is not viable in $\Lambda$CDM because TDGs are already
included in hydrodynamical cosmological simulations like Illustris (Haslbauer
et al., 2019). Their very low frequency makes this an unlikely explanation.
Even if TDGs were more common, the high $\sigma_{{}_{i}}$ of the MW and M31
satellites rules out this possibility if gravity is Newtonian at low
accelerations (Table 2).

Both objections to the TDG scenario disappear in MOND. It therefore offers the
possibility that the highly flattened LG satellite planes formed as TDGs,
which we know from observations can form a flattened satellite system (Mirabel
et al., 1992; Bournaud et al., 2004). We have seen that it is extremely
difficult to come up with an alternative explanation. Moreover, MOND provides
a very natural answer to the question of which galactic interaction(s) formed
the LG satellite planes. This is because the stronger long-range gravity
between the MW and M31 acting on their nearly radial orbit (van der Marel et
al., 2012, 2019; Salomon et al., 2021) implies that they must have experienced
a close encounter $9\pm 2$ Gyr ago (Zhao et al., 2013). This is based on the
two-body force law in the deep-MOND limit (Equation 47). For the MW-M31
system, the MW mass fraction is close to 0.3, in which case the factor
$Q=0.7937$ (Banik et al., 2018a).

The MW and M31 are so widely separated that the EFE on the whole LG must be
considered. If the EFE dominates, the point mass potentials can be superposed
because the gravitational field is linear in the matter distribution (Equation
41). However, neither the isolated nor the EFE-dominated approximation is
valid because the LG is in the intermediate regime. This requires numerical
calculations of $g_{\text{rel}}$, the results of which can be fit fairly well
analytically (Banik et al., 2018a). In addition, tides from M81, IC 342, and
Centaurus A should also be considered, though fortunately the EFE must be
provided by more distant sources and can therefore be treated as acting on all
these objects in a similar manner. By taking these factors into account, Banik
et al. (2018a) showed that MOND is consistent with the classical LG timing
argument (Section 6.1) for a timing argument mass compatible with baryons
alone. Their work reaffirmed that a past MW-M31 flyby is required by MOND, or
else the timing argument mass would become extremely small and fall below the
directly observed baryonic mass within the LG (see also Benisty & Guendelman,
2020; McLeod & Lahav, 2020).

$N$-body simulations of the MW-M31 flyby in MOND produce a tidal tail between
the galaxies, suggesting that the LG satellite planes may have formed due to
the flyby (Bílek et al., 2018). This problem was also investigated using
restricted $N$-body models where the MW and M31 were treated as point masses
surrounded by test particle discs (Banik et al., 2018a). The tidal debris
around each galaxy generally had a preferred orbital pole. In some models,
this aligned with the observed satellite plane orbital pole for both the MW
and M31 (see their figure 5). This work has recently been followed up with
hydrodynamical MOND simulations of a flyby between two disc galaxies carefully
chosen to resemble the MW and M31 as they might have looked 9 Gyr ago (Banik
et al., 2022a). Their major result is the orbital pole distribution of the
material in the satellite region of each galaxy, i.e. within 250 kpc but
outside the disc (see their figures 7 and 8, the important aspects of which
are reproduced in our Figure 24). Taken in combination with several other
results from their paper which are not shown here, the model provides a good
description of the MW and M31 discs and satellite planes, including their
present-day orientations, disc scale lengths, and the MW-M31 separation.
Despite not being considered when selecting the best model, it is consistent
with the observed M31 proper motion (Salomon et al., 2021). Note that since
the models have a rather high temperature floor and do not allow star
formation, individual TDGs do not form, so the main aim was to compare the
preferred orbital pole of the tidal debris with the corresponding observations
(Pawlowski et al., 2013; Pawlowski & Kroupa, 2013). It seems very likely that
satellite planes of some sort would be produced from a past MW-M31 encounter,
especially as this could also have triggered the rapid formation of the
Galactic bulge (Ballero et al., 2007) and thick disc (Kilic et al., 2017),
which could be related to the formation and subsequent buckling of the
Galactic bar due to tidal torques (Grady et al., 2020). The flyby might also
explain the secondary peak in the age distribution of young halo globular
clusters in the Galactic halo (figure 9 of Mackey & Gilmore, 2004). In MOND,
TDGs would naturally have high $\sigma_{{}_{i}}$, just like any other galaxy
in the weak-field regime. Besides explaining their high observed
$\sigma_{{}_{\text{LOS}}}$, their enhanced self-gravity would make them more
robust to tides and feedback than TDGs in $\Lambda$CDM (Renaud et al., 2016).
On a larger scale, the flyby is a promising explanation for why the RVs of
galaxies in the NGC 3109 association exceed the $\Lambda$CDM prediction by
$\approx 100$ km/s (Section 6.1).

Figure 24: Results from a hydrodynamical MOND simulation of the MW and M31
discs interacting with each other 7.65 Gyr ago, showing the present-day
orbital pole distribution of material within 250 kpc but outside the disc of
the MW (left) and M31 (right). We show the mass in $6^{\circ}\times 6^{\circ}$
squares in Galactic coordinates, as indicated in the colour bar. The total
mass in the considered region is indicated at the top, along with the fraction
of this in stars. The MW disc spin vector is at the south Galactic pole by
convention (not shown). The preferred angular momentum of its satellite
galaxies is indicated with a magenta star in the left panel. The M31 disc spin
vector is shown with a black star in the right panel, while the preferred
angular momentum of the galaxies in its satellite plane is shown with a
magenta star. The lack of material whose angular momentum closely aligns or
anti-aligns with that of each disc is caused by the disc-shaped excluded
region having its symmetry axis aligned with that of the corresponding LG
galaxy. For both the MW and M31, notice the good agreement between the
preferred tidal debris orbital pole in the simulation and the preferred
orbital pole of its observed satellite galaxies, which were obtained from
section 3 of Pawlowski & Kroupa (2013) and section 4 of Pawlowski et al.
(2013), respectively. Adapted from figures 7 and 8 of Banik et al. (2022a) by
its lead author.

Regardless of the specific MOND scenario of a past MW-M31 flyby, the LG
satellite planes are strongly suggestive of a tidal origin. The high
$\sigma_{{}_{i}}$ of their members then requires a departure from GR. Given
the large body of prior work on the satellite plane issue (reviewed in
Pawlowski, 2018, 2021b), this reasoning is likely to withstand the test of
time. Moreover, any scenario for the LG satellite planes should not be too
contrived because satellite galaxy planes are not unique to the LG.

###  5.7 Satellite planes beyond the Local Group

One of the closest major galaxies beyond the LG is Centaurus A. Its satellites
were initially thought to define two planes (Tully et al., 2015), but later
observations clarified that both are part of one thicker plane which includes
Centaurus A (Müller et al., 2016). Subsequent results confirm the existence of
a thin kinematically coherent satellite plane around Centaurus A (Müller et
al., 2018a). Such a structure is expected in $\Lambda$CDM in only 0.2% of
cases (Müller et al., 2021). This work also showed the found analogues are
just chance alignments, presumably because TDGs are very rare in $\Lambda$CDM
(see figure 16 of Haslbauer et al., 2019). Therefore, the Centaurus A
satellite plane poses a significant challenge to $\Lambda$CDM despite the
present lack of $\sigma_{{}_{i}}$ measurements.

If this structure is interpreted as caused by a past interaction, then the
perturber may have already merged with Centaurus A. Alternatively, its
satellite plane might have formed due to a past interaction with M83, around
which there is some evidence for a satellite plane (Müller et al., 2018b).

Some hints for anisotropy in the satellite distribution are also evident
around NGC 2750 based on 7 satellites, which may form part of a lower mass
analogue to the satellite plane around Centaurus A (Paudel et al., 2021).
There is also evidence for a satellite plane around M81 (Chiboucas et al.,
2013) and NGC 253 (Martínez-Delgado et al., 2021; Mutlu-Pakdil et al., 2022).

It is difficult to constrain the 3D shape of the satellite distribution around
a more distant host galaxy. But if we are fortunate enough to view a satellite
plane close to edge-on, then satellites on opposite sides of the host galaxy
would have RVs with opposite signs once the host RV is subtracted. Such an
anti-correlation has been reported (Ibata et al., 2014a). It was later claimed
to not be a significant detection (Cautun et al., 2015), but the concerns
raised in this work were later addressed (Ibata et al., 2015). Those authors
used a large galaxy sample from SDSS to confirm the signal and its unexpected
nature in $\Lambda$CDM. While this remains a matter of conjecture, there is
certainly at least one satellite plane beyond the two in the LG, so these
structures should be reasonably common in a viable cosmological model.

##  6 Galaxy groups

MOND was originally motivated by the flat outer RCs of disc galaxies (Milgrom,
1983a). By definition, these regions must lie well into the MOND regime. Since
gravity declines with distance even in MOND, the gravity between galaxies
should lie deep in the MOND regime, as is the case for the MW and M31 (Banik
et al., 2018a). If the mutual acceleration between galaxies can be measured,
it might help to distinguish Newtonian from Milgromian gravity. For this, we
should primarily consider the phase space coordinates of galaxy barycentres,
thus treating them as point masses. This is especially helpful when studying
groups of galaxies, which can be analysed in detail if they are sufficiently
nearby to know the internal structure of the group. We describe some previous
work in this regard for the LG (Section 6.1) and other galaxy groups (Section
6.2). Further afield, binary galaxies offer a rather direct way to test
Equation 47 (see Section 6.3). It is also possible to analyse a distant galaxy
group statistically by applying the virial theorem, provided it has enough
members (Section 6.4).

###  6.1 The Local Group and the NGC 3109 association

One of the oldest yet still valid lines of evidence for the missing gravity
problem is the LG timing argument (Kahn & Woltjer, 1959). The basic idea is
that the gravity between the MW and M31 must explain how they turned around
from their presumed post-Big Bang recession as they are currently approaching
each other (e.g. van der Marel et al., 2012). We first discuss the timing
argument in $\Lambda$CDM to highlight certain problems, before explaining how
these might be resolved in MOND.

If we consider the MW and M31 as the only two point masses in an otherwise
homogeneous $\Lambda$CDM universe and assume a radial orbit, it is
straightforward to show that their total mass is $M_{\text{LG}}\approx 5\times
10^{12}M_{\odot}$ (e.g. Li & White, 2008). This estimate relies on the usual
assumption of no past close MW-M31 flyby because a past flyby would cause a
much higher $M_{\text{LG}}$ (Benisty et al., 2021) and lead to significant
dynamical friction during the encounter, inevitably causing a merger rather
than subsequent separation to large distance (e.g. Privon et al., 2013;
Kroupa, 2015). Since $M_{\text{LG}}$ and the unknown initial separation can
both be adjusted, we have two model parameters. These are sufficient to match
the two observational constraints, namely the Galactocentric distance and RV
of M31. Clearly, this is not a strong test of $\Lambda$CDM $-$ agreement with
observations is guaranteed in a model where galaxies can have arbitrarily
large amounts of invisible mass. Moreover, the mass ratio between the MW and
M31 is impossible to constrain based purely on the mutual gravity between
them.

Fortunately, we can consider other LG galaxies, including especially dwarf
galaxies that can be used as tracers of the velocity field. Such analyses have
a long history, with an early attempt by Sandage (1986) finding it difficult
to simultaneously explain all the data available at that time. Improvements in
survey capabilities led to an updated catalogue of LG dwarfs (McConnachie,
2012), which typically forms the basis of more recent timing argument
analyses. One such example is the study of Peñarrubia et al. (2014), which
found that matching observations requires an additional source of uncertainty
with magnitude ${35\pm 5}$ km/s. Both these analyses treat the LG potential as
spherically symmetric, thus assuming the tracer dwarf galaxies are much
further from the MW and M31 than their mutual separation of $783\pm 25$ kpc
(McConnachie, 2012). This is only approximately correct for dwarf galaxies
within 3 Mpc.

The first axisymmetric dynamical timing argument analysis of the LG was
conducted by Banik & Zhao (2016), who showed how the equations governing the
Newtonian timing argument can be derived from GR. By then, it was clear that
the MW-M31 orbit is indeed almost radial (van der Marel et al., 2012), which
was later confirmed (van der Marel et al., 2019; Salomon et al., 2021).
Centaurus A can also be included in an axisymmetric model of the LG because it
lies close to the MW-M31 line (Ma et al., 1998). This allowed Banik & Zhao
(2016) to consider LG dwarfs as test particles in the gravitational field of
these three massive moving objects. The computational costs were kept low,
allowing them to investigate a wide range of model parameters using a full
grid search. Despite this, a good fit was never obtained, even when
considering that the model is not a perfect representation of $\Lambda$CDM $-$
the expected model uncertainty is only 30 km/s based on the scatter about the
Hubble flow in detailed $N$-body simulations (Aragon-Calvo et al., 2011).
However, some LG dwarfs have a receding RV $>100$ km/s in excess of the best-
fitting model prediction (Banik & Zhao, 2016). While several LG dwarfs have
unexpectedly high RVs, there are very few cases where the RV is unexpectedly
low (see their figure 9). The typical mismatch between predicted and observed
RVs was estimated to be $45^{+7}_{-6}$ km/s.

To check whether this apparently serious problem for $\Lambda$CDM might be
alleviated in a more advanced 3D model of the LG, Banik & Zhao (2017) applied
the algorithm described in Shaya & Olling (2011). In the best-fitting model of
Banik & Zhao (2017), the typical RV mismatch was slightly higher than in the
2D case, with a clear asymmetry evident in their figures 7 and 9 such that an
unexpectedly high RV occurs much more often than an unexpectedly low RV. This
is despite including all the major galaxies in and around the LG out to almost
8 Mpc, as listed in their table 3. Though phrased more positively for
$\Lambda$CDM, similar results were obtained by Peebles (2017b) using a similar
algorithm. Banik & Zhao (2018c) borrowed this algorithm from Peebles and
improved it in certain important respects to give it the best possible chance
of finding dwarf galaxy trajectories compatible with the $\Lambda$CDM timing
argument (see their section 4.1). However, this had only a small impact on the
results.

The high RVs of some LG dwarfs outside the MW and M31 virial volumes can
partly be understood simply by reducing the LG timing argument mass. This is
of course not a solution to the problem as the masses of LG galaxies are
already treated as free parameters in timing argument analyses, but it is
helpful to bear this in mind when considering other works. In particular,
attempts to constrain $M_{\text{LG}}$ by considering all LG dwarfs can be
expected to yield a lower $M_{\text{LG}}$ than estimates considering only the
MW and M31. This could explain the unusually low $M_{\text{LG}}=\left(1.6\pm
0.2\right)\times 10^{12}M_{\odot}$ inferred in the former manner by Kashibadze
& Karachentsev (2018). However, Zhai et al. (2020) obtained a much higher
timing argument mass of $M_{\text{LG}}=4.4^{+2.4}_{-1.5}\times
10^{12}M_{\odot}$ by searching cosmological $\Lambda$CDM simulations for LG
analogues based on properties of the MW and M31 alone, especially with regards
to their relative separation and velocity. This estimate is in line with
earlier results and simple analytic estimates that consider only the MW and
M31 (Li & White, 2008). The mass of M31 alone has been estimated at
$1.9_{-0.4}^{+0.5}\times 10^{12}M_{\odot}$ based on its giant southern tidal
stream (Fardal et al., 2013). $M_{\text{LG}}$ must be higher as we also need
to include the MW and material outside the major LG galaxies. Thus, several
timing argument analyses of the whole LG have found that the RVs of some LG
dwarfs are too high to easily explain, though the tension is sometimes phrased
as an anomalously low $M_{\text{LG}}$ rather than as a problem for the
currently popular $\Lambda$CDM paradigm.

In addition to the timing argument, $\Lambda$CDM also faces a problem with the
abundance matching relation between the stellar and halo masses of nearby
galaxies. The prescription is that the observed number density of dwarf
galaxies per unit logarithmic interval in $M_{\star}$ tells us their halo mass
based on which DM halos in a cosmological simulation have the same number
density per unit logarithmic interval in the halo mass (Behroozi et al., 2013;
Moster et al., 2013). This simple idea encounters difficulties in the LG
because the Newtonian dynamical mass of the MW (Wang et al., 2022b) is lower
than its abundance matching mass given its $M_{\star}$ (McGaugh & van Dokkum,
2021). Those authors showed that the same problem arises in a more significant
way around M31, whose dynamical mass can be estimated from modelling its giant
southern tidal stream (Fardal et al., 2013). Lower values arise from
considering the M31 RC, worsening the problem (Chemin et al., 2009; Kafle et
al., 2018). Another way to look at the discrepancy is that given the combined
dynamical mass of the MW and M31, the stellar mass in the LG should be just
that of the MW, not the MW and M31. Reducing the LG mass to better fit the
kinematics of non-satellite LG dwarfs would worsen this problem.

The high-velocity galaxies (HVGs) identified by Banik & Zhao (2018c) are all
part of the NGC 3109 association, at least if we exclude HIZSS 3 due to its
very low Galactic latitude of $0.09^{\circ}$ and KKH 98 due to the RV mismatch
being only $65.5\pm 9.1$ km/s (see their table 3). This almost linear
association (Bellazzini et al., 2013) consists of several dwarf galaxies
moving away from the LG with rather high RVs. The NGC 3109 association was
studied in more detail by Pawlowski & McGaugh (2014a) in light of other LG
structures like the satellite planes (Section 5.6). Their work reached a
similar conclusion regarding the high RVs. Viewing the problem backwards in
time, the association should have been close to the MW in the past. This
suggests that an earlier gravitationally bound NGC 3109 galaxy group came
close to the MW and gained orbital energy at that time, which is most easily
understood as caused by a three-body interaction. Indeed, timing argument
calculations like that of Banik & Zhao (2018c) do not perfectly represent
$\Lambda$CDM because they lack dynamical friction and thus galaxy mergers,
which are however expected in $\Lambda$CDM. During such mergers, galaxies can
temporarily have high relative velocities. A dwarf near the spacetime location
of the merger could then gain orbital energy from the time-dependent LG
potential. This leads to ‘backsplash galaxies’ or backsplashers, defined as
objects on rather extreme orbits that were once within the virial radius of
their host but were subsequently carried outside of it. This backsplash
process was studied in detail by Sales et al. (2007), who found it very
difficult to get backsplashers at the $1.30\pm 0.02$ Mpc distance of NGC 3109
(Soszyński et al., 2006). This is also evident in figure 3 of Teyssier et al.
(2012), which additionally shows that backsplashers are very rarely more
massive than $10^{10}M_{\odot}$ regardless of their present position. In a
$\Lambda$CDM context, NGC 3109 has a mass of at least $4.0\times
10^{10}M_{\odot}$ (Li et al., 2020), so both its mass and its Galactocentric
distance are unusually high for a backsplasher. The distance from M31 is even
higher at almost 2 Mpc, so backsplash from M31 appears even less plausible.

Figure 25: The distribution of total mass and distance from the host for
backsplashers in the TNG300 $\Lambda$CDM cosmological simulation. The colour
is used to indicate the cosmic time when the backsplasher interacted with its
host, from which it may well be unbound today. Results are based on 13225 host
galaxies similar to the MW or M31, not all of which have an associated
backsplasher. The blue stars show different estimates for the mass of NGC 3109
in this framework. Its present Galactocentric distance is well known
(Soszyński et al., 2006; Dalcanton et al., 2009), with the $5\sigma$ lower
limit used here. The backsplashers are required to have not lost energy during
the interaction with their host, as this would make it even more difficult to
explain the observed high RV of NGC 3109 compared to $\Lambda$CDM expectations
(see the text). There are no backsplashers with a higher mass and distance
from the host than NGC 3109 is from the MW. The distance from M31 is $\approx
2.0$ Mpc (not shown). Reproduced from figure 8 of Banik et al. (2021a).

Early studies of the backsplash process considered only a small number of
analogues to the LG in $\Lambda$CDM simulations, e.g. Teyssier et al. (2012)
considered just one. Recently, it has become possible to revisit this issue
with the advent of high-resolution hydrodynamical cosmological $\Lambda$CDM
simulations like IllustrisTNG (Nelson et al., 2019). TNG300 is well suited for
this because its volume is large enough to capture many thousand galaxies
similar to the MW. Banik et al. (2021a) used TNG300 to quantify the
distribution of backsplashers in terms of present mass and distance from the
host within whose virial radius the backsplasher once was. The results are
shown in our Figure 25, which reproduces their figure 8. To account for the
high RV of NGC 3109, it is necessary to gain energy during a past interaction
with the host, so the backsplasher considered as a point mass is required to
leave the host with at least as much energy as it came in with (see their
section 5.3). The end result is that there are no backsplashers with similar
properties to NGC 3109 around any of the 13225 identified analogues to the MW
or M31, of which 640 are in 320 LG-like paired configurations. The results
remain similar when considering the DM-only version of TNG300, demonstrating
their robustness to how baryonic physics is handled. The issue with the
backsplash scenario is that dynamical friction between the CDM halo of NGC
3109 and the MW or M31 makes it likely that orbital energy would be lost
rather than gained during any close interaction. However, a more distant
interaction would be rather weak, so its effect should be correctly handled in
a few-body timing argument analysis of the LG. It is therefore very difficult
to understand the anomalously high RVs of galaxies in the NGC 3109 association
compared to $\Lambda$CDM expectations (Peebles, 2017b; Banik & Zhao, 2017,
2018c). The problem is evident even without doing a detailed timing argument
calculation and thereby dropping the requirement for NGC 3109 to be a
backsplasher: its observed RV is higher than the typical RVs of galaxies with
a similar mass and at a similar distance to a host resembling the MW or M31
($P=1.09\%$; see section 4.1 of Banik et al., 2021a). If the NGC 3109 analogue
is also required to be rather isolated (as is the case for the actual NGC 3109
group), then the tension worsens ($P=0.72\%$).

In the MOND context, the lack of CDM halos around galaxies removes the
objection regarding dynamical friction. Moreover, MOND requires a past close
MW-M31 flyby (Section 5.6), during which their high relative velocity would
have created substantial time variation of the LG potential. This would have
substantial effects on any dwarf near the spacetime location of the flyby.
Using a restricted $N$-body model of the LG in MOND, it was shown that some LG
dwarfs might well have been gravitationally slingshot outwards at high speed
(Banik & Zhao, 2018c). These dwarfs could certainly reach the distance of NGC
3109 (and even Leo P) relative to the LG barycentre, with the gravitational
slingshot causing a much higher present RV than galaxies at the same position
that never closely approached the MW or M31. Moreover, their figure 6 shows
that it would not be necessary for a dwarf galaxy that experienced such an
encounter to have passed very close to the MW or M31 in order to reach a large
present distance with a high RV, thus limiting damage to e.g. the disc of NGC
3109. Even so, it seems inevitable that a galaxy group would be tidally
disrupted, perhaps explaining the filamentary nature of the NGC 3109
association (Bellazzini et al., 2013) as an imprint of a past physical
association that is gravitationally unbound today (Kourkchi & Tully, 2017). It
could be a structure analogous to a tidal stream, but traced by galaxies
rather than stars. MOND therefore appears to provide a more holistic
explanation for the LG satellite planes and the larger scale velocity field
out to $\approx 2$ Mpc, all of which are very difficult to understand in
$\Lambda$CDM despite much effort by the community over many decades.

###  6.2 M81 and Hickson Compact Groups

One of the nearest galaxy groups is the M81 group, which lies at a distance of
only 3.6 Mpc (Gerke et al., 2011). Its properties are well explained by an
interaction between its three main members (M81, M82, and NGC 3077) in a
Newtonian context (Yun, 1999), with good agreement between the times of close
interactions and properties of the star formation histories (Rieke et al.,
1980). However, the models of Yun (1999) neglected dynamical friction, which
is certainly not valid in $\Lambda$CDM as including this process leads to the
rapid merging of all three galaxies (Thomson, Laine & Turnbull, 1999).

This problem has been revisited more recently, but the conclusions remain
similar (Oehm, Thies & Kroupa, 2017). The only possible scenario is if all
three galaxies recently encountered each other from large distances, which
seems quite unlikely $-$ yet needs to have occurred in one of the nearest
galaxy groups beyond the LG. Moreover, galaxies should have been receding
after the Big Bang, but the models require them to have been approaching each
other for several Gyr at nearly constant velocity (see their figure 14). This
implies a significant velocity perturbation at very early times, but cosmic
structures and peculiar velocities then were much less pronounced than today,
posing a similar problem to NGC 3109 and the LG (Section 6.1). In other words,
the timing argument is likely to be violated in a scenario where three
galaxies are currently rapidly approaching each other for the first time, as
would be required to allow for dynamical friction.

The existence of compact galaxy groups may pose problems for the $\Lambda$CDM
scenario in which galaxies are surrounded by extended CDM halos, since
dynamical friction between these halos would cause a rapid merger (Privon et
al., 2013; Kroupa, 2015). We would therefore need to observe the group
immediately before it merges, possibly posing a fine-tuning problem similar to
the M81 group. However, many so-called Hickson Compact Groups have been
identified (Hickson, 1982; Sohn et al., 2015, 2016). It has been argued that
their observed frequency can be explained in $\Lambda$CDM because some are
only chance alignments (Hartsuiker & Ploeckinger, 2020). This is surely not
always the case, and genuine compact groups should merge within a few Gyr.
However, those authors showed that it is possible for new groups to form out
of accretion along filaments.

This highlights the degeneracy between scenarios like $\Lambda$CDM with
frequent mergers that are only detectable as such for a brief period, and
scenarios where galaxy interactions are rare but would be detectable over most
of a Hubble time due to very weak dynamical friction (Renaud et al., 2016).
The observed frequency of interacting galaxies may be similar in both cases.
Further work is required to break this degeneracy. One possibility is that
frequent mergers between galaxies would lead to a low proportion of fragile
thin discs, so their high observed fraction argues against this scenario
(Section 4.1).

###  6.3 Binary galaxies

An isolated binary galaxy offers a geometrically simple configuration in which
the accelerations are typically deep into the MOND regime. Since we also know
that MOND was not originally formulated with binary galaxies in mind, these
can offer a strong test of the paradigm. The important MOND prediction is that
the mutual gravity $g_{\text{rel}}$ between two otherwise isolated point
masses in the deep-MOND regime is given by Equation 47. This can be translated
directly into a prediction for the relative velocity
$v_{c}=\sqrt{rg_{\text{rel}}}$ if the galaxies are on a circular orbit with
separation $r$. The predicted $v_{c}$ is independent of $r$. We expect the
actual distribution of the relative velocity $v_{\text{rel}}$ to be somewhat
broader than a $\delta$-function at $v_{c}$ due to orbital eccentricities and
inaccuracies in estimating $M_{s}$. However, these are rather modest effects
compared to the much larger variations possible in Newtonian gravity due to
scatter in the relation between baryonic and halo masses and the fact that
$v_{c}$ follows a Keplerian decline, so a more widely separated pair should
have a lower $v_{\text{rel}}$.

One complication is that observations from afar are sensitive to only the sky-
projected relative velocity. However, this can be accounted for statistically
because the relative LOS velocity should be uniformly distributed over the
range $0-v_{\text{rel}}$, allowing the distribution of $v_{\text{rel}}$ to be
recovered using a large sample (Nottale & Chamaraux, 2018b). Those authors
applied this technique to the HyperLEDA database (Makarov et al., 2014) to
obtain the Isolated Galaxy Pairs Catalog (IGPC; Nottale & Chamaraux, 2018a).
The deprojection algorithm is discussed in more detail in Nottale & Chamaraux
(2020). Their main result is that the $v_{\text{rel}}$ distribution shows a
clear peak at $\approx 150$ km/s. The reason for this was unclear $-$ their
section 5 goes off on a tangent about how exoplanets have a similar peak in
their distribution of Keplerian velocities, but this is completely irrelevant
to understanding galaxy dynamics. The section also states that the results
support the equivalence principle, even though this leads to GR and its
consequent prediction that $v_{\text{rel}}$ should undergo a Keplerian decline
with increasing separation, contrary to the observations.

To take advantage of subsequent additions to the HyperLEDA database, Scarpa et
al. (2022a) reanalysed it to identify galaxy pairs using a similar technique
to Nottale & Chamaraux (2020). This confirmed the existence of a peak in the
$v_{\text{rel}}$ distribution at $\approx 150$ km/s (see figure 1 of Scarpa et
al., 2022a). Those authors pointed out that such a clear peak is rather
unlikely in $\Lambda$CDM according to the prior CDM-only Millennium simulation
(Moreno et al., 2013), most likely due to the reasons discussed above.
However, the peak is in line with MOND expectations due to the rather narrow
range in $M_{s}$ and the fact that the predicted
$v_{\text{rel}}\propto\sqrt[{}^{4}]{M_{s}}$ (see figure 7 of Scarpa et al.,
2022a). An updated version of this analysis has recently been published
(Scarpa et al., 2022b). Its figure 1 (reproduced here as our Figure 26) shows
the distribution of $v_{\text{rel}}$, which has a clear peak at $\approx 130$
km/s for all three choices of width for the velocity bins used in the
reconstruction.

Figure 26: Deprojected distribution of $v_{\text{rel}}$ between 8571 isolated
galaxy pairs. Different colours show different velocity bin widths used in the
reconstruction. The dashed vertical line at 132 km/s marks the location of the
very clear peak in the distribution, which is at a similar location in all
three cases. Reproduced from figure 1 of Scarpa et al. (2022b).

To assess whether these results are consistent with MOND, it is necessary to
compare them with the circular velocities predicted by Equation 47. Even for
the case of purely circular orbits and a perfect reconstruction of the
$v_{\text{rel}}$ distribution, we still expect a spread in values because of
differences in the total mass and mass ratio between galaxy pairs. The
predicted distribution of $v_{c}$ is shown in Figure 27, which reproduces
figure 5 of Scarpa et al. (2022b). The black line shows a Gaussian fit to the
peak region in the observationally reconstructed $v_{\text{rel}}$
distribution. The red histogram shows the predicted distribution of $v_{c}$.
This provides a rather good match to the observations. There is an excess of
systems with a lower velocity than the predicted $v_{c}$, which could well be
due to the fact that masses on an eccentric orbit spend more time near
apocentre, where they are moving slower. In general, the clear peak in the
$v_{\text{rel}}$ distribution is well understood in MOND as a consequence of a
narrow range of total $M_{s}$ for each pair combined with only a weak
dependence on this and little sensitivity to the mass ratio (Equation 47).

Figure 27: The red histogram shows the distribution of the MOND-predicted
circular velocity (Equation 47) between the binary galaxies considered in
Figure 26. The black curve shows the Gaussian fit to the deprojected
$v_{\text{rel}}$ distribution shown in that figure after binning more finely
in the peak region. The best-fitting $M_{\star}/L_{B}$ is shown at the top
left. Reproduced from figure 5 of Scarpa et al. (2022b).

These results only consider a narrow range of galaxy luminosity, so the galaxy
pairs all have a similar total $M_{s}$. This makes it difficult to test the
predicted scaling behaviour $v_{c}\propto\sqrt[{}^{4}]{M_{s}}$. We can test
this using galaxy pairs in the GAMA survey (Driver et al., 2011) thanks to the
work of Loveday et al. (2018). Their abstract indicates that the typical
difference in the LOS peculiar velocity between two galaxies increases by a
factor of 3 for a $100\times$ increase in luminosity, which translates to
about the same increase in $M_{s}$. This would imply that the logarithmic
slope of the relation is close to $\ln 3/\ln 100=0.24$, which is almost
identical to the MOND prediction of $1/4$. In a $\Lambda$CDM context, there
appears to be a problem at the low-mass end due to the abundance matching halo
masses being too large, a problem which is also evident in the LG (Section
6.1).

###  6.4 Virial analysis of galaxy groups

In more distant galaxy groups, the data quality becomes poor such that only
the most basic analyses remain possible. Assuming a galaxy group is in
equilibrium, the simplest dynamical analysis is an application of the virial
theorem to deduce the group dynamical mass. This can be compared with the mass
estimated in other ways, including especially the photometric $M_{s}$.

Galaxy groups are generally in the deep-MOND limit because the gravity at e.g.
the Solar circle of the MW is comparable to $a_{{}_{0}}$ (Gaia Collaboration,
2021b), so the gravity on neighbouring galaxies is much smaller still.
Therefore, MOND analyses of galaxy groups generally make use of Equation 56,
which is valid for an isolated spherically symmetric system in the deep-MOND
limit. The equation assumes virial equilibrium and an isotropic velocity
dispersion tensor. These assumptions may not hold exactly, partly because the
longer dynamical timescales of galaxy groups compared to individual galaxies
can make the assumption of equilibrium less accurate. Even so, Equation 56
provides a good starting point to check whether MOND gives reasonable results
in systems where the departure from Newtonian gravity is expected to be large.

This equation was first applied to galaxy groups by Milgrom (1998), who found
a dynamical $M/L<10$ that might be explained using baryons alone. Milgrom
(2002) analysed 9 galaxy groups within 5 Mpc which exhibit often very large
amounts of missing gravity in a conventional context. The study found that the
MOND dynamical $M/L$ ratios were order unity. This remained true for groups
with and without “luminous” galaxies, i.e. dwarf-only groups and groups
containing more massive galaxies. However, these two types of group have
significantly different proportions of missing gravity (Tully et al., 2002),
which in a MOND context can be understood due to the dwarf-only groups having
lower internal accelerations.

This topic was revisited in a MOND analysis of 53 galaxy groups (Milgrom,
2018a), which was later extended to 56 medium-richness groups despite strict
quality cuts (Milgrom, 2019a). The MOND dynamical $M/L$ ratios were found to
be order unity in the near-infrared $K$ band within uncertainties, suggesting
reasonable agreement. These results should be considered in light of dwarf
galaxies also following Equation 56 (e.g. McGaugh et al., 2021). We provide a
holistic picture in our Figure 28, which reproduces figure 6 of Milgrom
(2019a). It is clear that Equation 56 provides a good description for the
dynamics of pressure-supported systems across 7 dex in $M_{\star}$, with the
low-mass end probed by stars in dwarf galaxies and the high-mass end by
galaxies in galaxy groups.

Figure 28: The $K$-band luminosity of galaxy groups in Solar units as a
function of their $\sigma_{{}_{\text{LOS}}}$ in km/s, with both shown on
logarithmic axes (solid black dots towards the upper right). The open circles
towards the lower left show dwarf galaxies, whose observed $V$-band
luminosities have been scaled by $\alpha=2/0.7$ to account for the higher
expected $M_{\star}/L$ in the $V$-band, the goal being to correlate
$\sigma_{{}_{\text{LOS}}}$ with mass. The solid blue line is the prediction of
Equation 56 for $M_{\star}/L_{K}=0.7$ (McGaugh & Schombert, 2014), with the
dashed (dotted) lines showing a factor of 2 (4) variation on either side.
Notice the huge dynamic range in luminosity shown here. Reproduced from figure
6 of Milgrom (2019a) by permission of Moti Milgrom and the American Physical
Society.

Milgrom (2019a) also attempted to detect the EFE in terms of groups with a
very small internal acceleration having a lower $\sigma_{{}_{\text{LOS}}}$
than implied by Equation 56, which is valid only for isolated systems.
However, no strong evidence was found for or against the EFE assuming a
typical external field strength of $0.03\,a_{{}_{0}}$ (e.g. Famaey et al.,
2007; Banik et al., 2018a; Haslbauer et al., 2020). One reason could be that
groups would experience a different EFE depending on their environment, so
assuming the same EFE on all groups is inaccurate. At the galaxy scale, the
EFE preferred by MOND RC fits depends on the environment of the galaxy (Chae
et al., 2020b, 2021), so it would be natural if the same applies to galaxy
groups. Checking this would require a large-scale map of the gravitational
field (e.g. Desmond et al., 2018).

While the 3D structure of a galaxy group is generally unknown, recent
observations have clarified the 3D structure of the Centaurus A galaxy group
and allowed a kinematic analysis (Müller et al., 2022). Allowing for both
rotation and random motion, those authors detected both at high significance
based on the RVs of 28 member galaxies and distances for 27 of them. Both
contribute similarly to supporting the system against gravity, whose strength
can thus be inferred from the kinematics. The estimated group $v_{c}$ of
$258\pm 57$ km/s is consistent with the BTFR (Equation 53). Since the internal
gravity should dominate over the EFE out to 500 kpc (Oria et al., 2021) and
most of the used tracers are much closer in, the kinematics of the Centaurus A
group are consistent with MOND expectations (Müller et al., 2022). This
success is on a much larger scale than Centaurus A itself, which is also
consistent with MOND (Samurović, 2016a). In a $\Lambda$CDM context, the large
size of this group and the low accelerations at its outskirts mean that a
significant proportion of DM is required, much more than the cosmic fraction
inferred from e.g. fits to the CMB power spectrum. As a result, there should
be $\approx 8\times 10^{11}M_{\odot}$ of hot gas in this very nearby galaxy
group (Müller et al., 2022). X-ray observations indicate that the actual
amount is $\approx 4\%$ of this out to the virial radius of 300 kpc (Gaspari
et al., 2019).

In addition to roughly spherical galaxy groups, it is also possible to test
MOND with filamentary structures like the Perseus-Pisces superfilament (Haynes
& Giovanelli, 1986). The observations were found to be broadly consistent with
MOND, though the uncertainties were large because it was necessary to make
various simplifying assumptions (Milgrom, 1997). Nonetheless, the study is
interesting because it probed very low accelerations (only a few percent of
$a_{{}_{0}}$), which in a Newtonian analysis leads to a very large proportion
of missing gravity (Eisenstein et al., 1997). Moreover, the linear mass
distribution is rather unique $-$ typical RC analyses only probe a disc
geometry, and even then the critical outer parts of the RC are mainly
sensitive to the total baryonic mass. While the inner parts of the RC are
somewhat more sensitive to the disc geometry, such analyses mainly test the
force law for an isolated point mass (Equation 9), especially if the focus is
on the outer regions (e.g. Lelli et al., 2019).

In the $\Lambda$CDM paradigm, the above-mentioned observational results merely
indicate how much DM each system should have. It is difficult to judge whether
this is reasonable. As with galaxies, in principle a wide range of DM
fractions would be possible, so the agreement with MOND is not a priori
expected. We conclude that this agreement could pose a fine-tuning problem
similar to that faced by $\Lambda$CDM in galaxies (Section 3). The problem is
less severe in galaxy groups due to the possibility of interlopers and the
smaller number of tracers in each group, among other difficulties.

##  7 Galaxy clusters

The largest gravitationally bound structures in the Universe are galaxy
clusters. Their importance to the missing gravity problem goes back almost a
century (Zwicky, 1933, 1937). The very large amounts of missing gravity
claimed in the earlier works were mostly due to an erroneous assumption that
most of the baryonic mass lies in the galaxies. Galaxy clusters actually
contain a very significant component of intracluster gas, but this only became
detectable much later with the advent of space-based X-ray telescopes
(Sarazin, 1986). However, this very substantial addition to the baryonic mass
budget of galaxy clusters is still not sufficient to resolve the missing
gravity problem on this scale (for a review, see Allen et al., 2011). In a
standard context, the missing gravity is assumed to come from CDM. As we
discuss below, this problem is reduced but not completely alleviated in MOND.
We will argue that this does not rule out the MOND scenario, but it does place
important constraints on how MOND can be made to work in a broader
cosmological context.

###  7.1 Internal dynamics

Galaxy clusters have internal accelerations that are typically of order
$a_{{}_{0}}$ or slightly above, so MOND can provide only a small amount of
extra gravity. It has long been known that this is insufficient to completely
solve the missing gravity problem in galaxy clusters (Sanders, 1999; Aguirre
et al., 2001; Sanders, 2003). It is sometimes argued that the MOND dynamical
mass is only twice the detectable $M_{s}$, making for a much smaller
discrepancy than with Newtonian gravity (e.g. Milgrom, 2015b). While this may
be correct, the discrepancy in MOND must still be understood. We argue it is
unlikely that the very substantial mismatch between the MOND dynamical mass
and the directly detected baryonic mass in galaxy clusters arises from some
form of baryonic matter that has thus far evaded detection (though see
Milgrom, 2015b).

An important clue is that if we assume the MOND dynamical mass is correct and
also that our inventory of baryons in galaxy clusters is nearly complete, then
the implied baryon fraction matches that required by standard fits to the CMB
power spectrum and the acoustic oscillations therein (Pointecouteau & Silk,
2005). Part of the cluster baryon budget is the intracluster light due to
stars outside galaxies (Gonzalez et al., 2007). Including this component,
Gonzalez et al. (2013) obtained baryon fractions consistent with the Planck
determination (Planck Collaboration VI, 2020). Other studies also assign
galaxy clusters a baryon fraction quite close to that implied by the CMB,
especially at the high-mass end (COSMOS Collaboration, 2009; McGaugh et al.,
2010). Since the accelerations there are generally large enough for MOND to
have only a small impact, it seems logical that the Universe as a whole
contains a dominant non-baryonic component.

Colliding galaxy clusters provide important clues regarding where the extra
gravity comes from, mainly because they allow for more detailed studies beyond
just the radial profiles of the gravitational field and the enclosed baryonic
mass. A crucial observation in this regard is the Bullet Cluster (discovered
by Tucker et al., 1995). Its interacting nature was noted by Tucker et al.
(1998), who identified the Bullet sub-component and thereby revealed the
ongoing interaction. This is also apparent in Chandra X-ray observations of a
clear bow shock (Markevitch et al., 2002). The importance of the Bullet
Cluster stems from the offset between its X-ray and weak lensing peaks, which
has been detected at high significance (Clowe et al., 2004, 2006).

The conventional interpretation is that the clusters consisted of hot gas
embedded in DM halos. When the clusters collided at high velocity, the gas
experienced substantial hydrodynamical drag, causing it to slow down. However,
the collisionless DM halos were not decelerated in this manner. As a result,
the DM halos reached a larger post-interaction separation than the gas,
explaining why the weak lensing peaks are more widely separated than the X-ray
peaks. This morphology can be reproduced in idealized non-cosmological
Newtonian simulations of the merger (Lage & Farrar, 2014).

It is of course not wise to place undue emphasis on any individual system,
especially when the Train Wreck cluster seems to have the opposite properties
to the Bullet Cluster (as briefly reviewed in Deshev et al., 2017). One reason
is that filaments near a system in projection can have non-trivial effects on
the weak lensing signature in MOND due to its non-linearity (Feix et al.,
2008). In any case, structures along the LOS can also create a weak lensing
peak, even if this may seem a lucky coincidence that is not the most favoured
interpretation. However, a chance alignment of this sort is often advocated
for the Train Wreck cluster because a central concentration of DM would not be
expected. In general, disagreement of a theory with observations can often be
attributed to errors with the observations or their interpretation, but this
logic is rarely applied when a popular theory agrees with observations. Just
as a claimed falsification of a theory might be a mistake, so too can a
seemingly correct prediction be a mistake caused by a coincidence between
incorrect observations and incorrect predictions. Nonetheless, it could be
argued that only one chance alignment is needed in the Train Wreck cluster but
two are needed in the Bullet Cluster to explain its two observed weak lensing
peaks. It is therefore probable that there is indeed an offset between its
X-ray gas and $\nabla\cdot\bm{g}$, the source of the gravitational field. Such
offsets have also been detected in 37 other clusters (Shan et al., 2010).

Figure 29: Weak lensing convergence map of the Bullet Cluster, shown using
contours in a MOND model with 2 eV/$c^{2}$ sterile neutrinos (solid black
lines) and as observed (dotted red lines). A four-centred potential is used,
with two centres corresponding to the gas (red stars labelled XR) and two
centres corresponding to the neutrinos (red stars labelled CM). The inset
shows the gas surface density in this model. Notice the good match to the weak
lensing data (Clowe et al., 2006). Reproduced from figure 1 of Angus et al.
(2007) using ApJ author republication rights.

Mainly due to this offset, the Bullet Cluster is quite difficult to reproduce
in MOND using only its directly detected baryonic matter (Knebe et al., 2009).
However, the Bullet Cluster can be reconciled with MOND using extra
collisionless matter in the form of neutrinos with a mass of 2 eV/$c^{2}$
(Angus et al., 2007). A large offset between the weak lensing and X-ray peaks
is clearly evident in their figure 1, which we reproduce as our Figure 29.
Angus et al. (2007) suggested that ordinary neutrinos could be sufficient, but
2 eV/$c^{2}$ is now known to be above the Katrin cosmology-independent upper
limit on the ordinary neutrino mass, which is only 0.8 eV/$c^{2}$ at 90%
confidence (Katrin Collaboration, 2019, 2022). Moreover, even neutrinos with a
mass of 7 eV/$c^{2}$ would be insufficient to explain the lensing caused by
some galaxy clusters in a MOND context (Natarajan & Zhao, 2008).

These difficulties led to the hypothesis of an undiscovered sterile neutrino
with a mass of 11 eV/$c^{2}$ (Angus, 2009). While significantly more massive
than the ordinary neutrinos, this is still so low that galaxy RC fits would
not be affected (Angus, 2010). This is not due to the high velocities of the
neutrinos $-$ the CMB temperature of $T=2.725$ K corresponds to an energy of
$kT=2.35\times 10^{-4}$ eV, implying that the Universe has expanded
$47000\times$ since any such neutrinos became non-relativistic. In the absence
of structure, their typical velocities today would therefore be about
$c/47000=6.4$ km/s. This is only $\approx 1\%$ of the local escape velocity
from the MW (Williams et al., 2017; Monari et al., 2018). Sterile neutrinos
with a mass of 11 eV/$c^{2}$ are nonetheless unable to form substantial
galactic halos around MW-like discs due to the Tremaine-Gunn limit (Tremaine &
Gunn, 1979). The important constraint is the phase space density of the
neutrinos, which are subject to the Pauli Exclusion Principle. As a result,
halos of gravitationally bound sterile neutrinos are actually expected around
galaxies, but they would be dynamically irrelevant even if all available phase
space states were filled. Galaxy clusters are not only larger, they also have
a deeper potential and thus a higher velocity dispersion, increasing the
available phase space volume. Therefore, it is possible for sterile neutrinos
to be dynamically relevant in galaxy clusters but not in galaxies.

Assuming MOND is the correct theory, we can infer the mass distribution and
subtract the baryonic component, thereby identifying where the sterile
neutrinos must be hiding. This was attempted by Angus et al. (2010), who then
addressed whether the inferred density distribution of the sterile neutrinos
is consistent with the Tremaine-Gunn limit. Strikingly, the neutrinos just
reach this limit at the centre of the cluster in all 30 cases considered, at
least if plausible assumptions are made regarding the $M_{\star}/L$ of the
brightest cluster galaxy. This suggests that the sterile neutrinos accreted
onto a halo until they reached the Tremaine-Gunn limit, preventing further
accretion (at least in the central regions). While some other interpretation
could be advanced, we argue that it would be a rather unlikely coincidence if
the observations are caused by CDM particles with a vastly higher Tremaine-
Gunn limit to their mass density.

The sterile neutrino solution to the MOND cluster problem is also suggested by
the MOND galaxy cluster analysis of Ettori et al. (2019). Though the overall
discrepancy (at large radii) is smaller than in previous studies, their figure
6 shows that the extra dynamical mass required in MOND is in a centrally
concentrated component, as might be expected for a dense sterile neutrino core
at the Tremaine-Gunn limit. Moreover, if sterile neutrinos with a mass of 11
eV/$c^{2}$ were thermalized in the early universe, they would have the same
relic abundance as the CDM in the $\Lambda$CDM paradigm, thus explaining the
acoustic oscillations in the CMB (Section 9.2.1; see also Angus, 2009;
Haslbauer et al., 2020).

There are also some hints for sterile neutrinos in terrestrial experiments
(for a review, see Giunti & Lasserre, 2019). While the latest results
disfavour the sterile neutrino interpretation of the anomalies in short-
baseline experiments (MicroBooNE collaboration, 2022), this is far from ruled
out at present (Argüelles et al., 2022a). One reason for the apparent
inconsistencies between different terrestrial experiments could be the
approximation of the neutrinos as a plane wave (Argüelles et al., 2022b). A
recent paper considering a wide range of experiments found $>5\sigma$ evidence
for sterile neutrinos (Berryman et al., 2022). This is largely due to neutrino
flux measurements at two distances from the radioactive chromium-51 source in
the Baksan Experiment on Sterile Transitions (BEST; Barinov & Gorbunov, 2022;
Barinov et al., 2022a, b). These experiments favour sterile neutrino masses
significantly above 1 eV/$c^{2}$, a region of parameter space which is often
not analysed due to the difficulties this would pose for $\Lambda$CDM.
However, it is important to conduct particle physics experiments without being
biased by prior beliefs regarding the law of gravity at low accelerations.

Substructures identified in galaxy cluster lenses could pose a tight
constraint on sterile neutrinos due to the Tremaine-Gunn limit. In this
regard, the straight arc around Abell 2390 was found to be consistent with 11
eV/$c^{2}$ sterile neutrinos in a MOND context (Feix et al., 2010). Other
cluster-scale lenses could pose further constraints. In principle, the sterile
neutrinos could have a higher mass, but their mass should be $\la 300$
eV/$c^{2}$ to avoid disrupting MOND fits to the internal dynamics of galaxies
(Angus, 2010).

If MOND is correct, some explanation must be advanced for why galaxies obey it
but galaxy clusters seemingly do not. We have argued above that this is
because galaxy clusters have much deeper potentials and larger sizes, both of
which increase the phase space volume available to sterile neutrinos. As a
result, sterile neutrinos at the Tremaine-Gunn limit would significantly
affect galaxy clusters but not galaxies.

The deeper potentials in galaxy cluster environments may lead to a new regime
of gravity that is not apparent in galaxies. This has given rise to the
hypothesis known as extended MOND (EMOND; Zhao & Famaey, 2012), which promotes
$a_{{}_{0}}$ to a dynamical variable that takes on higher values in regions
with a deeper potential. This idea was explored further by Hodson & Zhao
(2017a). While initial results were somewhat promising, we prefer the sterile
neutrino interpretation because it would also explain the high third peak in
the CMB power spectrum (Section 9.2.1).

We should of course bear in mind the possibility that some modification to
gravity can account for these observations using only the directly observed
$M_{s}$. For example, MOG can account for the weak lensing properties of the
Bullet Cluster using only its observed mass (Brownstein & Moffat, 2007). While
MOG has now been excluded at high significance by many interlocking lines of
evidence (Section 3.6), it remains possible that some other modification to GR
could satisfy the relevant constraints (e.g. Acedo, 2017). This is suggested
by the fact that the total matter density in $\Lambda$CDM exceeds the baryonic
density by a factor very close to $2\mathrm{\pi}$ rather than by many orders
of magnitude, as might be expected if the mean baryon and CDM densities were
set by very different physics (Milgrom, 2020).

The cluster-scale difficulties with modified gravity approaches can be
interpreted as a sign that we should return to the $\Lambda$CDM picture, which
can statistically account for the CMB anisotropies (Planck Collaboration VI,
2020) and seems able to account for the internal dynamics of galaxy clusters
thanks to the presence of CDM in large amounts. The problem then becomes how
to explain the successes of MOND (e.g. Famaey & McGaugh, 2012; Lelli et al.,
2017b; Shelest & Lelli, 2020) and the failures of $\Lambda$CDM (e.g. Kroupa,
2012, 2015; Peebles, 2020). A hybrid approach (elaborated in Section 9.2)
seems to be the best way to explain available observations across all scales.

###  7.2 Probing structure formation

The distribution of galaxy clusters is sensitive to the growth of structure
from the small observed density perturbations in the CMB. The observed low-
redshift cluster mass function seems broadly consistent with $\Lambda$CDM over
the mass range $3\times 10^{12}M_{\odot}-5\times 10^{14}M_{\odot}$ (Böhringer
et al., 2017), though this is subject to the normalization of the relation
between cluster mass and X-ray luminosity (Böhringer et al., 2014; Schrabback
et al., 2021a). The general expectation in MOND is that there should be
deviations from $\Lambda$CDM in terms of galaxy clusters forming with too high
a mass too early in cosmic history (Sanders, 2001). This is also apparent from
$N$-body simulations (Angus et al., 2013; Katz et al., 2013) using a MOND
cosmology that we discuss later (Section 9.2.2).

The Bullet Cluster was thought to provide just such a challenge. Initial
modelling attempts implied a very high collision velocity (Milosavljević et
al., 2007; Springel & Farrar, 2007; Mastropietro & Burkert, 2008). This leads
to significant tension with $\Lambda$CDM expectations (Lee & Komatsu, 2010;
Thompson & Nagamine, 2012). It was later shown that the morphology can be
reproduced in idealized non-cosmological $\Lambda$CDM simulations of the
merger with a lower infall velocity (Lage & Farrar, 2014). The parameters of
their simulation are not thought to be seriously problematic for $\Lambda$CDM
(Watson et al., 2014; Thompson et al., 2015). Indeed, it has been estimated
that 0.1 analogues are expected over the full sky out to the Bullet Cluster
redshift of 0.30 (Kraljic & Sarkar, 2015).

A much more problematic galaxy cluster collision is El Gordo (Menanteau et
al., 2010, 2012). This has a higher redshift, mass (Jee et al., 2014), and
ratio of infall velocity to escape velocity (Molnar & Broadhurst, 2015; Ng et
al., 2015; Zhang et al., 2015). Collisions between individually rare massive
galaxy clusters should be quite unlikely as it is necessary to form two
massive clusters within striking distance. By considering the properties of
galaxy cluster pairs that have turned around from the cosmic expansion in the
Hubble-volume Jubilee cosmological $\Lambda$CDM simulation (Watson et al.,
2013), it was shown that El Gordo rules out $\Lambda$CDM at $6.16\sigma$
significance (Asencio, Banik & Kroupa, 2021). The falsification is illustrated
in their figure 7, which we reproduce as our Figure 30. This shows the mass-
redshift distribution along our past lightcone of galaxy cluster pairs that
have turned around from the cosmic expansion with similar dimensionless
parameters to El Gordo. The solid red contour passes through the parameters
for the pre-merger configuration of El Gordo, shown with a red cross.

The pre-merger configuration was obtained from non-cosmological hydrodynamical
simulations of the merger (Zhang et al., 2015), with their best-fitting Model
B used in the main analysis. Those authors were aware that their preferred
pre-merger parameters would cause tension with $\Lambda$CDM, motivating them
to consider a wide range of model parameters and check which ones might be
compatible with observations of El Gordo. They ran 123 models altogether, with
their preferred parameters in agreement with previous studies. But some other
models also gave plausible fits to the observations at some epoch if viewed
from some direction. Consequently, it is important to look beyond their best-
fitting model and consider a wider range of possible pre-merger configurations
that evolve into something resembling the observed morphology of El Gordo. To
this end, the parameters of the Zhang et al. (2015) Model A were also
considered by Asencio et al. (2021), with the lower mass making the pre-merger
configuration more likely in Jubilee. However, the tension in this case is
still $5.14\sigma$, so the initial configuration of the model could not
plausibly arise in $\Lambda$CDM cosmology. If the pre-merger parameters are
pushed into the range where consistency is gained, then it becomes impossible
to explain the observed morphology of El Gordo. For instance, using too low a
mass would cause disagreement with the observed X-ray flux, while too low a
velocity would be inconsistent with the observed two-tailed morphology. It
appears impossible for a Newtonian hydrodynamical simulation to reproduce the
observed morphology of El Gordo using a pre-merger configuration that is
plausible in $\Lambda$CDM.

Figure 30: The mass-redshift distribution in $\Lambda$CDM of El Gordo-like
galaxy cluster collisions in a 755 deg2 region of sky along our past
lightcone. The probability density is shown using the indicated colour scale
in the parameter space of the cosmic scale factor $a$ when the collision
occurred and the total virial mass $M$ of the clusters, expressed here as
$\widetilde{M}\equiv\log_{10}\left(M/M_{\odot}\right)$. The red cross shows
the El Gordo parameters, with the error bar showing a 20% uncertainty on its
$M$. The total probability outside each contour is used to find its level of
significance, which is indicated on the figure. The contour passing through
the El Gordo parameters is shown with a solid red line representing the locus
of $6.16\sigma$ outliers to $\Lambda$CDM expectations, the probability beyond
which is $7.51\times 10^{-10}$. Reproduced from figure 7 of Asencio et al.
(2021).

The high El Gordo mass was recently confirmed in a detailed $\Lambda$CDM weak
lensing analysis that gives it a total mass of $2.13^{+0.25}_{-0.23}\times
10^{15}M_{\odot}$ (Kim et al., 2021b). While slightly lower than the
$3.2\times 10^{15}M_{\odot}$ assumed in Model B of Zhang et al. (2015), it is
slightly higher than the $1.95\times 10^{15}M_{\odot}$ assumed in their Model
A, which as mentioned above has initial conditions incompatible with
$\Lambda$CDM cosmology. It is therefore very surprising that Kim et al.
(2021b) claim their weak lensing measurements reconcile El Gordo with
$\Lambda$CDM expectations. A closer examination of their work reveals that
this is certainly not the case because the slightly lower weak lensing mass
only slightly reduces the tension for $\Lambda$CDM. Instead, the main reason
for this claim is the use of a very low pre-merger infall velocity, which is
not directly measured in a weak lensing analysis. The authors attempted to fit
some aspects of the El Gordo morphology using a simplified analysis that does
not include dynamical friction, which would be quite considerable for galaxy
clusters as massive as the El Gordo progenitors. This might explain why Kim et
al. (2021b) were unable to reproduce the observed locations of the radio
relics. It is possible that including dynamical friction would lead to a
better match, but then the preferred model parameters might differ greatly
from what they obtained. A detailed parameter study of El Gordo using
hydrodynamical simulations was indeed conducted previously by Zhang et al.
(2015). The only problem identified by Kim et al. (2021b) with this study was
that the collision velocities were too high for $\Lambda$CDM, which might
naturally arise in the MOND scenario due to an isolated point mass having a
logarithmic potential in the deep-MOND limit. A Newtonian simulation of the
interaction should still be quite accurate because the accelerations are above
$a_{{}_{0}}$ during the interaction phase and because galaxy clusters are
expected to contain a dominant dark matter component even in MOND (Section
9.2). Zhang et al. (2015) were aware that their model parameters would be in
tension with $\Lambda$CDM expectations, so they already tried lowering the
mass and infall velocity of the El Gordo progenitors. This is possible to only
a limited extent while remaining consistent with the observed morphology of El
Gordo. Therefore, we reject the claim of Kim et al. (2021b) that El Gordo is
compatible with $\Lambda$CDM $-$ there are currently no Newtonian
hydrodynamical simulations of it that match its observed properties and
morphology with initial conditions that might plausibly arise in the
$\Lambda$CDM framework. Rather, the great many hydrodynamical simulations that
have been done so far indicate that the observed properties of El Gordo can
arise only for a pre-merger configuration which is not feasible in
$\Lambda$CDM. Some of these simulations use a slightly lower El Gordo mass as
favoured by Kim et al. (2021b). Moreover, those authors seem to favour the
returning scenario for El Gordo whereby the clusters are observed when
approaching each other for a second time. This would increase the elapsed time
since the clusters first came into contact with each other, so the pre-merger
configuration would need to exist at an even earlier stage in cosmic history,
worsening the problem for $\Lambda$CDM (Figure 30).

In addition to El Gordo, there are several other galaxy clusters which might
be problematic for $\Lambda$CDM (Asencio et al., 2021). For instance, scaling
the Bullet Cluster results of Kraljic & Sarkar (2015) to the sky area of its
discovery survey reveals that it presents $2.78\sigma$ tension (Asencio et
al., 2021). Combined with El Gordo, the tension rises to $6.43\sigma$ (see
their section 3.4). The tension is so serious that it cannot be resolved by
the discovery of no further problematic clusters for $\Lambda$CDM in the rest
of the sky.

An important aspect of Asencio et al. (2021) was the use of previously
conducted $N$-body cosmological MOND simulations (Katz et al., 2013) to
estimate the occurrence rate of El Gordo analogues in MOND. Scaling their
results to the effective volume of the survey in which El Gordo was
discovered, we expect $\approx 1.2$ El Gordo analogues in MOND cosmology
(Asencio et al., 2021). This is quite consistent with observations (see their
section 4.3). However, the expected number of El Gordo analogues in
$\Lambda$CDM is only $7.51\times 10^{-10}$ in the survey region (see their
section 3.3).

This demonstrates that the successes of MOND and the failures of $\Lambda$CDM
extend well beyond the low-redshift Universe and the galaxy RC data which
originally motivated the MOND hypothesis (Milgrom, 1983a). El Gordo is on a
large enough scale that it should be well resolved in cosmological
simulations. Moreover, the large-scale shock features would be little affected
by galaxy-scale baryonic feedback processes. Given the many detailed
hydrodynamical simulations of El Gordo done in the past (Donnert, 2014; Molnar
& Broadhurst, 2015; Ng et al., 2015; Zhang et al., 2015), it seems extremely
unlikely that a Newtonian simulation of it will be found to reproduce its
detailed morphological properties while also using a pre-merger configuration
consistent with $\Lambda$CDM cosmology. The very large size of the Hubble-
volume Jubilee cosmological simulation and the large number of cluster pairs
used in the analysis of Asencio et al. (2021) argues against uncertainties on
the $\Lambda$CDM side. Indeed, figure 8 of Böhringer et al. (2017) shows that
very few galaxy clusters with the El Gordo mass are expected at its redshift
of 0.87 (Menanteau et al., 2012). Since El Gordo is actually two interacting
clusters, this would be even less likely, and is perhaps the main reason why
the Bullet Cluster and especially El Gordo are so problematic for
$\Lambda$CDM.

##  8 Large-scale structure

Since disc galaxy outskirts are already in the low-acceleration regime, tests
of MOND should be possible on scales larger than individual galaxies. In
addition to virialized systems, the large-scale structure of the Universe
could provide important constraints. We defer a discussion of the CMB to
Section 9 because this does not presently help to distinguish between the
$\Lambda$CDM and MOND models. For the moment, we focus on later epochs.

It was argued that MOND has difficulty reproducing the internal properties of
spectroscopically detected Lyman-$\alpha$ (Ly-$\alpha$) absorbers (the
Ly-$\alpha$ forest; Aguirre et al., 2001). However, those authors noted that
consistency could be gained if the EFE is considered. In general, the EFE
should have been stronger at higher redshift because the less pronounced
cosmic structures at earlier times would be more than compensated by a smaller
Universe (e.g. Asencio et al., 2021). Those authors estimated that
$g_{{}_{e}}\mathrel{\vbox{ \offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern