across almost six orders of magnitude in $M_{s}$ (Wu & Kroupa, 2015).

In a $\Lambda$CDM context, the tightness of the RAR is in $3.6\sigma$ tension
with expectations (Desmond, 2017b) based on a previously developed statistical
framework (Desmond, 2017a). A puzzling aspect is the diversity of RC shapes at
fixed $v_{{}_{f}}$, which essentially fixes halo properties like the virial
mass (Oman et al., 2015). In dwarf galaxies with a large proportion of missing
gravity, the baryons are sub-dominant tracers of the potential. Consequently,
they should reveal the same RC in galaxies with the same $v_{{}_{f}}$, even if
the baryonic distribution differs. Since this expectation is not correct, one
can postulate that baryonic feedback processes drove large fluctuations in the
potential depth and that this led to redistribution of the dominant DM
component (Pontzen & Governato, 2012). However, the process would be
stochastic, especially in dwarf galaxies where e.g. only a few extra
supernovae could make a large difference. As a result, it is very difficult to
obtain a tight RAR.

In this regard, it is surprising that the RAR can apparently be reproduced in
$\Lambda$CDM (Dutton et al., 2019). However, those authors did not consider
the above-mentioned diversity problem. Section 4.2.1 of Ghari et al. (2019b)
demonstrates that it is extremely difficult to get a tight RAR simultaneously
with a diversity of dwarf galaxy RCs at fixed $v_{{}_{f}}$. This diversity can
readily be understood in MOND as a consequence of differences in the baryonic
surface density relative to the MOND threshold (Equation 52), as discussed
further in McGaugh (2020). Essentially, the MOND RC of a deep-MOND galaxy has
the same characteristic length scale as the baryon distribution, so
differences in this at fixed $M_{s}$ (and thus $v_{{}_{f}}$) nicely account
for differences in the observed RC. For theories where a tight RAR is not
obvious, claims to reproduce the RAR need to consider a wide range in surface
brightness at fixed $M_{s}$.

Dwarf galaxies can be challenging to resolve in a large cosmological
simulation. Results should be more reliable for larger galaxies, especially if
the galaxies are resimulated with better resolution. One such attempt is the
zoom-in $\Lambda$CDM resimulation project known as AURIGA, which considers 30
MW-like galaxies (Grand et al., 2017). This has difficulty reproducing the
BTFR, as is evident in their figure 11. In addition to only covering the range
$v_{{}_{f}}\ga 160$ km/s and thus not addressing lower mass galaxies, the
figure reveals that $v_{{}_{f}}$ rises more steeply with $M_{s}$ than the
observed slope, which is very close to the MOND expectation of $1/4$ (McGaugh
& Schombert, 2015). The sample size is perhaps too small to reliably address
the question of intrinsic scatter, though that should be checked in subsequent
work.

Despite a vast amount of effort by the $\Lambda$CDM community, it remains
extremely challenging to match the small intrinsic scatter of the RAR over a
very wide range of galaxy properties including mass, size, surface density,
and gas fraction (figure 4 of Lelli et al., 2017b). The latter is certainly a
concern if the RAR is explained by loss of baryons from galaxies due to
stellar feedback, because then one might expect galaxies that formed more
stars to have lost a larger fraction of their baryons, thus deviating from the
BTFR defined by gas-rich galaxies. However, gas-rich and gas-poor galaxies lie
at the same location on a BTFR diagram (see e.g. figure 3 of Famaey & McGaugh,
2012). More generally, a tight RAR implies that the fraction of the total
gravity that is missing in the outer region of a high-mass galaxy is the same
as in the inner region of a low-mass galaxy provided the two locations have
the same $g_{r}$. It remains very difficult to understand this if the
proportion of missing gravity depends strongly on how many baryons were
ejected from the galaxy by stochastic feedback processes like individual
supernovae and major mergers. Other workers like Salucci (2019) have also
pointed out the tight relation between the luminous and dark components in
galaxies if their dynamical discrepancies are attributed to CDM particles. The
RAR is much easier to understand as a fundamental relation not contingent on
the formation history of a galaxy, much like Newtonian gravity accounts for
the motions of planets and moons in the Solar System across a vast range of
mass and size without knowing how these objects formed. It would be unusual if
the generally favoured interpretation of Solar System motions was that they
were governed by an inverse cube gravity law supplemented by a distribution of
invisible mass, the end result being identical to the action of an inverse
square gravity law sourced only by the known masses.

####  3.1.2 Vertical dynamics

In principle, the $\Lambda$CDM paradigm can match the RC of a galaxy
satisfying the RAR if its DM halo is tuned to have a particular density
profile. The RAR is also naturally explained in MOND. We can break this
degeneracy by considering the vertical gravity $g_{z}$ felt by a star outside
the disc mid-plane in the direction parallel to the disc spin axis. The basic
principle is that very little of the mass in any CDM halo would lie close to
the disc mid-plane because the DM is expected to be collisionless. Given the
significant amount of DM required around the MW, postulating it to be
collisional would create a DM disc, causing tension with observations (Buch,
Leung & Fan, 2019). Therefore, measurements of $g_{z}$ close to the disc mid-
plane should reveal only a small amount of missing gravity commensurate with
the amount of DM located even closer to the mid-plane. Since the scale height
of the disc should be far smaller than the extent of the DM halo, it is
possible for a star to lie outside the disc and yet have its $g_{z}$ little
affected by the DM halo, even when this dominates $g_{r}$.

In MOND, $g_{z}$ would be enhanced over Newtonian expectations by the local
value of $\nu$, which can greatly exceed 1 if the surface density is
sufficiently small. Consequently, MOND typically predicts higher $g_{z}$ just
outside the baryonic disc, as calculated for the MW by Bienaymé et al. (2009).
The higher predicted $g_{z}$ in MOND implies a higher divergence to the
gravitational field, which can be attributed to a disc of PDM. The PDM
distribution of thin exponential discs was visualized in e.g. Lüghausen et al.
(2013); Lüghausen et al. (2015).

Observationally, $g_{z}$ can be inferred from the vertical profile of the
vertical velocity dispersion $\sigma_{z}$. Combined with measurements of the
tracer density $\rho_{t}$, the vertical momentum flux through a surface of
constant height is $\rho_{t}{\sigma_{z}}^{2}$ per unit area. Considering a
slab of finite thickness, the momentum flux though the top and bottom surfaces
will differ slightly. If the galaxy is in equilibrium, this difference should
be balanced by the slab’s thickness multiplied by $\rho_{t}g_{z}$, which can
therefore be determined from the observed kinematics. If we are considering
the Solar neighbourhood of the MW, a simplification is that
$\sigma_{z}\approx\sigma_{{}_{\text{LOS}}}$ because observations towards the
Galactic poles will have the LOS aligned similarly to the desired component of
the velocity dispersion. Even if this is not exactly correct, the corrections
will be small provided the LOS is nearly along or opposite the disc spin axis,
reducing sensitivity to proper motion uncertainties.

Inferring $g_{z}$ in this manner fundamentally relies on measuring gradients
in $\rho_{t}$ and $\sigma_{z}$. Differentiating the data necessarily increases
the uncertainties. Moreover, the assumption of dynamical equilibrium may not
be accurate. Nonetheless, there have been several attempts to test the MOND-
predicted enhancement to $\sigma_{z}$. There is currently no evidence for this
enhancement, with the data actually supporting the $\Lambda$CDM model (Lisanti
et al., 2019). Their figure 1 shows that MOND overpredicts $g_{z}$ by
$2\sigma$, so this scenario cannot be ruled out either due to the still
significant uncertainty. Other works also reported difficulties when
attempting to measure $g_{z}$ (Hagen & Helmi, 2018), partly because of an
asymmetry between determinations towards the north and south (Guo et al.,
2020) that indicates non-equilibrium processes like the passage of Sagittarius
through the Galactic disc (see their section 5.3.3).

Another consequence of the enhanced $g_{z}$ in MOND is that the velocity
dispersion tensor would appear to tilt to a different extent as one moves away
from the disc mid-plane. In general, the stronger the disc self-gravity, the
closer the potential to cylindrical symmetry. If the disc self-gravity is
weaker, then the geometry is closer to spherical symmetry. This has been
proposed as another related test of MOND (Bienaymé et al., 2009).
Observational measurements of the velocity ellipsoid’s orientation remain
uncertain (Everall et al., 2019), partly because of sensitivity to the Gaia
parallax zero-point offset (Hagen et al., 2019). This should be clarified with
future Gaia data releases, though departure from equilibrium would remain a
concern. It may be necessary to address this using simulations of Sagittarius
crossing the disc mid-plane in $\Lambda$CDM and MOND, with the disc response
perhaps being a better discriminant (Laporte et al., 2018, 2019; Bennett et
al., 2022).

Looking beyond the MW, the DiskMass survey (Bershady et al., 2010) attempted
to measure $g_{z}$ statistically by combining $\sigma_{{}_{\text{LOS}}}$
measurements from face-on galaxies with scale height measurements from edge-on
galaxies. The idea was that the scale heights of the face-on galaxies would be
known statistically based on the relation between scale length and scale
height in edge-on galaxies. The analysis gave very low values for $g_{z}$
(Angus et al., 2015), though this could be due to underestimation of
$\sigma_{{}_{\text{LOS}}}$ (Angus, 2016). It was later shown that this is
entirely possible because luminosity-weighted $\sigma_{{}_{\text{LOS}}}$
measurements give a higher weighting to more massive stars, but the mass-
weighted $\sigma_{{}_{\text{LOS}}}$ entering a dynamical analysis is more
sensitive to less massive stars (Aniyan et al., 2016; Aniyan et al., 2018). In
a galaxy like NGC 6946 that is still forming stars, this would cause a
mismatch between the typical ages of stars used to determine
$\sigma_{{}_{\text{LOS}}}$ and the generally older stars that comprise the
bulk of the mass (Aniyan et al., 2021). As the stars in disc galaxies should
be getting dynamically hotter with time due to interactions with molecular
clouds (Aumer et al., 2016) and other processes, the mass-weighted
$\sigma_{{}_{\text{LOS}}}$ could indeed be higher than the luminosity-weighted
$\sigma_{{}_{\text{LOS}}}$ if appropriate precautions are not taken. As a
result, it has been argued that the DiskMass results are quite consistent with
MOND (Milgrom, 2018b).

It is also possible to measure $\sigma_{z}$ from the gas. This would be
simpler in some ways as the velocity dispersion tensor should be isotropic, so
the HI velocity dispersion $\sigma_{{}_{\text{HI}}}\approx\sigma_{z}$. Using
this approach, Sánchez-Salcedo et al. (2008) found plausible agreement between
MOND and the observed flaring of the MW disc over Galactocentric distances of
$R=16-40$ kpc, though the observed thickness exceeds the predicted value by
$\approx 70\%$ at $R=10-16$ kpc. Those authors attributed the mismatch to non-
thermal sources of pressure support, which they argued would be quite
plausible.

Gas kinematic and scale height measurements are also possible in external
galaxies viewed close to edge-on. This technique was attempted by Das et al.
(2020) in four large galaxies and three dwarfs, with the latter probing the
low-acceleration regime. The velocity dispersions in these cases are higher
than expected in $\Lambda$CDM, leading the authors to conclude that a dark
disc is present. This is qualitatively consistent with the phantom disc
predicted by MOND, but more detailed analysis is required to check if MOND can
explain the vertical dynamics of these galaxies. Even if $\Lambda$CDM
underpredicts their $g_{z}$, it is possible that MOND will overpredict their
$g_{z}$.

We conclude that the vertical dynamics of disc galaxies do not currently break
the CDM-MOND degeneracy, but this remains a promising avenue for further
investigation.

###  3.2 Elliptical galaxies and dwarf spheroidals

Disc galaxies have historically been the focus of attempts to test MOND due to
their relative simplicity from a theoretical perspective and the relative ease
of interpreting the observations (Section 3.1). In principle, the
gravitational fields of elliptical galaxies can be calculated in a similar
manner to thin discs (Equation 21 or 24). The task is actually simpler for
spherical galaxies. Considering them makes it possible to study the missing
gravity problem in dwarf galaxies, which tend to be pressure-supported and
often have quite low internal accelerations. In this section, we consider what
can be learned from the internal dynamics of elliptical and dwarf spheroidal
galaxies using a variety of tracers to probe their potential.

####  3.2.1 Velocity dispersion

A test particle on a near-circular orbit has a well-known relation between
rotational velocity $v_{c}$ and gravity $g$, namely $g={v_{c}}^{2}/r$. The
relation is more complicated for a pressure-supported system, where the
typically eccentric orbit of each particle means that its instantaneous
velocity is not a meaningful constraint on the potential. Instead, one must
consider the problem statistically. This can be done using the MOND
generalization of the virial theorem (Milgrom, 1994b). Its equation 14 states
that for an isolated system in the deep-MOND limit with an isotropic velocity
dispersion tensor, the mass-weighted LOS velocity dispersion

$\displaystyle\sigma_{{}_{\text{LOS}}}~{}=~{}\sqrt[{}^{4}]{\frac{4GMa_{{}_{0}}}{81}}\,.$
(56)

This applies to any modified gravity theory of MOND (Milgrom, 2014b), thereby
providing a good starting point for estimating what MOND predicts about an
isolated pressure-supported low-acceleration system. The 3D velocity
dispersion is $\sigma_{{}_{\text{LOS}}}\sqrt{3}$, so formulae in the
literature may differ depending on which measure of velocity dispersion is
intended. Observational studies generally focus on $\sigma_{{}_{\text{LOS}}}$.
As with Newtonian gravity, the application of the virial theorem is subject to
the usual caveats, for instance that the system should be in equilibrium.
Moreover, the velocity dispersion tensor could be anisotropic, which would
require a more complicated treatment such as solving the Jeans equations
(Equation 60 in spherical symmetry). In this case, the 3D velocity dispersion
would still be $\sqrt[{}^{4}]{4GMa_{{}_{0}}/9}$.

An important aspect of Equation 56 is that the MOND dynamical mass
$\propto{\sigma_{{}_{\text{LOS}}}}^{4}$, which is much steeper than the
Newtonian scaling of ${\sigma_{{}_{\text{LOS}}}}^{2}$. As a result, the early
claim of a severe discrepancy between MOND and the observed
$\sigma_{{}_{\text{LOS}}}$ for some dwarfs (Gerhard & Spergel, 1992) was later
resolved through better measurements and by properly taking into account their
uncertainties (Milgrom, 1995). This is just one of many examples where claims
to have falsified MOND were later shown to be premature, with better data
actually agreeing quite well with MOND.

The premature claims sometimes rely on an incomplete treatment of MOND, which
reduces to Equation 56 only for isolated deep-MOND systems in virial
equilibrium with an isotropic velocity dispersion tensor. Haghi et al. (2019a)
used $N$-body simulations conducted by Haghi et al. (2009) to overcome two
important limitations of Equation 56, namely the assumption of having an
isolated system in the deep-MOND limit. Haghi et al. (2019a) analytically fit
the earlier $N$-body results in their section 2, coming up with a very useful
series of formulae to predict $\sigma_{{}_{\text{LOS}}}$ for a system where
the external field is non-negligible and where the gravity is a possibly
significant fraction of $a_{{}_{0}}$. Their formulae yield the correct
asymptotic limits. Though the underlying $N$-body simulations are based on an
interpolating function with a sharper transition between the Newtonian and
Milgromian regimes than Equation 19, this should have only a modest effect on
the results, and even then only if the typical gravity is close to
$a_{{}_{0}}$ (Haghi et al., 2009). In this case, we expect
$\sigma_{{}_{\text{LOS}}}$ to be slightly higher with a more realistic
interpolating function.

It is also possible to estimate how the $\sigma_{{}_{\text{LOS}}}$ implied by
Equation 56 should be altered to include departures from isolation and the
deep-MOND limit, but without performing $N$-body simulations (Famaey et al.,
2018). While their approach is less rigorous, it perhaps gives a more
intuitive understanding of the corrections, which in their work are based on
equation 59 of Famaey & McGaugh (2012). This uses the AQUAL formulation
(Equation 21) with the simple interpolating function (Equation 18). Their
prescription is to solve for the internal gravity $g_{{}_{i}}$ implicitly
using

$\displaystyle\left(g_{{}_{i}}+g_{{}_{e}}\right)\mu\left(\frac{g_{{}_{i}}+g_{{}_{e}}}{a_{{}_{0}}}\right)~{}=~{}g_{{}_{N,i}}+g_{{}_{e}}\mu_{{}_{e}}\,,$
(57)

with $i$ and $e$ subscripts denoting internal and external contributions to
the gravity as before.

Figure 5: The $\sigma_{{}_{\text{LOS}}}$ of a stellar Plummer sphere (Plummer,
1911) with mass $M=2\times 10^{8}M_{\odot}$ for a range of internal Newtonian
gravity $a_{{}_{N}}\equiv GM/\left(2{r_{h}}^{2}\right)$, which is varied by
altering the half-mass radius $r_{h}$. Results for different external field
strengths are shown with different coloured lines, as indicated in the legend.
The solid lines show a fit to numerical AQUAL results (Haghi et al., 2009),
while the dotted lines are based on the ansatz proposed by Famaey & McGaugh
(2012) in their equation 59. The dashed grey horizontal line shows the
isolated deep-MOND result (Equation 56). For a different mass galaxy, the
values shown here should be scaled $\propto\sqrt[{}^{4}]{M}$. Notice how the
ansatz of Famaey & McGaugh (2012) systematically underpredicts the numerical
results when the internal and external gravity are both weak. The reason for
this deficiency is explained analytically in the text. Reproduced from figure
3 of Haghi et al. (2019a).

The numerical values given by this approach are very similar to those obtained
by Haghi et al. (2019a) based on fitting $N$-body results, as demonstrated in
their figure 3 (reproduced here as our Figure 5). This shows
$\sigma_{{}_{\text{LOS}}}$ for a spherical Plummer distribution of stars with
$M=2\times 10^{8}\,M_{\odot}$ for different internal and external
gravitational fields, considering both Equation 57 and the fit of Haghi et al.
(2019a) to the numerical AQUAL simulations of Haghi et al. (2009), which used
the standard interpolating function $\mu=x/\sqrt{1+x^{2}}$. The approaches
broadly agree, but a systematic offset is apparent when the internal and
external gravity are both very weak. The reason is that the angle-averaged
radial gravity should be boosted by a factor of
$\mathrm{\pi}/\left(4\mu_{{}_{e}}\right)$ if the EFE dominates, a consequence
of the AQUAL EFE-dominated potential being (Milgrom, 1986a):

$\displaystyle\Phi~{}=~{}-\frac{GM}{\mu_{{}_{e}}r\sqrt{1+L_{{}_{e}}\sin^{2}\theta}}\,.$
(58)

This is analogous to the QUMOND result in Equation 41, with $L_{{}_{e}}$ being
the logarithmic derivative of $\mu$ with respect to its argument. Since
$\mu_{{}_{e}}\nu_{{}_{e}}\equiv 1$, there are only slight differences between
the two approaches, e.g. the factor of $\mathrm{\pi}/{4}$ in AQUAL becomes
$5/6$ in QUMOND. Another useful relation is
$\left(1+L_{{}_{e}}\right)\left(1+K_{{}_{e}}\right)=1$ (equation 38 of Banik &
Zhao, 2018d). The EFE-dominated solutions in QUMOND and AQUAL were also
compared in Milgrom (2010). However, the factor of $\mathrm{\pi}/4$ is not
correctly obtained from Equation 57 in the appropriate asymptotic limit
$g_{{}_{i}}\ll g_{{}_{e}}$. By linearizing the $\mu$ function, it is
straightforward to show that Equation 57 reduces to

$\displaystyle\mu_{{}_{e}}\left(1+L_{{}_{e}}\right)g_{{}_{i}}~{}=~{}g_{{}_{N,i}}\,.$
(59)

Since $L_{{}_{e}}\to 1$ in the deep-MOND limit, the approach of Famaey &
McGaugh (2012) and Famaey et al. (2018) implies
$g_{{}_{i}}=g_{{}_{N,i}}/\left(2\mu_{{}_{e}}\right)$, yielding a factor of 0.5
where there should be $\approx 0.8$ in AQUAL or QUMOND. An extra factor of
$0.8/0.5=1.6$ in the gravity should translate to a factor of $\sqrt{1.6}$ in
the velocity dispersion, so we expect $\sigma_{{}_{\text{LOS}}}$ to be
$\approx 30\%$ higher in numerical simulations. This is indeed roughly the
case in Figure 5. We therefore recommend using the fit of Haghi et al. (2019a)
to numerical determinations of $\sigma_{{}_{\text{LOS}}}$ as shown in this
figure. Its results can be used directly if a curve is available for the
appropriate $g_{{}_{e}}$, bearing in mind that the plotted
$\sigma_{{}_{\text{LOS}}}$ should be scaled by $\sqrt[{}^{4}]{M/\left(2\times
10^{8}M_{\odot}\right)}$. The MOND boost to the Newtonian
$\sigma_{{}_{\text{LOS}}}$ can also be estimated by applying Equation 45 to
find $\overline{\nu}\equiv\overline{g}_{r}/g_{{}_{N,r}}$, the average factor
by which the Newtonian radial gravity is enhanced in MOND. The enhancement to
$\sigma_{{}_{\text{LOS}}}$ should then be $\approx\sqrt{\overline{\nu}}$, but
this approach has never been demonstrated.

MOND fares well with $\sigma_{{}_{\text{LOS}}}$ measurements of M31
satellites, where the EFE sometimes plays an important role (McGaugh &
Milgrom, 2013a, b). For the classical MW satellites, MOND agrees reasonably
well in the sense that the luminosity profile and global velocity dispersion
are mutually consistent (Sanders, 2021). For fainter satellites, care must be
taken to ensure that the analysed object would be tidally stable in a MOND
context (Section 5.1). Restricting attention to galaxies which should be
tidally stable in MOND, it agrees reasonably well with the observed
$\sigma_{{}_{\text{LOS}}}$ values of the classical MW satellites (Angus, 2008)
and dwarf spheroidal satellites in the LG more generally (McGaugh & Wolf,
2010). $N$-body simulations show that the slightly problematic MW satellite
Carina is in only mild tension with MOND because it requires
$M_{\star}/L_{V}=5.3-5.7$, somewhat higher than expected from stellar
population synthesis modelling (Angus et al., 2014).

In addition to the satellite galaxies of the MW and M31, the LG also contains
non-satellite dwarf galaxies that can be used to test MOND. Pawlowski &
McGaugh (2014a) predicted $\sigma_{{}_{\text{LOS}}}$ for the isolated LG
dwarfs Perseus I, Cetus, and Tucana to be 6.5 km/s, 8.2 km/s, and 5.5 km/s,
respectively, with an uncertainty close to 1 km/s in all cases due to the
uncertain $M_{\star}/L$ (see their section 3). Observationally,
$\sigma_{{}_{\text{LOS}}}$ of Perseus I is constrained to $4.2^{+3.6}_{-4.2}$
km/s, with a 90% confidence level upper limit of 10 km/s (Martin et al.,
2014). The reported $\sigma_{{}_{\text{LOS}}}$ for Cetus is $8.3\pm 1.0$ km/s
(Kirby et al., 2014) or $11.0^{+1.6}_{-1.3}$ km/s (Taibi et al., 2018). The
latest measurements for Tucana indicate that its
$\sigma_{{}_{\text{LOS}}}=6.2^{+1.6}_{-1.3}$ km/s (Taibi et al., 2020), with
the more careful analysis and larger sample size allowing the authors to rule
out earlier claims that $\sigma_{{}_{\text{LOS}}}>10$ km/s (Fraternali et al.,
2009; Gregory et al., 2019). In all three cases, there is good agreement with
the a priori MOND prediction, which is just Equation 56 as these dwarfs are
quite isolated (Pawlowski & McGaugh, 2014a).

A more recent study considered a much larger sample of LG dwarfs, though again
restricting to those which should be immune to tides (McGaugh et al., 2021).
Those authors defined a parameter $\beta_{\star}$ by which the
$\sigma_{{}_{\text{LOS}}}$ of a dwarf galaxy would need to be scaled up for it
to fall on the BTFR. In MOND, comparison of Equations 53 and 56 shows that we
expect $v_{{}_{f}}=\sigma_{{}_{\text{LOS}}}\sqrt[{}^{4}]{81/4}$ for an
isolated dwarf at low acceleration, which implies that scaling the velocity
dispersions by $\beta_{\star}=2.12$ should reconcile them with the BTFR. After
first establishing that the BTFR holds for 9 rotationally supported LG
galaxies not in the SPARC sample that underlies Figure 3, the authors then
obtained $\beta_{\star}$ observationally. The values scatter in a narrow range
around a median value of 2 if we assume that $M_{\star}/L_{V}=2$. The median
$\beta_{\star}$ reaches the MOND prediction of 2.12 if we assume a slightly
higher $M_{\star}/L_{V}=2.5$, which is quite reasonable for the typically old
dwarf galaxies in the LG. Alternatively, the velocity dispersions could be
higher than reported by 6% on average, which is also reasonable given the
typical uncertainties and the sample size. Therefore, the main conclusion of
McGaugh et al. (2021) was that MOND is able to account for the velocity
dispersions of LG dwarfs that should be immune to tides, which is necessary to
allow an equilibrium virial analysis. Since their sample generally avoids
dwarfs near the major LG galaxies, the EFE should not have significantly
influenced their results.

The application of Equation 56 to the galaxy known as NGC 1052-DF2 (hereafter
DF2) has received considerable attention because the observational estimate is
significantly lower (van Dokkum et al., 2018). However, their claim to have
falsified MOND was swiftly rebutted in a brief commentary on the work because
it did not consider the EFE from NGC 1052, which is at a projected separation
of only 80 kpc (Kroupa et al., 2018b) for a distance to both galaxies of 20
Mpc (Shen et al., 2021). In addition to this deficiency, several choices were
made by van Dokkum et al. (2018) which pushed the observational estimate of
$\sigma_{{}_{\text{LOS}}}$ to very low values, worsening the discrepancy with
Equation 56. One of the most serious problems unrelated to the
misunderstanding of MOND is that the statistical methods used to infer
$\sigma_{{}_{\text{LOS}}}$ were not well suited to the problem. Using instead
basic Gaussian statistics returns a higher $\sigma_{{}_{\text{LOS}}}$ (Martin
et al., 2018; Haghi et al., 2019a), though still below the result of Equation
56. Another issue with the work of van Dokkum et al. (2018) is their claim to
have discovered DF2 in an earlier work, even though it was clearly marked in
plate 1 of Fosbury et al. (1978) and had the alternative designation
[KKS2000]04. It is therefore clear that care must be taken before testing MOND
with the photometry and observed $\sigma_{{}_{\text{LOS}}}$ of a dwarf galaxy.

Figure 6: Results of numerical MOND simulations of DF2, showing
$\sigma_{{}_{\text{LOS}}}$ (top), half-mass radius $r_{h}$ (middle), and the
virial ratio $2K/\left|W\right|$ between twice the kinetic energy $K$ and the
magnitude of the virial
$W\equiv\sum_{i}m_{i}\left(\bm{r}_{i}\cdot\bm{g}_{i}\right)$ in the dwarf
reference frame (bottom), where $i$ labels different particles with mass
$m_{i}$ at position $\bm{r}_{i}$ experiencing gravity $\bm{g}_{i}$ in excess
of the external field from NGC 1052. The initial $r_{h}$ is 3 kpc (left) or 4
kpc (right). The vertical green shading in the top panels shows the range of
likely distances between DF2 and NGC 1052 given their projected separation is
80 kpc at a heliocentric distance of $\approx 20$ Mpc (Shen et al., 2021). The
horizontal shaded bands in the top panels show the observed
$\sigma_{{}_{\text{LOS}}}$ and its uncertainty based on ten globular clusters,
while the shaded region in the middle panels covers the range between their
$r_{h}$ and the slightly smaller stellar $r_{h}$. The dotted black line in
each top panel shows the semi-analytic prediction based on the fitting formula
illustrated in Figure 5. The other black lines show results of n-mody
simulations that consider the time-varying EFE but neglect tides on the dwarf,
which is started at a pericentre of 20 kpc and launched with the velocity
indicated in the middle left panel’s legend. The dashed red lines show por
models launched outwards from a pericentre of 40 kpc. The golden models in the
right panel show the most realistic por models started at 200 kpc and launched
inwards towards a pericentre of 64 kpc (dot-dashed) or 28 kpc (dashed), with
the latter providing a good fit to the observables. In all cases, only the
outgoing portion of the trajectory after pericentre passage is shown, so time
flows from left to right. Reproduced from figure 5 of Haghi et al. (2019a).

DF2 is very gas-poor (Chowdhury, 2019; Sardone et al., 2019), thereby
providing model-independent evidence that it feels a non-negligible external
gravitational field from a massive host galaxy. This is almost certainly the
nearby NGC 1052, whose projected separation of only 80 kpc implies an actual
separation of $\approx 100$ kpc. This makes the orbital period rather short,
causing the external field to be time-dependent on a scale not much longer
than the internal dynamical time of DF2. Moreover, tides on DF2 should in
principle also be considered as the virial theorem cannot be applied to an
object undergoing tidal disruption, as seems to be the case for NGC 1052-DF4
(Montes et al., 2020; Keim et al., 2022). These complications were handled
using fully self-consistent $N$-body simulations with por (Lüghausen et al.,
2015) in which DF2 was put on an orbit around NGC 1052 (Haghi et al., 2019a).
Those authors also conducted simulations with n-mody (Ciotti et al., 2006)
that include the EFE but not tides. The results of these simulations are shown
in Figure 6, which reproduces figure 5 of Haghi et al. (2019a). Their work
clarified that non-equilibrium memory effects would be quite small in the case
of DF2 for a wide range of plausible assumptions regarding its orbit around
NGC 1052 (compare the por and analytic results in the top panels). Such
effects might be more significant elsewhere, but DF2 can be approximated as
being in equilibrium with the present external field from NGC 1052. After
pericentre passage, $\sigma_{{}_{\text{LOS}}}$ of the dwarf takes time to rise
towards its equilibrium value. This memory effect (Haghi et al., 2009; Wu &
Kroupa, 2013) is partially counteracted by tides inflating
$\sigma_{{}_{\text{LOS}}}$ (compare por and n-mody results in Figure 6).
Before pericentre passage, both effects are expected to work in tandem to
inflate $\sigma_{{}_{\text{LOS}}}$ above the equilibrium for the local EFE.
The simulations of Haghi et al. (2019a) therefore support earlier analytic and
semi-analytic estimates that show the observed $\sigma_{{}_{\text{LOS}}}$ of
DF2 is well accounted for in MOND given its observed luminosity, size, and
position close to NGC 1052 (Famaey et al., 2018; Kroupa et al., 2018b). The
subsequently measured stellar body $\sigma_{{}_{\text{LOS}}}$ of
$8.5^{+2.3}_{-3.1}$ km/s (Danieli et al., 2019) is consistent with the
globular cluster-based estimate of $8.0^{+4.3}_{-3.0}$ km/s shown in Figure 6.

Equation 56 and its generalization to higher acceleration systems feeling a
non-negligible EFE (Figure 5) only pertain to the system as a whole. This is
suitable for unresolved observations, but in general we expect that
$\sigma_{{}_{\text{LOS}}}$ varies within a system. Assuming that it is
spherically symmetric, collisionless, and in equilibrium, it would satisfy the
Jeans equation

$\displaystyle\frac{d\left(\rho{\sigma_{r}}^{2}\right)}{dr}+\frac{2\beta\rho{\sigma_{r}}^{2}}{r}\,=\,-\rho
g\,,\quad\text{where}~{}\beta~{}\equiv~{}1-\frac{{\sigma_{t}}^{2}}{{2\sigma_{r}}^{2}}\,,$
(60)

$\rho$ is the tracer density, $r$ is the radius, $\sigma_{r}$ is the velocity
dispersion in the radial direction, $\sigma_{t}$ is the total velocity
dispersion in the two orthogonal (tangential) directions (hence the factor of
2), $\beta$ is the anisotropy parameter, and $g>0$ is the inward gravity. The
results are insensitive to an arbitrary rescaling of $\rho$.

Similarly to RCs, we expect the velocity dispersion profile to flatten at
large radii in MOND rather than continue on a Keplerian decline, as would
occur in Newtonian gravity. Thus, the flattened outer velocity dispersion
profiles of some Galactic globular clusters were used to argue in favour of
MOND (Hernandez et al., 2013; Durazo et al., 2017). However, the outskirts of
globular clusters would contain stars that are no longer gravitationally bound
to the cluster but continue to follow its Galactocentric orbit. These unbound
‘potential escaper’ stars can mimic an outer flattening of the velocity
dispersion profile (Claydon et al., 2017), so MOND might not be the only
explanation. Moreover, we do not expect significant MOND effects in a nearby
globular cluster due to the Galactic EFE (Section 2.4).

A more distant globular cluster would be less affected by the EFE. If in
addition its internal acceleration is small, then the cluster could serve as
an important test of MOND (Haghi et al., 2009; Haghi et al., 2011). One such
example is NGC 2419, which is 87.5 kpc away (Di Criscienzo et al., 2011).
Analysis of its internal kinematics led to the conclusion that it poses severe
problems for MOND (Ibata et al., 2011a, b). However, it was later shown that
allowing a radially varying polytropic equation of state yields quite good
agreement with the observed luminosity and $\sigma_{{}_{\text{LOS}}}$ profiles
(Sanders, 2012a, b). More importantly, it is very misleading to consider the
formal uncertainties alone because there is always the possibility of small
systematic uncertainties such as rotation within the sky plane.

Another system with known $\sigma_{{}_{\text{LOS}}}$ profile is the ultra-
diffuse dwarf galaxy DF44. Its global $\sigma_{{}_{\text{LOS}}}$ was claimed
to exceed the MOND prediction (Hodson & Zhao, 2017b) based on earlier
observations (van Dokkum et al., 2016). However, subsequent observations
reduced the estimated $\sigma_{{}_{\text{LOS}}}$ and revealed its radial
profile (van Dokkum et al., 2019). This is reasonably consistent with MOND
expectations (Bílek et al., 2019). Joint modelling of its internal dynamics
and its star formation history in the IGIMF context showed DF44 to be in mild
$2.40\sigma$ tension with MOND (Haghi et al., 2019b). Its internal
acceleration of $\approx 0.2\,a_{{}_{0}}$ (see their section 5) makes this a
non-trivial success, though still somewhat subject to the mass-anisotropy
degeneracy (Equation 60).

The agreement of MOND with both the low $\sigma_{{}_{\text{LOS}}}$ of DF2 and
the high $\sigma_{{}_{\text{LOS}}}$ of DF44 despite a similar baryonic content
in both cases is due to the EFE playing an important role in DF2 but not in
DF44. The stronger external field on DF2 can be deduced independently of MOND
based on comparing their environments. We therefore see that dwarf spheroidal
galaxies can provide strong tests of MOND due to their weak internal gravity,
but this makes them more susceptible to the EFE and to tides. The lower
surface density also makes accurate observations more challenging, but it
should eventually be possible to test the predicted velocity dispersions of
e.g. the 22 ultra-diffuse dwarf galaxies considered in Müller et al. (2019).
Note that since they used Equation 57, the predictions are likely slight
underestimates if the EFE is important (Figure 5), which applies to the vast
majority of their sample.

Figure 7: Left: The RAR defined by disc galaxies and strong gravitational
lenses (Section 3.4; see also Tian & Ko, 2017), shown here on logarithmic axes
as the relation between $g_{{}_{N}}$ of the baryons and the amount by which
the observationally inferred gravity exceeds this (solid black points with
uncertainties). The fit to the data is shown as a solid black line
(subtracting 1 from Equation 55). The light purple squares show results from
500 Monte Carlo trials of randomly selected points at radii of
$\left(0.01-2\right)$ effective radii in galaxies with a virial mass of
$\left(10^{11}-10^{12}\right)M_{\odot}$. A Burkert profile is assumed for the
halo, with properties drawn from a $\Lambda$CDM distribution. The baryonic
component is assigned using abundance matching (Behroozi et al., 2013). The
dotted curves show the radial run for three individual galaxies. The binned
Monte Carlo results (purple squares with error bars) indicate that a
significant contribution from the halo is expected at high $g_{{}_{N}}$, at
odds with observations. This discrepancy is more severe for the NFW profile
(not shown here, but shown in original publication). Reproduced from figure 3
of Tian & Ko (2019). Right: Similar to the left panel, but now with the
$\Lambda$CDM results (Equation 61; solid red line) obtained from a different
semi-analytic prescription (Paranjape & Sheth, 2021). The dotted red lines
show their estimated intrinsic dispersion of 0.2 dex, which corresponds to a
smaller uncertainty of 0.075 dex in the total $g$ (see their section 4.1). The
mean simulated relation is known far more precisely due to the large number of
mock galaxies they considered. The mean observed relation (black points with
error bars) is from table 1 of Tian & Ko (2019).

Observations are easier for a massive elliptical galaxy, but the central
region tends to be in the Newtonian regime. Even so, MOND does still affect
the overall size of the system. This is because if we assume that it is close
to isothermal and in the Newtonian regime, then we run into the problem that a
Newtonian isothermal sphere has a divergent mass distribution (Binney &
Tremaine, 2008), which is also the case for the NFW profile expected in
$\Lambda$CDM (Navarro, Frenk & White, 1997). However, isothermal spheres in
MOND have a finite total mass (Milgrom, 1984), as do all polytropes in the
deep-MOND regime (Milgrom, 2021a). Thus, a nearly isothermal system that
formed with a radius initially much smaller than its MOND radius (Equation 22)
should expand until it reaches its MOND radius. At that point, further
expansion would be difficult due to the change in the gravity law. This
argument does not apply to low-density systems initially larger than their
MOND radius, but we might expect that a wide variety of systems were initially
smaller. These would end up on a well-defined mass-radius relation. Elliptical
galaxies do seem to follow such a relation, which perhaps explains why their
internal accelerations only exceed $a_{{}_{0}}$ by at most a factor of order
unity (Sanders, 2000). Moreover, the DM halos of elliptical galaxies inferred
in a Newtonian gravity context have rather similar scaling relations to the
PDM halos predicted by MOND (Richtler et al., 2011).

The stars in an elliptical galaxy can serve as tracers of its potential even
if they cannot be resolved individually. This has led to observational
projects such as ATLAS3D (Cappellari et al., 2013) and SDSS-IV MaNGA (Bundy et
al., 2015), an extension of the original Sloan Digital Sky Survey (SDSS; SDSS
Collaboration, 2000). The so-obtained stellar $\sigma_{{}_{\text{LOS}}}$
profiles of 19 galaxies reveal a characteristic acceleration scale consistent
with $a_{{}_{0}}$ as determined using spirals (Chae et al., 2020a). While the
uncertainties are larger and their study consists of 387 data points rather
than the 2693 in SPARC (Lelli et al., 2016), it is still important to note
that MOND works fairly well in elliptical galaxies.

This success may seem trivial due to the high accelerations. However, there is
no reason for physical DM particles to avoid high-acceleration regions. As a
result, it is entirely possible that the central regions of elliptical
galaxies should contain more DM in the $\Lambda$CDM picture than PDM in the
MOND picture. Indeed, the rather small amounts of missing gravity in the
central regions of elliptical galaxies are actually problematic for
$\Lambda$CDM, but this can be explained naturally in MOND if the interpolating
function is chosen appropriately (Tian & Ko, 2019). This is shown in their
figure 3, which we reproduce here as the left panel of our Figure 7. When
$g_{{}_{N}}\approx 10\,a_{{}_{0}}$, the predicted extra gravity from the DM
component is too large to sit comfortably with the data if a Burkert profile
is assumed for the halo. The other panels in figure 3 of Tian & Ko (2019) show
that the discrepancy is worse for the NFW profile predicted in $\Lambda$CDM
(Navarro et al., 1997). In MOND, it is possible to accommodate the
observations with a sufficiently sharp transition from the Milgromian to the
Newtonian regime at these high accelerations, e.g. with an exponential cutoff
to the MOND corrections (Equation 55).

That $\Lambda$CDM yields too much DM in the central regions of massive
galaxies is also evident from a different semi-analytic prescription to assign
a baryonic component to purely DM halos (Paranjape & Sheth, 2021). As shown in
their figure 2 and stated in their equation 19, their simulated RAR is well
fit by the interpolating function

$\displaystyle\nu~{}=~{}\left[\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{1}{y^{0.8}}}\right]^{1/0.8}\,,\quad
y~{}\equiv~{}\frac{g_{{}_{N}}}{a_{{}_{0}}}\,.$ (61)

We use the right panel of Figure 7 to show this as a relation between
$g_{{}_{N}}$ and the halo contribution $g_{{}_{h}}\equiv g-g_{{}_{N}}$. This
allows a comparison with the binned observational results from table 1 of Tian
& Ko (2019). The dotted red lines show the estimated intrinsic dispersion of
0.2 dex in the simulated $g_{{}_{h}}$ (section 4.1 of Paranjape & Sheth,
2021). The uncertainty in the mean simulated relation is much smaller due to a
very large sample size. As a result, the simulated RAR is discrepant with
observations at the high-acceleration end. In particular, the extra halo
contribution must peak at $g_{{}_{N}}\approx 3\,a_{{}_{0}}$ and start dropping
thereafter, but the simulation results indicate that $g_{{}_{h}}$ continues
rising with $g_{{}_{N}}$. This is essentially the same result as shown in the
left panel of Figure 7, but with a different “baryonification” scheme
(Paranjape et al., 2021). Another issue is that the intrinsic dispersion in
the total $g$ at fixed $g_{{}_{N}}$ is expected to be 0.075 dex (Paranjape &
Sheth, 2021), but observationally “the intrinsic scatter in the RAR must be
smaller than 0.057 dex” (section 4.1 of Li et al., 2018). It is therefore
clear that the results obtained by Paranjape & Sheth (2021) are actually very
problematic for $\Lambda$CDM, contrary to what is stated in their abstract.
The main reason seems to be that they only considered data on galaxy
kinematics (Chae et al., 2019), but strong lenses also provide important
constraints, especially at the high-acceleration end (Section 3.4; see also
section 2 of Tian & Ko, 2019). The discrepancy is related to the expected
adiabatic compression of CDM halos due to the gravity from a centrally
concentrated baryonic component (Li et al., 2022; Forouhar Moreno et al.,
2022). If $M_{\star}/L$ is higher in elliptical galaxies than typically
assumed due to a higher proportion of stellar remnants (as expected with the
IGIMF theory), then $g_{{}_{h}}$ would have to be smaller still, worsening the
discrepancy with $\Lambda$CDM expectations.

Planetary nebulae can also serve as bright tracers that should have similar
kinematics to the stars. This technique was used in three intermediate
luminosity elliptical galaxies, revealing an unexpected lack of DM when
analysed with Newtonian gravity (Romanowsky et al., 2003). These observations
can be naturally understood in MOND due to the high acceleration (Milgrom &
Sanders, 2003). The problem was revisited more recently with better data on
those three galaxies and newly acquired data on four more galaxies (Tian & Ko,
2016). Their analysis confirmed the earlier finding that MOND provides a good
description of the observed kinematics.

Though fewer in number than stars, satellite galaxies must also trace the
gravitational field modulo the mass-anisotropy degeneracy (Equation 60). This
technique was attempted by Angus (2008), who found that it is possible to fit
the observations of Klypin & Prada (2009) if the velocity dispersion tensor
transitions from tangentially biased ($\beta<0$) in the central regions to
radially biased ($\beta>0$) in the outskirts. This was expected a priori from
MOND simulations of dissipationless collapse (Nipoti et al., 2007). However,
the mass-anisotropy degeneracy and possible LOS contamination mean that this
is not a very sensitive test of the gravity law.

Figure 8: The elliptical galaxy RAR defined by the relation between
$g_{{}_{N}}$ and the gravity required for hydrostatic equilibrium of the X-ray
halo (Equation 60 with $\beta=0$). Results are shown as filled circles, with
the colour indicating the galaxy (see the legend). The RAR of disc galaxies
(Figure 3) is indicated by the shade of blue in the background, with the fit
to it (Equation 55; solid black line) and the width of the distribution at
fixed $g_{{}_{N}}$ (dotted red lines) shown similarly to that figure. Notice
that disc galaxies define the same RAR as X-ray halos around elliptical
galaxies. Reproduced from figure 8 of Lelli et al. (2017b) by permission of
Federico Lelli and the American Astronomical Society.

The satellite region can sometimes be probed using an X-ray halo, whose
temperature and density profile allow for a determination of the gravitational
field assuming hydrostatic equilibrium. The uncertainty is reduced somewhat as
we expect $\beta=0$ for collisional gas. This technique has been applied to
NGC 720 and NGC 1521, yielding good agreement with MOND over the acceleration
range $\left(0.1-10\right)a_{{}_{0}}$ (Milgrom, 2012b). The fit is less good
in the central few kpc of NGC 720, but it is quite likely that the assumption
of hydrostatic equilibrium breaks down here (Diehl & Statler, 2007), possibly
due to feedback from stars and active galactic nuclei. Additional studies of
X-ray halos around elliptical galaxies are reviewed in Buote & Humphrey
(2012), with the results for 9 well-observed galaxies summarized in figure 8
of Lelli et al. (2017b). We reproduce this as our Figure 8, showing how the
results fall on the RAR traced by the RCs of spiral galaxies. The RAR is thus
not unique to rotationally supported systems or to rotational motion, nor is
it confined to thin disc galaxies.

####  3.2.2 Rotation of a sub-dominant component

The mass-anisotropy degeneracy can be broken by considering a thin rotating
gas disc. By definition, this can only be a sub-dominant component of an
elliptical galaxy. Even so, such a disc would provide an excellent way to
measure $g_{r}$. This is possible in 16 early-type galaxies from the ATLAS3D
survey (den Heijer et al., 2015). Using 21 cm HI observations extending beyond
5 effective radii, those authors showed that the analysed galaxies very
closely follow the BTFR defined by spiral galaxies. The BTFR was again the
most fundamental relation (smallest scatter), not e.g. the relation between
$M_{\star}$ and $v_{{}_{f}}$. The results for these 16 galaxies are shown on
an RAR diagram in figure 8 of Lelli et al. (2017b), revealing good agreement
with the RAR defined by disc galaxies.

More recently, HI RCs have been obtained for three lenticular galaxies
extending out to $10-20$ effective radii (Shelest & Lelli, 2020). Those
authors used 3.6 $\mu$m photometry to reduce uncertainty on $M_{\star}/L$.
Their figure 8 (reproduced here as our Figure 9) shows that the 3 analysed
galaxies fall on the spiral RAR, with the dataset probing $g_{{}_{N}}$ out to
$\approx 1.5$ dex either side of $a_{{}_{0}}$ (the range in $g$ is slightly
smaller).

Figure 9: The RAR defined by the gas RC in three lenticular galaxies (points
with error bars), as determined from CO in the inner regions and deep HI
observations in the outer regions. The colour of each data point indicates
which galaxy it comes from (see the legend). Notice the similarity to the mean
RAR of disc galaxies (solid black line) and its standard deviation (dashed
black lines). Reproduced from figure 8 of Shelest & Lelli (2020) with
permission from Astronomy & Astrophysics.

Perhaps because these results are still very recent, not much attention has
been paid to whether $\Lambda$CDM can get the same tight RAR in elliptical and
lenticular galaxies as that observed in spirals. Since the star formation
histories are likely to have been different, the proportion of baryons lost by
feedback should also be different. This makes it difficult to understand how
the same tight RAR arises in galaxies of different Hubble types. However, this
was the a priori prediction of MOND.

###  3.3 Observational signatures of the EFE

In this section, we review the observational evidence for the EFE, which
theoretically is an integral part of MOND (Section 2.4). A classic example of
the EFE is when a satellite galaxy orbits around its host. The EFE was
included in an early analysis of $\sigma_{{}_{\text{LOS}}}$ in 7 Galactic
satellites, some of which are dominated by the EFE (Milgrom, 1995). Neglecting
the EFE would enhance the self-gravity and thus reduce the MOND dynamical
$M/L$, which could be problematic in the case of Fornax because the required
$M/L$ is already on the low side.

The EFE should also be present in satellites of M31, which is near enough that
some fairly detailed observations are now available. McGaugh & Milgrom (2013b)
considered the dynamics of three matched pairs of satellites which are
photometrically similar but differ in their distance from M31 or in their
size, leading to the EFE being less important in one satellite than for the
other (see their section 2.4). Since the considered dwarfs are deep in the
MOND regime, neglecting the EFE would mean that $\sigma_{{}_{\text{LOS}}}$
depends only on the total luminosity, so it should be the same for both
galaxies in each matched pair (Equation 56). This is no longer the case once
the EFE is considered. The difference between the predicted
$\sigma_{{}_{\text{LOS}}}$ with or without the EFE was too small to be tested
by observations in one matched pair, but for the other two pairs, the data
seemed to indicate that the satellite more dominated by the EFE does indeed
have a lower $\sigma_{{}_{\text{LOS}}}$.

Another strong hint for the EFE comes from the internal dynamics of the MW
satellite Crater II at a Galactocentric distance of $d\approx 120$ kpc
(Torrealba et al., 2016). Neglecting the Galactic EFE and assuming that Crater
II has $M_{\star}/L_{V}=2$ with this being uncertain by a factor of 2, its
predicted $\sigma_{{}_{\text{LOS}}}^{\text{vir}}=4.0_{-0.6}^{+0.8}$ km/s, but
this decreases to $2.1_{-0.6}^{+0.9}$ km/s once the EFE is included (McGaugh,
2016). The ‘vir’ superscript indicates that their calculation neglects the
role of tides, which can be important for an ultra-diffuse galaxy such as
Crater II because its half-mass radius $r_{h}=1.07\pm 0.08$ kpc (table 4 of
Torrealba et al., 2016) despite a very low $M_{s}=3.2\times 10^{5}\,M_{\odot}$
for $M_{\star}/L_{V}=2$. Tides (discussed further in Section 5.1) can
generally be expected to inflate $\sigma_{{}_{\text{LOS}}}$ further (figure 4
of Brada & Milgrom, 2000b). The amount by which this occurs can be estimated
by considering the opposite limit in which self-gravity is negligible and the
observed stars are considered to be travelling along independent
Galactocentric orbits. $\sigma_{{}_{\text{LOS}}}$ would then be set by the
Galactic RC and the range of orbital phases that we observe, which is roughly
given by the angular size $r_{h}/d$. Using this to scale the Galactic
$v_{c}\approx 200$ km/s, we get that
$\sigma_{{}_{\text{LOS}}}^{\text{tid}}\approx v_{c}r_{h}/d=1.8$ km/s, though
this would vary somewhat around the orbit given the high eccentricity (Li et
al., 2021b). If we assume that the actual $\sigma_{{}_{\text{LOS}}}$ should be
found by adding $\sigma_{{}_{\text{LOS}}}^{\text{vir}}$ and
$\sigma_{{}_{\text{LOS}}}^{\text{tid}}$ in quadrature, we get that
$\sigma_{{}_{\text{LOS}}}$ should be at least 3.9 km/s if the Galactic EFE has
no effect on the self-gravity of a satellite. In this case, our simple
estimate shows that tides would be relatively unimportant due to the strong
self-gravity of Crater II. However, including the EFE leads to its predicted
$\sigma_{{}_{\text{LOS}}}$ dropping to 2.8 km/s. Observationally, Crater II
has $\sigma_{{}_{\text{LOS}}}^{\text{obs}}=2.7\pm 0.3$ km/s (Caldwell et al.,
2017). This shows that in MOND, Crater II must feel the Galactic EFE. However,
its internal kinematics are not a strong test of MOND because we can achieve
plausible agreement with observations even if we neglect self-gravity
altogether. The observations can thus be considered as putting a useful upper
limit on the amount of self-gravity in Crater II.

Detailed predictions for the internal kinematics of dwarfs are difficult to
obtain in $\Lambda$CDM, but it should still be possible to statistically
predict roughly what $\sigma_{{}_{\text{LOS}}}$ Crater II should have. McGaugh
(2016) predicted in their section 2 that $\sigma_{{}_{\text{LOS}}}=17.5$ km/s
based on an empirical scaling relation (Walker et al., 2010). The observed
properties of Crater II remain difficult to reconcile with $\Lambda$CDM even
if significant tidal disruption is considered, at least if it formed inside a
CDM halo (Borukhovetskaya et al., 2022; Errani et al., 2022). This is strongly
suggested by the fact that its Newtonian dynamical $M_{\star}/L_{V}\approx 50$
(Caldwell et al., 2017; Fu et al., 2019). If instead Crater II is purely
baryonic like a globular cluster, then its large size and high observed
$\sigma_{{}_{\text{LOS}}}$ imply that it must now be past perigalacticon.
Proper motion data indicates that this would have been at $33\pm 8$ kpc (table
4 of Li et al., 2021b), which combined with the very large mass ratio between
Crater II and a CDM-rich MW would have caused complete tidal destruction at
that point (Equation 76). This problem is somewhat alleviated in MOND because
a purely baryonic Crater II is much more resilient to tides if gravity is
Milgromian rather than Newtonian, making it more plausible that a very diffuse
remnant would still be recognisable.

These issues with Crater II highlight that it is generally difficult to
identify satellite galaxies that might be significantly affected by the EFE
but not by tides, thus remaining amenable to an equilibrium virial analysis. A
more promising approach might be to exploit the differing levels of tidal
stability in $\Lambda$CDM and MOND by comparing with some model-independent
measure of whether a galaxy is affected by tides. This would test the EFE
because the tidal radius is necessarily in the EFE-dominated regime. We
explore such ideas further in Section 5.1. Another possibility is to construct
detailed numerical simulations of tidal streams and compare with observations
in order to search for specific signatures of the EFE, in particular an
asymmetry between the leading and trailing arms (Section 5.2; see also Thomas
et al., 2018).

Turning now to disc galaxies, the EFE would cause the outer RC to decline,
even though a flat RC is expected for an isolated system. Observational
evidence for this was found by Haghi et al. (2016). In a major breakthrough,
Chae et al. (2020b) extended their work by cross-correlating the high-quality
SPARC dataset with the large-scale gravitational field (Desmond et al., 2018).
The analysis of Chae et al. (2020b) considered whether including the EFE leads
to better RC fits, bearing in mind that allowing an extra model parameter
inevitably improves the fit to some extent. Even so, they found very strong
evidence for the EFE. This is clearly evident in their figure 2, which reveals
the expected flat outer RC for relatively isolated galaxies, but galaxies that
should experience a stronger EFE clearly have a declining outer RC. These
observations cannot easily be fit by instead altering conventional parameters
like the distance and inclination, leading Chae et al. (2020b) to conclude
that they detected the EFE at $8\sigma-11\sigma$ in some galaxies but not at
all in others. The EFE strength inferred from the RC of a galaxy in a more
isolated environment is typically much less than for a galaxy in a more
crowded environment (Chae et al., 2021). These pioneering studies test for the
first time whether the strong equivalence principle at the heart of GR remains
valid at the low accelerations typical of galactic outskirts. A significant
failure of the principle could be analogous to the breakdown of energy
equipartition at low temperatures in condensed matter physics (Section 2.2).
Departures from GR are inevitable in some cases, so empirically finding
precisely which ones could provide important clues to a better theory of
gravity that incorporates quantum effects.

It is interesting to consider whether $\Lambda$CDM predicts a correlation
between the RC of a galaxy and the external field on its barycentre from
large-scale structure. While this cannot arise as a fundamental consequence of
the gravity law, we might expect that in a high-density region, a CDM halo
would typically have a higher density by virtue of forming earlier. At fixed
$g_{{}_{N}}$, this would make the observed $g$ higher, a result which was
recently demonstrated using semi-analytic models (Paranjape & Sheth, 2022).
The predicted sign of the EFE is thus opposite to observations (Chae et al.,
2020b, 2021).

In MOND, an important consequence of the EFE is that the point mass
gravitational field transitions to an inverse square law once
$\bm{g}_{{}_{e}}$ dominates (Equation 44). Without the EFE, Equation 22 would
apply, creating a logarithmically divergent potential (e.g. Equation 23).
Regardless of theoretical issues, an observational consequence of this would
be an extremely deep Galactic potential from which escape is virtually
impossible. Including even a weak EFE would substantially alter this picture.
These scenarios can be distinguished using the escape velocity from the Solar
circle of the MW (Famaey et al., 2007; Wu et al., 2008) $-$ and indeed the
escape velocity curve more generally (Banik & Zhao, 2018b). This agrees quite
well with observations near the Solar circle (Monari et al., 2018) and over
Galactocentric radii of $8-50$ kpc (Williams et al., 2017), which is another
indication that the EFE from more distant masses like M31 does affect the
internal dynamics of the MW in a MOND context.

The EFE can have important effects even in situations where it is sub-dominant
to the internal gravity. For instance, the central regions of M33 would not be
much affected by the EFE from M31, but even so a hydrodynamical MOND
simulation of M33 agrees better with some observables once the EFE is
approximately included (Banik et al., 2020). A similar situation is apparent
in AGC 114905 (Banik et al., 2022b). In addition to weakening the self-gravity
of a system, the EFE also causes the point mass potential to become
anisotropic (Equation 41). If the internal and external gravity are
comparable, the potential can even become asymmetric with respect to
$\pm\bm{g}_{{}_{e}}$. This can leave a characteristic imprint on tidal
streams, as discussed further in Section 5.2. The asymmetry could also cause
discs to become warped in an external field (Banik et al., 2020). Indeed, the
warp of the Galactic disc might be due to the EFE from the LMC (Brada &
Milgrom, 2000a).

Figure 10: Histogram of the sky-projected relative velocity of wide binary
stars in the Solar neighbourhood divided by the predicted Keplerian velocity
at their sky-projected separation $r_{p}$, shown for different $r_{p}$ bins as
indicated at the top of each panel. Observational results from Gaia data
release 2 (red) are compared with the prediction of Newtonian gravity (blue)
and MOND (magenta) for $g_{{}_{N,e}}=1.2\,a_{{}_{0}}$ as appropriate for the
Solar neighbourhood. The green distribution is for a toy version of MOND which
neglects the Galactic EFE. The peak of this distribution shifts to the right
at larger $r_{p}$ (bottom panel) because the binary circular velocity no
longer declines with separation, analogously to the flat RCs of isolated
galaxies. The predicted peak shift is strongly excluded by the data. Note that
MOND is not excluded here as the EFE is an integral part of it. Including the
EFE considerably narrows the gap between the Newtonian and Milgromian
predictions (Section 2.4), causing both to be in reasonable agreement with
Gaia data. Adapted from figure 11 of Pittordis & Sutherland (2019) by its lead
author.

Besides the open cluster data discussed in Milgrom (1983a), the properties of
Ly-$\alpha$ absorbers also require the EFE in a MOND context (Section 8; see
also Aguirre et al., 2001). On large scales, the MOND scenario for the KBC
void and Hubble tension requires the EFE to make the void dynamics consistent
with observations (Section 9.2.2; see also the inference on the EFE in figure
4 of Haslbauer et al., 2020).

Similarly to Galactic star clusters in the Solar neighbourhood of the MW,
local wide binary stars would also experience too much self-gravity in MOND if
the Galactic EFE is neglected. Wide binaries provide a promising test of MOND
(Section 11.3). By definition, these systems have low internal accelerations,
which would cause a large departure from Newtonian expectations if the
binaries were isolated (Hernandez et al., 2012). However, the difference is
much smaller once the Galactic EFE is considered (Banik & Zhao, 2018d). The
sky-projected separation $r_{\text{sky}}$ of each binary is very small
compared to its Galactocentric radius, so considering the binary in isolation
would indeed involve neglecting the Galactic EFE for a very wide range of
plausible choices on where to put the system’s boundary. Pittordis &
Sutherland (2019) showed that Gaia data are strongly inconsistent with the
large ‘MOND’ signal that would be expected if we neglect the Galactic EFE,
especially for the more widely separated binaries where the internal
accelerations are lower (see their figure 11, reproduced here as our Figure
10). Including the EFE makes the observations consistent with MOND, though
they also agree with Newtonian expectations. There have also been many
terrestrial experiments at low internal accelerations that show no sign of
departure from Newtonian dynamics (Gundlach et al., 2007). We therefore
conclude that there are compelling theoretical and empirical reasons for
supposing that the EFE is an important aspect of MOND $-$ and probably of the
real Universe as well.

###  3.4 Strong gravitational lensing

A gravitational lens is said to be strong if it leads to multiple resolved
images of the same background source. This definition is somewhat technology-
dependent because gravitational microlensing events (Refsdal, 1966) should
yield multiple images which could in principle be resolved with future
technology. Nonetheless, the definition is useful because important
constraints can be obtained from the angular separation between the multiple
images in a strong gravitational lens.

Due to the cosmic coincidence of MOND (Equation 16), photons involved in
strong lensing are necessarily in the Newtonian regime at the point of closest
approach if the lens and source are at cosmological distances (Sanders, 2014).
This is not true for the entire trajectory, so some MOND enhancement is
expected. Moreover, we expect some departure from Newtonian expectations even
if $g>a_{{}_{0}}$ because the RAR (Figure 3) implies a fairly gradual
transition between the Newtonian and Milgromian regimes. The extra light
deflection due to the enhanced MOND gravity is expected to be $\approx 15\%$
for the simple interpolating function (Equation 19).

The best-known examples of gravitational lenses are Einstein rings, where the
lens and unlensed source are at almost the same location on our sky. This
makes the problem nearly axisymmetric, leading to a wide range of possible
photon paths that reach the Earth. If the lens is a massive galaxy, we can
constrain its potential using the angular size of the Einstein ring and the
redshifts of the source and lens, bearing in mind that the deflection angle
and non-relativistic $\bm{g}$ are related in the same way as in GR (Section
2.6). Due to the typically cosmological distances, it is necessary to assume a
background cosmology, with the choice generally being a standard flat
background cosmology because this should be appropriate to both $\Lambda$CDM
and MOND (Section 9). If the redshift is not too large, the angular diameter
distance can be determined with little dependence on cosmology provided the
Hubble constant is known.

MOND was shown to be consistent with 10 lenses in the Centre for Astrophysics-
Arizona Space Telescope Lens (CASTLES; Lehár et al., 2000) catalogue of strong
lenses (Chiu et al., 2011). A much larger sample was provided by the Sloan
Lens Advanced Camera for Surveys (SLACS) catalogue (Bolton et al., 2006; Auger
et al., 2009). MOND provides a good fit to all 36 (Sanders & Land, 2008), 65
(Sanders, 2014), or 57 (Tian & Ko, 2017) lenses analysed, with the studies
adopting different quality cuts. MOND has also been applied to the time delays
in three gravitational lenses (Tian et al., 2013), yielding broad agreement
for a standard background cosmology. As with the deflection angle, the MONDian
Shapiro delay (Shapiro, 1964) is expected to have a standard dependence on the
non-relativistic $\bm{g}$.

###  3.5 Weak gravitational lensing

The surface density rarely reaches the threshold required for strong
gravitational lensing. Nonetheless, any source is always distorted to some
extent by a foreground lens. This regime is known as weak lensing. Since the
true shape of the source is rarely known, it is necessary to use statistical
techniques to infer how much distortion has been induced by the lens, usually
by stacking many sources around many lenses. The idea is that the source
(typically a thin disc galaxy) would have its sky-projected major axis
oriented in a random direction. In contrast, a lens induces a subtle spacetime
distortion that causes the major axes of background galaxies to preferentially
align in a certain direction, which for a point mass lens would be the
tangential direction. One can imagine that the lens tries to distort
background galaxies into an Einstein ring, but is insufficiently powerful to
do so.

Weak lensing effects arise because the direction and magnitude of the
deflection angle vary over different parts of the source. In the deep-MOND
limit around an isolated point mass $M$, the deflection angle in the weak
deflection limit follows from Equations 3 and 50.

$\displaystyle\Delta\widehat{\bm{n}}~{}=~{}-\frac{2\mathrm{\pi}\sqrt{GMa_{{}_{0}}}}{c^{2}}\widehat{\bm{r}}_{p}~{}=~{}-2\mathrm{\pi}\left(\frac{v_{{}_{f}}}{c}\right)^{2}\widehat{\bm{r}}_{p}\,,$
(62)

where $\widehat{\bm{r}}_{p}$ is the direction from the lensing mass to the
photon at closest approach. We used Equation 53 to write the result in terms
of $v_{{}_{f}}$, which might be more readily observable. Unlike in the
Newtonian case, the deep-MOND deflection angle is independent of the impact
parameter. As a result, shear arises only due to changes in the direction of
the deflection caused by changes in $\widehat{\bm{r}}_{p}$, or equivalently
because the radial direction is in general different.

Galaxy-galaxy weak lensing was detected at high significance in a stacked
analysis of 12 million sources (Brimioulle et al., 2013). Their results agree
quite well with the MOND prediction (Equation 62) if the blue lens subsample
(presumably spiral galaxies) has $M_{\star}/L_{B}$ in the range $1-3$, which
covers the range expected from stellar population synthesis modelling
(Milgrom, 2013). The red lens subsample (presumably elliptical galaxies)
requires a higher $M_{\star}/L_{B}$ of $3-6$, in line with expectations for an
older population. The data cover out to projected separations of $\approx 280$
kpc, thus reaching accelerations down to only a few percent of $a_{{}_{0}}$.
This is close to the limit beyond which the assumption of isolation is
expected to break down due to the EFE from large-scale structure. The
consistency of MOND with galaxy-galaxy weak lensing was again demonstrated
using 33,613 isolated central galaxies (Brouwer et al., 2017). Recently, an
even larger sample was used to obtain similar results (Brouwer et al., 2021).
We reproduce their figure 4 in our Figure 11, showing that different surveys
and sample selections give fairly similar results. The very significant amount
of missing gravity evident at $g_{{}_{N}}\ll a_{{}_{0}}$ is correctly
reproduced in MOND down to the limit of the data. At the very low acceleration
end, there is a hint of a downturn from the RC-based RAR that might be due to
the EFE.

Figure 11: The RAR determined from stacked galaxy-galaxy weak gravitational
lensing using Equation 50 (points with error bars). Different surveys and
sample selections are used, as described in the legend. The relation for disc
galaxy RCs is shown towards the upper right using blue shading and red squares
for the binned results (similarly to Figure 3), with the fit (Equation 55)
shown as a solid grey line. The dashed grey curve shows the estimated
deviation due to the EFE assuming $g_{{}_{N,e}}=0.003\,a_{{}_{0}}$, reasonable
for the fairly isolated galaxies considered here. In MOND, the EFE is expected
to affect the results in the blue shaded region on the left (cf. Figure 1).
Reproduced from figure 4 of Brouwer et al. (2021) with permission from
Astronomy & Astrophysics.

It has been claimed that the weak lensing signal is anisotropic relative to
the lens, which in a Newtonian analysis implies a DM halo flattened similarly
to the baryonic disc (Schrabback et al., 2021b). The anisotropic signal has
not been conclusively detected, but there is strong evidence for it over the
radial range $45-200$ kpc. Since the disc is typically much smaller, this may
pose problems for the MOND scenario. However, there could be an extended halo
of hot gas, which we might expect would be flattened along the disc spin axis.
Moreover, the EFE from large-scale structure can also induce anisotropy in the
potential, even for a point mass lens (Equation 41; see also Oria et al.,
2021). Since the orientation of the central disc is likely to correlate with
large-scale structure due to the EFE and other effects, it is possible that
the PDM halo of a galaxy is typically flattened along the same axis as its
disc. Weak lensing is also sensitive to other structures along the LOS (Feix
et al., 2008), so the signal might not be directly associated with the lens
galaxy. Further work would help to clarify if the results are genuinely
problematic for MOND, though this is likely to require a cosmological
hydrodynamical simulation to assess the large-scale structure and the
possibility of a hot gas halo.

###  3.6 Implications for alternatives to $\Lambda$CDM and MOND

The results discussed thus far broadly favour MOND over $\Lambda$CDM, but it
is helpful to consider other solutions to the missing gravity problem. These
are generally much less theoretically mature, partly due to being proposed
more recently. Meaningful predictions can still be obtained in some systems,
so we briefly consider such ‘third ways’ in this section.

In the emergent gravity (EG) theory, gravity is an emergent phenomenon arising
from some currently not understood microphysical processes (Verlinde, 2017).
EG is currently not as well developed as MOND, but its application to
spherically symmetric non-relativistic systems should be clear. It predicts
the same gravitational field around a point mass, but departs from MOND inside
the mass distribution of a galaxy due to an extra dependence on the density.
Observationally, EG works fairly well with the dwarf spheroidal satellites of
the MW (Diez-Tejedor et al., 2018), but it fails when confronted with
kinematic measurements from a larger sample of dwarfs (Bradford et al., 2015;
Pardo, 2020). For plausible assumptions about how EG should be applied to disc
galaxies, it predicts a distinct hook-shaped feature in the RAR, leading to
significant tension with observations (Lelli, McGaugh & Schombert, 2017a). In
particular, the required $M_{\star}/L$ becomes much lower than expectations
from stellar population synthesis modelling. Once EG is developed further so
that a non-relativistic field equation is available (e.g. Hossenfelder, 2017),
this should be applied to the SPARC sample to check whether the above-
mentioned criticisms are valid. EG also faces significant problems explaining
the behaviour of ultracold neutrons in the terrestrial gravitational field
(Kobakhidze, 2011a). While these results have been questioned (Chaichian et
al., 2011), significant experimental problems for EG were later reaffirmed
(Kobakhidze, 2011b). The basic assumptions underpinning the thermodynamic
approach to gravity advocated in EG also appear not to hold outside of
spherical symmetry (Wang & Braunstein, 2018).

Figure 12: Joint inference on the $I$-band $M_{\star}/L$ of DF44 and the
duration $\Delta t$ over which it formed its stars, based on its kinematics
(van Dokkum et al., 2019) assuming MOG and the typical $\Delta t$ values for
other galaxies (Pflamm-Altenburg & Kroupa, 2009), respectively. The
probability density is shown using different shades of green, with different
confidence regions also plotted as indicated. The red points show the results
of stellar population synthesis models in the framework of the integrated
galactic initial mass function (Kroupa & Weidner, 2003), with the error bars
arising from the assumed 20% uncertainty on the luminosity of DF44. This is
already factored into the uncertainty on the kinematically inferred
$M_{\star}/L_{I}$. The lowest tension is obtained for the shortest considered
$\Delta t$, in which case there is still $5.49\sigma$ tension between the
photometrically expected $M_{\star}/L_{I}$ and that required to reproduce the
kinematics of DF44 in MOG, also bearing in mind that such a short $\Delta t$
is atypical. Reproduced from figure 3 of Haghi et al. (2019b) using ApJL
author republication rights.

Unlike EG, a relativistic modified gravity theory is Moffat gravity (MOG;
Moffat, 2006). This is ruled out by the fact that GWs travel at a speed close
to $c$ (Boran et al., 2018). Focusing on non-relativistic tests, the main
problem with MOG is the lack of a fundamental acceleration scale, making it
extremely difficult to match the observed RAR. As a result, the MW RC is in
strong tension with MOG $-$ for any plausible baryonic distribution, the
predicted circular velocity at the Solar circle is $v_{c,\odot}<200$ km/s
(Negrelli et al., 2018). However, the observed value is $>200$ km/s at
overwhelming significance (e.g. McMillan, 2017; McGaugh, 2018). Even the
directly measured acceleration of the Solar System barycentre relative to
distant quasars proves that $v_{c,\odot}>200$ km/s at almost $5\sigma$
confidence (Bovy, 2020; Gaia Collaboration, 2021b). Moreover, MOG is
incompatible with DF44 at $5.49\sigma$ based on a joint fit to its star
formation history and observed $\sigma_{{}_{\text{LOS}}}$ profile (Haghi et
al., 2019b). This is evident in their figure 3, which we reproduce as our
Figure 12. It is possible to justify a higher $M_{\star}/L_{I}$ with a very
short star formation timescale because this would raise $M_{\star}/L_{I}$ in
the IGIMF context. But the kinematics are such that an exceptionally short
timescale would be required that is very uncommon for dwarfs (Pflamm-Altenburg
& Kroupa, 2009), and even then only a poor fit can be obtained to the observed
$\sigma_{{}_{\text{LOS}}}$ profile. Several other very serious problems for
MOG were reviewed in section 5.2 of de Martino et al. (2020). Therefore, MOG
is currently ruled out as a viable theory on galactic scales.

More generally, figure 10 of Famaey & McGaugh (2012) shows that the proportion
of gravity that is unaccounted for on galaxy scales does not correlate with
the galactocentric radius of the data point. We have reproduced the relevant
panel in our Figure 13, highlighting the difficulty of explaining the tight
observed correlation with $g_{{}_{N}}$ in a theory that supposes a fixed
length $r_{0}$ beyond which the departure from Newtonian dynamics sets in.
Since a flat RC implies that the gravity at larger radii must be
$GM/\left(rr_{0}\right)$, such a theory would predict that
$M_{s}\propto{v_{{}_{f}}}^{2}$. However, observations indicate that
$M_{s}\propto{v_{{}_{f}}}^{3.85\pm 0.09}$, with systematic uncertainties
perhaps allowing exponents in the range $3.5-4$ (Lelli et al., 2019). This is
consistent with the MOND prediction of 4, but by now it is quite clear that
$M_{s}\propto{v_{{}_{f}}}^{2}$ is strongly incompatible with observations, a
fact which has been known for many decades (McGaugh et al., 2000). This also
rules out volume corrections to the Hawking-Bekenstein entropy relation as a
cause for the missing gravity problem (Pérez-Cuéllar & Sabido, 2021).

Figure 13: The factor by which the radial gravity inferred from a galaxy’s RC
exceeds the $g_{{}_{N}}$ of its detected baryons, shown as a function of the
galactocentric radius of the data point. There is no clear correlation, unlike
when the dynamical discrepancy is plotted against $g_{{}_{N}}$ (Figure 3).
Reproduced from figure 10 of Famaey & McGaugh (2012).

Another possibility is called superfluid dark matter (SFDM; Berezhiani &
Khoury, 2015, 2016). The basic idea is that galaxies reside in DM halos which
dominate the mass, but these condense into a superfluid state in the central
$\approx 100$ kpc, with the baryons dominating the mass in the region probed
by RCs (out to $\approx 30$ kpc). The gravity between baryons is supplemented
by extra forces from phonons propagating in the superfluid, so essentially the
missing gravity problem is solved using additional non-gravitational forces.
The SFDM parameters have been estimated using galaxy clusters (Hodson et al.,
2017) where the superfluid nature of the DM is not very relevant, making the
paradigm similar to $\Lambda$CDM. As discussed in section 5.6 of Roshan et al.
(2021a), the LG satellite planes extend beyond plausible estimates for the
superfluid core radii of the MW and M31 (Berezhiani & Khoury, 2016;
Hossenfelder & Mistele, 2020). Since the planar geometry and coherently
rotating velocity trend make a primordial origin unlikely, the satellite
galaxies in these structures are most likely tidal in origin regardless of the
gravitational framework (Section 5.6). In SFDM, the more distant members of
these structures would experience just the Newtonian gravity of their baryons,
contradicting their high observed $\sigma_{{}_{\text{LOS}}}$ (McGaugh & Wolf,
2010; McGaugh & Milgrom, 2013a, b). Another issue is that the phonon force
acts only on the baryons, so it is difficult to understand why the RAR is also
evident in weak gravitational lensing of background galaxies (Section 3.5). It
is possible that the gravity of the SFDM (excluding phonon forces) conspires
to yield the RC-based RAR at the much larger distances probed by weak lensing
(Brimioulle et al., 2013; Milgrom, 2013; Brouwer et al., 2017, 2021), but this
would undermine the aim of SFDM to reproduce the RAR naturally. It has also
been argued that since stars on circular orbits in the outer disc would
generally move slower than the sound speed of the superfluid, they would emit
Cherenkov radiation, with the resulting energy loss causing the orbit to
shrink in a small fraction of a Hubble time (Mistele, 2021). Moreover, the
fact that MOND effects would persist only in the superfluid region implies
that a galaxy does not feel the EFE from another galaxy hundreds of kpc away.
This makes it difficult to explain the detection of the EFE by comparing RCs
of galaxies in isolated and more crowded environments (Chae et al., 2020b,
2021). Even if these issues could be overcome, since the phonon force must be
added to the Newtonian gravity, the deep-MOND behaviour can only be reproduced
if SFDM reduces to a MOND-like description in which the interpolating function
is $\nu=1+1/\sqrt{y}$, where $y\equiv g_{{}_{N}}/a_{{}_{0}}$. This gives
$\nu\left(1\right)=2$, but the simple $\nu$ function (Equation 19) and the MLS
form based on fitting the RAR (Equation 55) both give
$\nu\left(1\right)\approx 1.6$. Given the small observational uncertainties
(Figure 3), the transition between the Newtonian and low-acceleration regimes
appears to be too gradual in the SFDM picture. The inferred value of the
acceleration constant could be reduced to address this, but we would need a
reduction by a factor of $0.6^{2}$. There is limited scope for this given that
the data go down to very low $g_{{}_{N}}$ and thus probe the low-acceleration
regime, with the BTFR normalization essentially fixing $a_{{}_{0}}$ (e.g.
Lelli et al., 2019). This argument also applies to EG, and more generally to
any theory that attempts to reproduce the RAR by simply adding an extra force
to the Newtonian gravity $-$ the extra force would need to have the deep-MOND
behaviour (Equation 3), leading to a MOND-like description with interpolating
function $\nu=1+1/\sqrt{y}$. This does not agree with the observed RAR,
perhaps explaining why SFDM has difficulty matching galaxy RCs with plausible
$M_{\star}/L$ ratios $-$ especially if lensing constraints are also included
(Mistele, McGaugh & Hossenfelder, 2022). In addition, serious theoretical
difficulties have been identified with SFDM when trying to make it covariant
(Hertzberg et al., 2021).

The observationally required boost factor to $g_{{}_{N}}$ depends only on
$g_{{}_{N}}$ across a wide range of galaxies and is large if $g_{{}_{N}}\ll
a_{{}_{0}}$. This regime can be probed either by considering large radii in a
high-mass galaxy or the inner regions of a dwarf galaxy with a typically low
surface brightness. The two situations are rather different in many ways, but
$g_{{}_{N}}$ falls short of the actual $g$ by the same factor. Whatever causes
this boost to $g_{{}_{N}}$ seems to care only about the strength of the
gravitational field. We have seen that alternative explanations for this “one
law to rule them all” (Lelli et al., 2017b) beyond $\Lambda$CDM and MOND face
serious difficulties given the diversity of the so-called MOND phenomenology,
e.g. the fact that it also applies to weak lensing (Section 3.5). In addition,
gravitational dipole dark matter (Hajdukovic, 2020) and scale-invariant
dynamics (Maeder & Gueorguiev, 2020) have been strongly excluded on the basis
of Solar System ephemerides (Banik & Kroupa, 2020a, b). We therefore focus on
MOND as the main alternative to $\Lambda$CDM when it comes to addressing the
missing gravity problem.

##  4 Disc galaxy stability and secular evolution

When we observe the RC of a galaxy, there is a direct degeneracy between its
matter distribution and the the gravity law. This degeneracy could be broken
by the vertical forces, which are expected to be stronger in MOND (Section
3.1.2). This increases the relative importance of disc self-gravity, which
would affect the disc’s dynamical stability. Consequently, we may gain
important insights from the number of spiral arms and properties of the
central bar, especially the bar strength and pattern speed. These differences
are related to the gravity law, but the other side of this coin is that the
matter distribution of a galaxy also differs substantially between the
$\Lambda$CDM and MOND scenarios. We will see that the combination of these
factors leads to distinct predicted outcomes.

###  4.1 Survival of thin disc galaxies

The continued survival of thin disc galaxies provided an important motivation
for the hypothesis of CDM halos around galaxies (Hohl, 1971; Ostriker &
Peebles, 1973). These massive extended halos would cause galaxies to have a
large collision cross-section, making mergers quite common due to dynamical
friction between overlapping halos (Privon et al., 2013; Kroupa, 2015). Minor
mergers with the predicted dark satellites would also be common. These mergers
can be very disruptive for any embedded disc galaxies, so their observed
frequency and morphology might provide clues to how often mergers actually
occur. In particular, the $\Lambda$CDM scenario seems to imply that relatively
fragile thin disc galaxies are rather rare. Even if gas is accreted after a
major merger and subsequently reforms a thin disc, there should still be a
central bulge.

Therefore, the high prevalence of bulgeless thin disc galaxies appears
difficult to explain (Kautsch et al., 2006; Graham & Worley, 2008). Locally,
only 4/19 large galaxies ($v_{c}>150$ km/s) are ellipticals or contain a
significant classical bulge (Kormendy et al., 2010). Moreover, galaxy
interactions would eject stars into the halo, but the very low stellar halo
mass of M94 (Smercina et al., 2018) and M101 (van Dokkum et al., 2014) argues
against this scenario. Similar results have been obtained around NGC 1042 and
NGC 3351 (Merritt et al., 2016). Even these low-mass halos might not have been
accreted. If they formed in situ by ejection of stars from near the galaxy
centre, then one possible way to confirm this is a steep radial metallicity
gradient (Elias et al., 2018).

A related problem is that if thin discs form in $\Lambda$CDM, they would
typically have a significant fraction of material on orbits that deviate
substantially from circular motion, e.g. in a pressure-supported bulge
(Peebles, 2020). In its own words, “there is reasonably persuasive evidence
that simulations of galaxy formation based on the $\Lambda$CDM theory with
Gaussian initial conditions produce unacceptably large fractions of stars in
hot distributions of orbits.” Importantly, the fact that this conclusion is
based on many different simulations “suggests this is characteristic of
$\Lambda$CDM in a considerable range of ways to treat the complexity” of
subgrid baryonic physics processes.

Figure 14: The sky-projected aspect ratio distribution of galaxies with
$10^{10.0}<M_{\star}/M_{\odot}<10^{11.65}$ according to the Galaxy and Mass
Assembly (GAMA; Driver et al., 2009, 2011) and SDSS surveys (black and grey
lines, respectively), with GAMA having a higher resolution. The Poisson error
bars are shown for both surveys, but these are very small for SDSS. The
galaxies have been statistically weighted to match the $M_{\star}$
distribution of TNG50-1, the highest-resolution hydrodynamical cosmological
$\Lambda$CDM simulation suitable for a comparison (solid blue line). The green
line shows the extrapolation to a simulation with $8\times$ better mass
resolution in the CDM, with corresponding improvements elsewhere. The red line
shows the estimated result of improving the TNG50-1 CDM mass resolution by a
factor of $8^{5}$, beyond which results should change little as such a
simulation should be numerically converged. The shaded bands show the Poisson
uncertainties on simulated distributions. Reproduced from figure 8 of
Haslbauer et al. (2022) using ApJ author republication rights.

Lower mass halos would be more common in $\Lambda$CDM, especially due to the
implied presence of many purely DM halos with a negligible amount of baryons
(e.g. Kim et al., 2018). Consequently, we expect there to have been many more
minor mergers than major mergers. Even minor mergers could be quite
disruptive, so it seems very difficult to avoid substantial dynamical heating
of thin galactic discs. One possibility is that if the perturber is gas-rich,
it is more likely for the remnant to contain a dominant thin disc (Hopkins et
al., 2009).

To check if the proportion of thin disc galaxies is correct in $\Lambda$CDM,
it is necessary to compare the distribution of the sky-projected aspect ratio
$q_{\mathrm{sky}}$ in a large cosmological hydrodynamical simulation with that
in actual observations. This reveals a very significant $12.52\sigma$
discrepancy between the GAMA survey and TNG50-1, the highest resolution
$\Lambda$CDM simulation suitable for such a comparison (Haslbauer et al.,
2022). The “-1” suffix distinguishes the simulation from TNG50-2, TNG50-3, and
TNG50-4, which have successively lower resolution. This allows for a test of
numerical convergence. The aspect ratio distribution differs somewhat
depending on the resolution because thin galactic discs are still somewhat
challenging to resolve in large simulations (Ludlow et al., 2021). By using
the different TNG50 runs, it is possible to extrapolate to the result of a
higher-resolution simulation than TNG50-1, as done by Haslbauer et al. (2022)
in their figure 8 (reproduced in our Figure 14). If the CDM mass resolution is
improved by a factor of $8^{5}$ with corresponding improvements elsewhere
(mimicking the improvements from TNG50-2 to TNG50-1 repeated five times), the
tension with observations is reduced to $5.58\sigma$ using the more accurate
quadratic extrapolation discussed in their work. It is therefore clear that
improvements to the resolution of the simulations are likely to prove
insufficient: $\Lambda$CDM faces a genuine problem reproducing the high
observed fraction of thin disc galaxies. Haslbauer et al. (2022) also
highlighted severe deficiencies with earlier studies that claimed good
agreement in this regard (Vogelsberger et al., 2014; Bottrell et al., 2017a,
b). The problem might be alleviated in MOND due to a reduced frequency of
major mergers and because of the stronger vertical gravity from the disc
making it thinner for the same $\sigma_{z}$ (McGaugh & de Blok, 1998).

Since MOND reduces to Newtonian dynamics at high surface density and
acceleration, MOND cannot explain the continued survival of thin disc galaxies
with a surface density much above the MOND threshold
$\Sigma_{M}=137\,M_{\odot}/\text{pc}^{2}$ (Equation 52). Such a galaxy would
behave like a self-gravitating Newtonian disc, which is known to be unstable
(Hohl, 1971). Fortunately for MOND, disc galaxies do not have central surface
densities very far above $\Sigma_{M}$, which was observationally first noticed
as a very narrow concentration in the distribution of surface brightness
(Freeman, 1970; van der Kruit, 1987). This ‘Freeman limit’ is almost certainly
caused by selection effects (Disney, 1976) since it later became clear that
many galaxies exist with a much lower surface brightness (McGaugh et al.,
1995b). Even so, the upper limit appears to be real (McGaugh, 1996; Zavala et
al., 2003; Fathi, 2010) and is not easily assigned to a selection effect. We
illustrate this in Figure 15, which reproduces figure 5 of McGaugh (1996). The
sharp break evident at $\approx 22\,B$ magnitudes per square arcsecond
corresponds to $100\,L_{\odot}/\text{pc}^{2}$ (Willmer, 2018). Assuming
$M_{\star}/L_{B}\approx 1$, this implies a transition in the behaviour of
galaxies at a critical surface density of $\approx
100\,M_{\odot}/\text{pc}^{2}$, which is very similar to the predicted
stability threshold in MOND. Note that disc galaxies somewhat above this
threshold $\Sigma$ can still be stable in MOND, as demonstrated in Roshan et
al. (2021a) who considered a case where the central $\Sigma=10\,\Sigma_{M}$.
This is partly because the dynamical stability of a galaxy would depend on the
average surface density within roughly the central scale length. For a typical
exponential profile, this average surface density is smaller than the central
value by a factor of $2\left(1-2/\mathrm{e}\right)\approx 0.5$. In a MOND
context, we therefore do not expect a complete lack of high surface brightness
galaxies (HSBs), but we do expect that galaxies with a central
$\Sigma\gg\Sigma_{M}$ are much more likely to have become unstable at some
point, thereby contributing little towards the overall disc galaxy population.

Figure 15: The logarithmic number density of disc galaxies as a function of
central surface brightness in $B$ magnitudes per square arcsecond, with error
bars showing the Poisson noise. The solid lines show a broken power-law fit,
with a sharp break clearly evident at a luminosity close to
$100\,L_{\odot}/\text{pc}^{2}$, the significance of which is discussed in the
text. Reproduced from figure 5 of McGaugh (1996).

In principle, MOND should also explain the formation of disc galaxies out of
collapsing gas clouds. This was investigated for the first time by Wittenburg
et al. (2020), leading to an exponential surface density profile after 10 Gyr
of evolution. However, the work was not done in a cosmological context.
Hydrodynamical cosmological MOND simulations are currently underway (N.
Wittenburg et al., in preparation) based on the cosmological paradigm
discussed in Section 9.2. Gas accretion in a cosmological context is expected
to yield different results compared to starting with all the gas in a rotating
spherical cloud. Prior MOND simulations indicate that gas-rich clumpy galaxies
in the early universe do not secularly form bulges from the coalescence of
clumps (Combes, 2014b), helping to explain the high fraction of discs with
little to no bulge (Kormendy et al., 2010). Hydrodynamical MOND simulations of
M33 were also able to avoid the formation of a significant central bulge
(Banik et al., 2020). MOND appears more consistent with the high prevalence of
thin disc galaxies due to a reduced rate of mergers and the lack of bar-halo
resonances that can cause a strong bar at late times (Section 4.3).

###  4.2 Number of spiral arms

The greater importance of disc self-gravity in MOND manifests itself in the
different Toomre disc stability condition (Toomre, 1964) for AQUAL (Milgrom,
1989) and QUMOND (Banik et al., 2018b). In what follows, we will focus on
QUMOND, but the results are very similar in AQUAL. The local stability of a
stellar disc is determined by three main factors:

  1. 1.

The radial velocity dispersion $\sigma_{r}$,

  2. 2.

Disc self-gravity, and

  3. 3.

Shear caused by differential rotation of the disc.

Disc galaxy RCs are well understood in MOND (Section 3.1.1). Since the CDM
halo is chosen to match the RC, this is the same in both paradigms by
construction. Therefore, differences in disc self-gravity feed through to
differences in the minimum $\sigma_{r}$ required for local stability.

The Toomre condition relies on the notion that there is a particular
wavelength $\lambda_{c}$ which is least stable to perturbations. This is
because very short wavelength modes are essentially just sound waves little
affected by gravity due to the short crossing time, while long wavelength
modes rapidly wind up in a differentially rotating disc. Intermediate length
modes are significantly affected by self-gravity, yet are not so long as to be
easily wound up. This is explained further in Binney & Tremaine (2008), whose
section 6.2.3 shows how these considerations lead to a parabolic dispersion
relation between the inverse length scale of perturbations and the square of
their oscillation frequency (negative values indicate the existence of an
unstable growing mode). Though their result is valid for an isothermal fluid,
the result for a collisionless stellar disc is very similar. According to
equation 22 of Toomre (1964), a Newtonian disc with the lowest $\sigma_{r}$
necessary for local stability has

$\displaystyle\lambda_{c}~{}=~{}0.55\times\frac{4\mathrm{\pi}^{2}G\Sigma}{{\Omega_{r}}^{2}}\,,$
(63)

where $\Sigma$ is the disc surface density and $\Omega_{r}$ is the radial
epicyclic frequency determined by the RC. The factor of 0.55 was derived
numerically in section 5c of Toomre (1964). According to its equation 65, this
mode is stable in the local approximation if

$\displaystyle\sigma_{r}~{}\geq~{}\frac{3.36\,G\Sigma}{\Omega_{r}}\,.$ (64)

The ratio of $\sigma_{r}$ to this critical threshold is known as the Toomre
parameter $Q$, which in Newtonian gravity is thus

$\displaystyle Q_{N}~{}\equiv~{}\frac{\sigma_{r}\Omega_{r}}{3.36\,G\Sigma}\,.$
(65)

Discs are locally stable if $Q>1$ everywhere, including near the centre, where
$Q$ is typically lower than in the outskirts. Note that Equation 63 assumes
$Q=1$, which we expect due to self-regulatory mechanisms (discs with $Q\ll 1$
should become dynamically hotter). In general, $\lambda_{c}\propto Q^{2}$.

As derived in Banik et al. (2018b), the above-mentioned classical Toomre
condition for a thin disc can be generalized to QUMOND by redefining $G\to
G\nu\left(1+K/2\right)$, where the MOND interpolating function $\nu$ (e.g.
Equation 19) and its logarithmic derivative $K$ (Equation 42) must be
evaluated just outside the disc plane so as to include the vertical component
of the Newtonian gravity. With this adjustment, the QUMOND version of Equation
65 is

$\displaystyle
Q_{M}~{}\equiv~{}\frac{\sigma_{r}\Omega_{r}}{3.36\,G\Sigma\nu\left(1+\frac{K}{2}\right)}\,.$
(66)

The extra factor of $\nu$ can be very large in low surface brightness (LSB)
discs, where $K$ is close to the deep-MOND value of $-1/2$. In the Newtonian
regime, $\nu\to 1$ and $K\to 0$, recovering the standard behaviour.

By considering QUMOND models of M33 with different values of the Toomre
parameter, Banik et al. (2020) demonstrated that using Equation 66 as the
QUMOND Toomre condition provides a good guide to where the transition lies
between stable and unstable discs (see their figure 24). This is reproduced
here as our Figure 16, which shows the effect of lowering the minimum allowed
$Q$ in the initial conditions (this parameter is called $Q_{\text{lim}}$ in
the code). When $Q_{\text{lim}}$ is reduced from 1.25 (black curve) to 1.1
(blue curve), the disc is initially dynamically colder and also ends up with a
lower $\sigma_{z}$. But if $Q_{\text{lim}}$ is reduced further from 1.1 to 1
(red curve), then $\sigma_{z}$ starts rising rapidly after $\approx 5$ Gyr,
ending up higher than in the other models despite being dynamically colder
initially. This is suggestive of a dynamical instability which is avoided by
the higher $Q_{\text{lim}}$ models. To ensure stability and leave some safety
margin, we recommend using $Q_{\text{lim}}=1.25$ in MOND simulations where the
goal is to have a stable disc.

Figure 16: Time evolution of the vertical velocity dispersion $\sigma_{z}$ of
stars at radii of $0.5-3$ kpc in a hydrodynamical MOND simulation of M33 that
includes a weak EFE representative of the gravity from M31. Different lines
show different values for $Q_{\text{lim}}$, the minimum Toomre parameter
(Equation 66) imposed on the initial conditions. Notice how reducing
$Q_{\text{lim}}$ from 1.25 (solid black) to 1.1 (dotted blue) reduces the
final $\sigma_{z}$, as expected for a dynamically colder initial disc.
However, $\sigma_{z}$ ends up higher when $Q_{\text{lim}}=1$ (solid red),
suggesting that the MOND generalization of the Toomre condition does indeed
provide a good guide to which models are stable and which are not. Reproduced
from figure 24 of Banik et al. (2020) using ApJ author republication rights.

The higher effective $G$ in MOND makes discs less stable for a given velocity
dispersion, so they must be dynamically hotter to remain locally stable
(McGaugh & de Blok, 1998). In practice, this means they would heat up to a
greater extent through dynamical instabilities until they become Toomre
stable. Therefore, the observation of an extremely low velocity dispersion
could challenge the MOND scenario. But if observations confirm the higher
velocity dispersion it requires, this could be interpreted in the standard
context as a disc which is dynamically hotter than required for local
stability. Since the Toomre condition only provides a lower limit on
$\sigma_{r}$, this is plausible, though other considerations might break the
degeneracy.

Interestingly, observed LSB galaxies (LSBs) have rather higher velocity
dispersions than required in the CDM picture $-$ they appear to be dynamically
overheated (Saburova, 2011). If so, it would be difficult for LSBs to have
ongoing star formation and sustain spiral density waves, the leading
explanation for observed spiral features in HSBs like the MW and M31 (Lin &
Shu, 1964). However, LSBs also have spiral features (Schombert et al., 1992;
de Blok et al., 1995; McGaugh et al., 1995a; Impey et al., 1996). This and
other features of LSBs suggest that their gravitating mass mostly resides in
their disc, contradicting $\Lambda$CDM expectations and instead confirming the
MOND prediction (section 3.3 of McGaugh & de Blok, 1998).

The extra MOND factors in Equation 66 also affect $\lambda_{c}$, which in
QUMOND becomes

$\displaystyle\lambda_{c}~{}=~{}\frac{2.2\mathrm{\pi}^{2}G\Sigma\nu\left(1+\frac{K}{2}\right)}{{\Omega_{r}}^{2}}\,.$
(67)

The least stable mode is thus longer in MOND (cf. Equation 63). Since the
amplitude of density fluctuations is presumably larger for a less stable mode,
we expect that the most prominent structures in a galaxy would differ
depending on the gravity law. In a disc galaxy, the most striking feature is
usually its pattern of spiral arms. Assuming the Lin-Shu density wave theory
for spiral structures (Lin & Shu, 1964), the number of spiral arms allows for
an estimate of $\lambda_{c}$ (Sellwood & Carlberg, 1984; Athanassoula et al.,
1987). Using this principle, D’Onghia (2015) analytically predicted the number
of spiral arms in galaxies observed as part of the DiskMass survey. She found
good agreement with observations for a reasonable $M_{\star}/L$. The DiskMass
survey “selects against LSB discs”, making the work of D’Onghia (2015) an
important check on the validity of the mode-counting technique in a regime
where $\Lambda$CDM and MOND give similar predictions.

The density wave theory should also apply in LSBs. Applying a similar
technique to LSBs under the assumption of Newtonian gravity, Fuchs (2003) and
Saburova & Zasov (2013) found that LSB discs would be much too stable to
permit the formation of their observed spiral arms for plausible $M_{\star}/L$
ratios consistent with stellar population synthesis models (e.g. Bell & de
Jong, 2001). The spiral structures seem to imply that much of the mass needed
to explain the elevated RCs with standard gravity resides not in a stabilizing
CDM halo but rather within the disc itself. This would make the discs very
massive, sometimes requiring ${M_{\star}/L_{R}>10\times}$ the Solar value.
Physically, this arises because density perturbations would wind up quickly in
a system where a dominant dark halo elevates the RC but not the disc self-
gravity. This scenario would make the disc very stable and thus unable to
easily form a grand-design spiral, which however is very common
observationally (Hart et al., 2017). While some cases like the well-known M51
are tidally induced (Dobbs et al., 2010), it seems unlikely that all grand-
design spirals are of tidal origin. Their high prevalence is more easily
understood as a sign that the least stable wavelength (Equation 63 or 67) is
rather long, favouring an enhancement to the effective value of $G$.

###  4.3 Bar strength

In the previous section, we assumed that CDM halos (if present) would act to
damp perturbations in the disc. This is reasonable because a physical halo is
distinct from the disc, unlike a phantom halo which would respond
instantaneously to changes in the disc (the response would occur within a few
light crossing times). However, a physical halo would also respond, albeit
over a longer timescale. An important characteristic of a real halo is that it
can absorb angular momentum from the disc, which is not possible for a phantom
halo. Moreover, a physical halo opens up the possibility of resonances with
the disc, i.e. that perturbations to the disc affect the halo which then
further affects the disc. We therefore discuss if the secular evolution of
disc galaxies is more easily understood with a real surrounding halo, or if
postulating such a halo would cause effects that are not observed. We begin by
describing $\Lambda$CDM expectations, then discuss the observations, and
finally consider the MOND scenario, focusing in particular on how the observed
morphology of M33 might be understood in both paradigms.

Early $N$-body simulations of self-gravitating Newtonian discs quickly became
dominated by bars (Miller & Prendergast, 1968; Hockney & Hohl, 1969; Hohl,
1971). These works suggested an almost unavoidable global instability of a
self-gravitating rotating stellar disc. It can be made more stable if it is
embedded within a dominant CDM halo (Ostriker & Peebles, 1973). However, their
equation 6 highlights an important deficiency in the analysis $-$ the halo was
assumed to provide a time-independent contribution to the potential. It was
not treated using live particles. This is sometimes known as a rigid halo
model, which can be thought of as including extra halo particles that provide
gravity but do not move.

Subsequent simulations with live halos have shown that over long periods, a
disc inside a massive halo is bar unstable due to a resonance between the
orbits of stars in the bar and DM particles in the inner halo (Athanassoula,
2002). The reduced role of disc self-gravity means that a massive halo can
delay the bar instability, but it is even stronger when it does occur, which
would still be well within the lifetime of the universe. This runs contrary to
the role of a halo in reducing the relative importance of disc self-gravity,
which might naively be expected to stabilize the disc (Section 4.2). While
this is true on short timescales, discs within a more centrally concentrated
halo end up with a stronger bar that is longer and thinner (Athanassoula &
Misiriotis, 2002).

For a bar to develop, it is necessary for stellar orbits to become
significantly non-circular. While the orbits need not be as elliptical as the
bar itself, some ellipticity is needed because otherwise any non-axisymmetric
feature would quickly wind up and disappear due to differential rotation of
the disc (Lin & Shu, 1964). Since elliptical orbits have less angular momentum
than a circular orbit for the same apocentre, it is clear that redistribution
of angular momentum is important to the bar instability. A critical issue is
how much angular momentum can be transferred by the bar, which in turn depends
on where it is being transferred to. If the sink is the outer disc where there
is relatively little material, the bar is unlikely to remain very strong over
an extended period $-$ the outer disc would expand further, making it less
tightly coupled to the bar region. But if the sink is a dominant halo, then a
very large amount of angular momentum can be transferred out of the inner
disc, leading to a rather strong bar (Athanassoula, 2003). Therefore, a halo
serves to increase the maximum possible amplitude of the bar to very high
values, even if it slows down the bar instability.

The halo would typically have a small amount of net rotation in the same sense
as the disc, because both would have been spun up by tidal torques from large-
scale structure (Doroshkevich, 1970; White, 1984). A rotating halo further
worsens the bar instability because it allows more DM to be on a resonant
orbit with the rotating bar (Saha & Naab, 2013). Those authors found that a
counter-rotating halo could suppress the bar instability, but this would be a
rare scenario as tidal torques from large-scale structure would typically spin
up the baryonic and DM components in the same sense. This is ultimately just a
consequence of the weak equivalence principle or the universality of free fall
(which remains valid even in MOND).

The $\Lambda$CDM expectation of typically rather strong bars was confirmed by
more recent $N$-body simulations (Sellwood, 2016; Berrier & Sellwood, 2016).
Scattering by giant molecular clouds delays the onset of the bar instability,
but is unable to prevent it (Aumer et al., 2016). This suggests that including
the gas component would have little impact on the results.

The larger scale cosmological environment could be important, though more
computationally demanding to include. This was attempted by Bauer & Widrow
(2019) using a technique to insert a disc into a DM-based cosmological
simulation (Bauer et al., 2018). Section 5.2 of Bauer & Widrow (2019)
describes the setup of their two cosmological disc simulations. Their figure
12 shows a rapid increase in the bar strength once the disc is switched from
rigid to live. Bars in cosmologically motivated disc environments are
therefore expected to develop in only a few Gyr (Bauer & Widrow, 2019),
somewhat contrary to what is stated in their abstract.

Since the bar is a self-gravitating instability, a dynamically hotter and
thicker disc could suppress the development of a bar. However, any bar in a
thicker disc would itself be thicker, making the bar less prone to buckle.
This leads to an even longer bar (Klypin et al., 2009). Indeed, their figure
12 shows that the bar in their thick disc model ends up covering the entire
disc. It therefore seems clear that if disc galaxies have a substantial amount
of CDM in a spheroidal halo, then gravitational interactions with the halo
would lead to a quite strong bar on a timescale much shorter than the age of
the universe. Differences in the timescale of this instability are less
important than what limits the magnitude of the bar, simply because the
universe is many dynamical timescales old.

Turning now to the real Universe, Masters et al. (2011) found that $\approx
30\%$ of disc galaxies have bars based on optical SDSS images (SDSS
Collaboration, 2000). The bar fraction doubles in the higher angular
resolution Spitzer Survey of Stellar Structure in Galaxies (S4G; Erwin, 2018).
Their estimated bar fraction of $\approx 60\%$ is in agreement with earlier
estimates using near-infrared data (Eskridge et al., 2000; Sheth et al.,
2008). Moreover, Eskridge et al. (2000) found that “strong bars are nearly
twice as prevalent in the near-infrared as in the optical.” Similar
conclusions were reached by Erwin (2018), who argued that the superior
resolution of the Spitzer Space Telescope allowed the detection of shorter
bars that might be missed in SDSS images. This may well explain the higher bar
fraction in higher resolution images, though the reduced dust extinction at
longer wavelengths might help as well. It is therefore clear that the majority
of disc galaxies have a bar, but $\approx 40\%$ are unbarred or have only a
very weak bar.

One such weakly barred galaxy is the rather isolated M33 (Corbelli &
Walterbos, 2007; Kam et al., 2015). Its RC requires a substantial amount of DM
in conventional gravity (Sanders, 1996; Famaey & McGaugh, 2012; Kam et al.,
2017), as expected in MOND for its surface brightness. However, unavoidable
gravitational interactions between the disc and the required dark halo make it
very difficult to explain the observed properties of M33 (Sellwood, Shen & Li,
2019). Those authors conducted a detailed investigation into the global
stability of M33 in the $\Lambda$CDM context, finding that a strong bar
develops in only a few Gyr. Their section 3.2 discusses how they obtained
similar results when they included gas with thermal or mechanical feedback
from the stars. The authors were able to suppress the bar instability by using
a rigid halo which provides a fixed potential, corresponding to a halo whose
particles do not move. Though unphysical, this demonstrated that the
instability is indeed due to the halo (see their section 4.1). A more
plausible solution is to reduce the disc mass, but this entails an
unrealistically low $M_{\star}/L$ of only 0.6 (0.23) in the $V$ ($K$) band
(Conroy, 2013; Schombert et al., 2019). Moreover, a lower mass disc would be
less maximal and require a larger halo contribution, reducing the least stable
wavelength (Equation 63). This caused Sellwood et al. (2019) to conclude in
their section 4.2.1 that “all the spiral patterns in our simulation of this
low-mass disc were multi-armed”, making it difficult to understand the
observed bisymmetric spiral in M33 (Smith et al., 1984; Elmegreen et al.,
1992; Regan & Vogel, 1994; Jarrett et al., 2003; Corbelli & Walterbos, 2007).
It could be made bar stable with a higher velocity dispersion and thicker disc
in which $Q=2$, showing that combining such a hot disc with a dominant DM halo
might perhaps give M33 the sought-after stability in a $\Lambda$CDM context.
However, it is difficult to understand how such a dynamically overheated disc
is currently forming stars (Verley et al., 2009). In addition, the simulation
does not yield a grand-design spiral because the least stable wavelength is
too short. The difficulties reproducing the observed properties of M33 led
Sellwood et al. (2019) to conclude that “none of the ideas proposed over the
subsequent years can account for the absence of a bar in M33” in the
$\Lambda$CDM context. Generalizing the implications to other weakly barred
galaxies, the authors went on to conclude that “it is shocking that we still
do not understand how the bar instability is avoided in real galaxies.”

In a MOND context, the enhanced role of disc self-gravity would make the most
unstable wavelength longer (Equation 67), helping to address the dominant two-
armed spiral. The lack of bar-halo angular momentum exchange would remove the
main mechanism causing difficulty for the simulations of Sellwood et al.
(2019). However, the enhanced disc self-gravity means that a MOND model of M33
would be much less stable at early times. MOND is therefore unable to suppress
the bar instability, potentially leading to a conflict with observations.
Consistency would require the bar strength to saturate at a rather low value
consistent with observations, or for the bar to become very strong and then
get destroyed. This is theoretically possible because if angular momentum can
only be redistributed within the disc, then the outer disc might not have
enough mass to sustain a very strong bar over a Hubble time. Clearly,
numerical simulations are required to address such subtle issues.

The first $N$-body MOND simulation of an isolated disc was the seminal work of
Brada & Milgrom (1999), who used a custom potential solver to advance
particles in 2D. This helped to clarify that MOND discs would be more stable
than a DM-free Newtonian disc, but less stable than a Newtonian disc in a DM
halo adjusted to reproduce the MOND RC. Perhaps the main outcome of the study
was that the stability problem of Newtonian self-gravitating discs (e.g. Hohl,
1971) might be solved in MOND, consistent with earlier analytic expectations
(Milgrom, 1989). Clearly, a dominant DM halo around galaxies is not the only
way to stabilize an embedded disc.

The 3D dynamics of disc galaxies were explored in a series of simulations by
Tiret & Combes (2007). These revealed the expected rapid strengthening of the
bar due to enhanced disc self-gravity, but the bar then weakened as
$\sigma_{z}$ increased (see their section 5.1). The same basic picture was
obtained when those authors included gas as sticky particles in their
simulations, though including dissipation in this way weakened the bar
somewhat (Tiret & Combes, 2008).

The stability of M33 was addressed in a recent hydrodynamical MOND simulation
(Banik et al., 2020). MOND is expected to have a significant effect on M33
because its surface density is below the critical MOND threshold (Equation 52)
outside the very central region (see their figure 1). The gravity law
governing M33 is thus rather different in MOND, with corresponding differences
in its total mass (i.e. the lack of a dominant dark halo). Despite very little
fine-tuning, figure 12 of Banik et al. (2020) demonstrates that the bar
strength has the sought-after behaviour of a sharp rise at early times
followed by a decline to a very low strength well within a Hubble time. We
show this in our Figure 17, which reproduces their figure 19 in which the EFE
from M31 is also approximately included. The bar strength evolved similarly in
earlier simulations at lower resolution, which also show a sharp initial rise
in the bar strength followed by a more gradual decline (Tiret & Combes, 2007,
2008).

Figure 17: Time evolution of the $\mathrm{m}=2$ azimuthal Fourier mode in a
hydrodynamical MOND simulation of M33 (Banik et al., 2020), found from a
Fourier decomposition of material at radii of $0.5-3$ kpc. Results are shown
for the stars (solid red), gas (dotted blue), and the total (solid black). For
comparison, the figure also shows a purely stellar $N$-body simulation (dotted
red) and the most realistic CDM simulation of Sellwood et al. (2019) in
Newtonian gravity (dot-dashed magenta), which they called their “full mass”
model. Notice that the CDM model develops a stronger bar than the
observationally estimated strength of 0.2 (Corbelli & Walterbos, 2007), shown
here as a grey line. The bar is expected to strengthen further with time. In
MOND, after a very strong bar instability at early times, the bar ends up
being rather weak. Reproduced from figure 19 of Banik et al. (2020) using ApJ
author republication rights.

In addition to a realistic bar strength at late times, figure 9 of Banik et
al. (2020) demonstrates that their model develops a bisymmetric spiral, which
is probably related to the most unstable wavelength being longer in MOND
(Equation 63). Interestingly, Banik et al. (2020) also considered the EFE from
M31, finding that it improves the agreement with some observables like the RC
in the central 2 kpc (see their figure 13). The more realistic models did not
develop a central bulge, with the EFE further suppressing the tendency for
material to concentrate towards the centre (see their figure 14). Therefore,
the global properties of M33 are reasonably well reproduced in MOND despite
the smaller theoretical flexibility it offers (videos of their isolated
simulations are available here:
https://seafile.unistra.fr/d/843b0b8ba5a648c2bd05/, accessed on 28/05/2022).
While this by no means proves the MOND model correct, various aspects of M33
like its weak bar argue against the presence of a dominant DM halo (Sellwood
et al., 2019).

###  4.4 Bar fraction

The properties of bars can be considered more generally, at least in the
$\Lambda$CDM paradigm where the galaxy population can be investigated in the
proper cosmological context. A major error with the evolution of bars should
be evident in the fraction of galaxies which have a bar, defined in Roshan et
al. (2021a) as those galaxies where the strength of the $\mathrm{m}=2$
azimuthal Fourier mode ($A_{2}$) in units of the azimuthal average ($A_{0}$)
is $A_{2}/A_{0}>0.2$ at the radius where this ratio peaks. The study made use
of the cosmological $\Lambda$CDM simulation known as Illustris The Next
Generation (hereafter TNG; Nelson et al., 2019), an improvement on the earlier
Illustris simulation (Vogelsberger et al., 2014). For the comparison, Roshan
et al. (2021a) considered TNG100 and TNG50, the highest resolution simulation
in the TNG suite. They also considered an earlier simulation called Evolution
and Assembly of Galaxies and their Environments (EAGLE; Crain et al., 2015;
Schaye et al., 2015), as well as observational results based on SDSS (Oh et
al., 2012) and a much more recent study called S4G that uses data from the
Spitzer Space Telescope (Erwin, 2018). As argued in Section 4.3, S4G should
provide a much more accurate estimate, especially at the low-mass end where
shorter bars can be difficult to resolve in SDSS. The bar fractions in
different numerical and observational studies are shown in Figure 18, which
reproduces figure 1 of Roshan et al. (2021b). It is clear that there is a very
significant disagreement between the latest observations and the bar fraction
expected in $\Lambda$CDM, which is similar between TNG100 and the much higher
resolution TNG50 (albeit slightly higher than in EAGLE100). The numbers refer
to the approximate side length of the cubic simulation volume, expressed in
co-moving Mpc. A similar analysis for TNG100 had previously been conducted by
Zhao et al. (2020) and gave similar results (the grey curve). The zoom-in
NewHorizon simulation (Dubois et al., 2021) also produces too few barred
galaxies at the low-mass end (Reddish et al., 2022) despite having a multi-
phase interstellar medium.

Figure 18: The fraction of disc galaxies with a bar ($A_{2}/A_{0}>0.2$) in the
surveys known as SDSS (cyan; Oh et al., 2012) and S4G (blue; Erwin, 2018),
with the latter having a much higher angular resolution due to the use of data
from beyond the atmosphere. These observational results are shown with dashed
lines, while solid lines are used for simulation results. The bar fractions
are similar in the TNG50 simulation (black) and TNG100 (red), with a similar
determination for TNG100 by a different team shown in grey (Zhao et al.,
2020). The TNG results appear to have converged, though the bar fraction is
somewhat higher than in the older EAGLE100 simulation (brown). Reproduced from
figure 1 of Roshan et al. (2021b).

The bar fraction decreases at high $M_{\star}$, but in $\Lambda$CDM it is
expected to rise. This could be due to lower mass galaxies typically having a
larger proportion of CDM to baryons, which overstabilizes the disc and
prevents bar formation. However, the results of Sellwood et al. (2019)
indicate that the DM halo can promote bar growth. Since the disc needs to be
quite cold dynamically to form a bar, the low bar fraction at low $M_{\star}$
might instead be related to such discs generally being thicker than in
reality, i.e. having a greater degree of pressure support. This may be related
to the underpredicted fraction of galaxies which are thin discs, possibly due
to baryonic feedback effects and heating by DM subhalos (Section 4.1). Lacking
subhalo heating, the isolated model of Sellwood et al. (2019) was able to
retain a rather thin disc in their M33 simulation $-$ and thus a strong bar.

At the high-mass end, the high bar fraction in $\Lambda$CDM simulations
(Figure 18) might be due to the disc being dynamically colder as subhalo
heating and feedback would be less disruptive for a more massive galaxy. A
major merger would thicken the disc, but in this case the galaxy might not be
considered a disc any more $-$ and thus removed from the sample used to
estimate the fraction of disc galaxies with a bar. Massive simulated disc
galaxies seem to have much more DM than they need to maintain consistency with
the observed RAR (Marasco et al., 2020), which is the disc galaxy analogue to
the similar problem for ellipticals (Figure 7). This makes it possible for the
disc to be Toomre stable with a rather low velocity dispersion, and thus able
to sustain a strong bar. Resonant bar-halo interactions can then strengthen
the bar further (Section 4.3).

The interpretation of Figure 18 is also far from clear in MOND, which did not
predict the observed trend. In MOND, we expect that bar properties would
differ between galaxies depending mainly on whether their central surface
density $\Sigma_{0}$ is below or above the critical threshold $\Sigma_{M}$
(Equation 52). Since higher mass galaxies generally have a higher
$\Sigma_{0}$, we interpret Figure 18 as showing a lower bar fraction at high
$\Sigma_{0}$ rather than at high $M_{\star}$. It would be useful to plot the
observed bar fraction in terms of central surface brightness rather than
$M_{\star}$, but we will assume in what follows that the trend would remain
similar. In the study of Banik et al. (2020) on M33,
$\Sigma_{0}\approx\Sigma_{M}$ (see their figure 1). The exponential disc
galaxy model considered in Roshan et al. (2021a) had
$\Sigma_{0}=10\,\Sigma_{M}$, which as discussed in their section 3.3 is
intermediate between the values for the MW and M31, thereby representing a
much higher mass galaxy than M33. In the lower $\Sigma_{0}$ model
representative of M33, the bar strength rises rapidly to $\approx 1$ (Figure
17), so there should be a high likelihood of observing a strong bar if we
consider the time evolution of a single model as giving similar statistics to
a galaxy population observed at the same epoch. But in the higher $\Sigma_{0}$
simulation representative of the major LG galaxies, the peak bar strength is
only $\approx 0.3$ (see figure 13 of Roshan et al., 2021a). The bar weakens
substantially by the end of both simulations to a similar final strength, but
there is a significant difference in the peak strength. For a comparison with
the purely stellar simulations in Roshan et al. (2021a), it is best to
consider the stellar-only model of M33 shown in figure 10 of Banik et al.
(2020), which gives an even stronger bar that is still well above the
typically used 0.2 threshold by the end of the simulation. As a result, we
might expect that we are much more likely to observe a bar in a lower
$\Sigma_{0}$ galaxy like M33. This is in line with the observed declining
trend in the bar fraction as a function of $M_{\star}$ (Figure 18).

As with the $\Lambda$CDM case, we assign the different behaviour to
differences in the degree of random motion. For a fair comparison, we first
need to consider the typical circular velocity of each galaxy, for which we
use Equation 53. The M33 model of Banik et al. (2020) has $M=6.5\times
10^{9}\,M_{\odot}$ ($v_{{}_{f}}=101$ km/s) while the MOND model of Roshan et
al. (2021a) has $M=8.57\times 10^{9}\,M_{\odot}$ ($v_{{}_{f}}=108$ km/s), so
both models have a similar RC amplitude. The random motions in the initial
conditions used by Roshan et al. (2021a) are given in their figure 5. We will
focus on $\sigma_{z}$, which is shown in the bottom right panel. The central
$\sigma_{z}\approx 30$ km/s, much higher than the $\approx 20$ km/s evident in
Figure 16 for the M33 model. This figure also shows that $\sigma_{z}$ rises by
a factor of $\approx 1.6$ over the full duration, with the factor being
$\approx 2$ in the stellar-only model shown in figure 6 of Banik et al.
(2020). The disc heating in the model of Roshan et al. (2021a) is shown in
their figure 26 in terms of the disc thickness $h$. We expect that
${\sigma_{z}}^{2}$ is proportional to the product of $h$ and the vertical
gravity. Assuming a thin disc, the vertical gravity remains similar as we move
up from the disc plane, so $\sigma_{z}\propto\sqrt{h}$. In reality, we expect
a slightly steeper scaling because at large radii, the vertical gravity is not
related to the local surface density but instead arises from the change in
angle towards the central regions of the galaxy, i.e. we can estimate the
result by considering the galaxy to be a point mass. In this regime, the
vertical gravity $\propto h$, so we expect that $\sigma_{z}\propto h$.
However, since the disc thickness shown in figure 26 of Roshan et al. (2021a)
is not measured that far out, it should be more accurate to use
$\sigma_{z}\propto\sqrt{h}$, which will give a conservative estimate of the
rise in $\sigma_{z}$ during the simulation. Since this figure shows $h$ rising
$\approx 4\times$, we can safely say that $\sigma_{z}$ doubles. This is a
similar rate of disc heating to the M33 model of Banik et al. (2020), if not