pdftexdestination with the same hyperrefOption hyperrefToken pdftex(dest)

# Stability of disks in quasilinear MOND

Indranil Banik1,2, Mordehai Milgrom3 and Hongsheng Zhao1
1Scottish Universities Physics Alliance, University of Saint Andrews, North
Haugh, Saint Andrews, Fife, KY16 9SS, UK
2Helmholtz-Institut für Strahlen und Kernphysik (HISKP), University of Bonn,
Nussallee 14$-$16, D-53115 Bonn, Germany
3Department of Particle Physics and Astrophysics, Weizmann Institute, Rehovot
7610001, Israel Email<EMAIL_ADDRESS>(Indranil Banik)
<EMAIL_ADDRESS>(Mordehai Milgrom)
<EMAIL_ADDRESS>(Hongsheng Zhao)

###### Abstract

We consider disk stability in the quasilinear formulation of MOND (QUMOND),
the basis for some $N$-body integrators. We derive the generalisation of the
Toomre criterion for the stability of disks to tightly wound, axisymmetric
perturbations. We apply this to a family of thin exponential disks with
different central surface densities. By numerically calculating their QUMOND
rotation curves, we obtain the minimum radial velocity dispersion required for
stability against self-gravitating collapse. MOND correctly predicts much
higher rotation speeds in low surface brightness galaxies (LSBs) than does
Newtonian dynamics without dark matter. Newtonian models thus require putative
very massive halos, whose inert nature implies they would strongly stabilize
the disk.

MOND also increases the stability of galactic disks, but in contradistinction
to Newtonian gravity, this extra stability is limited to a factor of 2. MOND
is thus rather more conducive to the formation of bars and spiral arms.
Therefore, observation of such features in LSBs could be problematic for
Newtonian galaxy models. This could constitute a crucial discriminating test.
We quantitatively account for these facts in QUMOND.

We also compare numerical QUMOND rotation curves of thin exponential disks to
those predicted by two algebraic expressions commonly used to calculate MOND
rotation curves. For the choice that best approximates QUMOND, we find the
circular velocities agree to within 1.5% beyond $\approx 0.5$ disk scale
lengths, regardless of the central surface density. The other expression can
underestimate the rotational speed by up to 12.5% at one scale length, though
rather less so at larger radii.

###### keywords:

gravitation – instabilities – dark matter – galaxies: kinematics and dynamics
– methods: analytical – methods: numerical

††pubyear: 2018††pagerange: Stability of disks in quasilinear MOND–A

##  1 Introduction

MOND (Milgrom, 1983a) is an alternative to the dark matter (DM) hypothesis in
accounting for the observed dynamical discrepancies in galactic systems,
especially between their measured rotation curves (RCs) and those predicted by
Newtonian gravity (e.g. Rubin & Ford, 1970; Rogstad & Shostak, 1972; Roberts &
Whitehurst, 1975). A discrepancy is also apparent in the ‘timing argument’: in
Newtonian gravity, the visible masses of the Galaxy and M31 are insufficient
to turn their initial post-Big Bang recession around to the observed extent
(Kahn & Woltjer, 1959). MOND posits that these acceleration
discrepancies111sometimes called mass discrepancies, though the reason for the
discrepancy may not be missing mass are not due to the presence of DM but
arise from a breakdown of Newtonian dynamics that is reckoned without in
analysing the dynamics of these systems. MOND is extensively reviewed in
Famaey & McGaugh (2012) and Milgrom (2014). The M31 timing argument in MOND
was discussed by Zhao et al. (2013) and detailed calculations were presented
in section 2 of Banik et al. (2018).

MOND introduces $a_{{}_{0}}$ as a fundamental acceleration scale of nature.
When the gravitational field strength $g\gg a_{{}_{0}}$, standard dynamics is
restored. However, when $g\ll a_{{}_{0}}$, the dynamical equations become
scale-invariant (Milgrom, 2009b). In this deep-MOND regime, an inevitable
consequence of scale invariance is that, asymptotically far outside a
distribution of total mass $M$, the rotational speed of a test particle
becomes independent of its distance $R$ from the mass. This occurs when $R\gg
r_{{}_{M}}\equiv\sqrt{GM/a_{{}_{0}}}$, where $r_{{}_{M}}$ is the MOND radius
of the mass $M$. Therefore, in a modified gravity formulation of MOND,
$g\propto R^{-1}$ beyond the MOND radius. Applying dimensional arguments to
the fact that $a_{{}_{0}}$ is the only additional constant in MOND, we must
have that

$\displaystyle g~{}\propto~{}\frac{\sqrt{GMa_{{}_{0}}}}{R}~{}~{}~{}\text{for
}~{}R~{}\gg~{}r_{{}_{M}}\equiv\sqrt{\frac{GM}{a_{{}_{0}}}}\,.$ (1)

The normalisation of $a_{{}_{0}}$ is taken so that the proportionality here
becomes an equality. Empirically, $a_{{}_{0}}\approx 1.2\times{10}^{-10}$ m/s2
to match galaxy RCs (e.g. McGaugh, 2011; Li et al., 2018).

It was noticed early on (e.g. Milgrom, 1983a, 1999) that, remarkably, this
value is similar to accelerations of cosmological significance. For example,

$\displaystyle
2\mathrm{\pi}a_{{}_{0}}~{}\approx~{}a_{{}_{H}}\left(0\right)~{}\equiv~{}cH_{0}~{}\approx~{}a_{{}_{\Lambda}}~{}\equiv~{}\frac{c^{2}}{\ell_{\Lambda}}\,$
(2)

where $H_{0}$ is the present-day value of the Hubble constant and
$\ell_{\Lambda}=\left(\Lambda/3\right)^{-1/2}$ is the de Sitter radius
corresponding to the observed value of the cosmological constant $\Lambda$
(Riess et al., 1998; Perlmutter et al., 1999).

Stated another way, $a_{{}_{0}}$ is similar to the acceleration at which the
classical energy density in a gravitational field (Peters, 1981, equation 9)
becomes comparable to the dark energy density
$u_{{}_{\Lambda}}\equiv\rho_{{}_{\Lambda}}c^{2}$ implied by $\Lambda$.

$\displaystyle\frac{g^{2}}{8\mathrm{\pi}G}~{}<~{}u_{{}_{\Lambda}}~{}~{}\Leftrightarrow~{}~{}g~{}\la~{}2\mathrm{\pi}a_{{}_{0}}\,.$
(3)

This association of local MOND with cosmology suggests that MOND may arise
from quantum gravity effects (e.g. Milgrom, 1999; Pazy, 2013; Verlinde, 2016;
Smolin, 2017).

Regardless of its underlying microphysical explanation, MOND correctly
predicted the RCs of a wide variety of both spiral and elliptical galaxies
across a vast range in mass, surface brightness and gas fraction (e.g.
Milgrom, 1988; Begeman et al., 1991; Sanders, 1996; Sanders & Verheijen, 1998;
Sanders & Noordermeer, 2007; Lelli et al., 2017; Li et al., 2018). Although
internal accelerations are harder to measure within ellipticals, this can
sometimes be done accurately when they have a surrounding $X$-ray emitting gas
envelope in hydrostatic equilibrium (Milgrom, 2012) or a thin rotation-
supported gas disk (e.g. Serra et al., 2012; den Heijer et al., 2015; Serra et
al., 2016). The success of MOND extends down to pressure-supported galaxies as
faint as the satellites of M31 (McGaugh & Milgrom, 2013) and the Milky Way,
though in the latter case one must be careful to exclude galaxies where MOND
predicts significant tidal distortion (McGaugh & Wolf, 2010). For a recent
overview of how well MOND works in several different types of galaxy across
the Hubble sequence, we refer the reader to Lelli et al. (2017).

It is worth emphasising that MOND does all this based solely on the
distribution of luminous matter. It is clear that these achievements are
successful a priori predictions because most of these RCs $-$ and sometimes
even the baryon distributions $-$ were measured in the decades after the first
MOND field equation was put forth (Bekenstein & Milgrom, 1984) and its new
fundamental constant $a_{{}_{0}}$ was determined (Milgrom, 1988; Begeman et
al., 1991). These predictions work due to underlying regularities in galaxy
RCs that are difficult to reconcile with the collisionless DM halos of the
$\Lambda$CDM paradigm (e.g. Salucci & Turini, 2017; Desmond, 2017a, b).

These halos were originally introduced to boost the RCs of disk galaxies. If
real, they would also endow their embedded disks with added stability because
the halo contribution to the gravitational field responds very little to
density perturbations in the disk, making the disk more like a set of test
particles. This fact forms the basis of another argument originally adduced
for the presence of DM halos around disk galaxies (Ostriker & Peebles, 1973)
$-$ without such halos, observed galactic disks would be deleteriously
unstable (Hohl, 1971).

A crucial prediction of MOND was that low surface brightness galaxies (LSBs)
would show large acceleration discrepancies at all radii because $g\ll
a_{{}_{0}}$ everywhere in the galaxy (Milgrom, 1983a, b). This prediction was
thoroughly vindicated by later observations (e.g. de Blok & McGaugh, 1997;
McGaugh & de Blok, 1998).

In $\Lambda$CDM, such LSBs must be assigned a halo much more massive than the
disk. Such a halo would cause the disk to become very stable, stymieing the
formation of bars and spiral arms which are believed to result from disk
instabilities (Lin & Shu, 1964).

Since MOND posits that dark halos are absent, the question arises regarding
the degree of stability of disks, especially LSB disks in which MOND is
predicted to have significant effects. This issue was discussed in some detail
by Milgrom (1989), who showed that MOND generally does add to the stability of
low-acceleration disks with a given mass distribution and velocity dispersion.
The reason is that instead of the Newtonian relation $g_{d}\propto\rho$
between mass density $\rho$ and the acceleration $g_{d}$ it produces in the
disk, the deep-MOND relation is $g_{d}\propto\sqrt{\rho}$. This means that in
Newtonian disks without a DM halo, density perturbations $\delta\rho$ produce
acceleration perturbations

$\displaystyle\frac{\delta g_{d}}{g_{d}}~{}\sim~{}\frac{\delta\rho}{\rho}\,.$
(4)

Adding a non-responsive DM halo with contribution $g_{h}$, we can write
$g_{d}=g_{b}+g_{h}$, where $g_{b}$ is the Newtonian contribution of the
baryonic disk that satisfies Equation 4. Because density perturbations in the
disk do not affect $g_{h}$, we have that

$\displaystyle\frac{\delta
g_{d}}{g_{d}}~{}\sim~{}\frac{\delta\rho}{\rho}\left(1+\frac{g_{h}}{g_{b}}\right)^{-1}.$
(5)

The reduced response of $g_{d}$ implies increased disk stability.

The analogous result in the deep-MOND regime is

$\displaystyle\frac{\delta
g_{d}}{g_{d}}~{}\sim~{}\frac{1}{2}\frac{\delta\rho}{\rho}$ (6)

because $g_{d}\propto\sqrt{\rho}$ and $g_{h}=0$. The added degree of stability
in deep MOND is thus similar to that endowed by a halo with $g_{h}\sim g_{b}$.
As MOND effects are strongest in this regime, it is clear that this is the
limit to how much MOND enhances the stability of disk galaxies, even those
with very low surface densities. However, the massive halos required by
$\Lambda$CDM would increase the stability indefinitely as the surface density
is reduced. This makes deep-MOND LSBs develop spiral arms and fast-rotating
bars more readily than according to $\Lambda$CDM (Tiret & Combes, 2007, 2008;
Combes, 2016).

Beyond such general semi-qualitative arguments, it is important to study disk
stability more quantitatively in specific, action-based MOND theories. This
not only gives a more accurate picture, it is also necessary for understanding
and avoiding the development of instabilities in numerical codes based on
these theories.

Milgrom (1989) studied disk stability analytically in the context of the
aquadratic Lagrangian formulation of MOND (AQUAL, Bekenstein & Milgrom, 1984).
This was followed by the numerical studies of Brada & Milgrom (1999), who
solved the AQUAL field equation using an expensive non-linear grid relaxation
stage. This is unavoidable in the context of AQUAL and remains an aspect of
more recent codes that solve it (Londrillo & Nipoti, 2009; Candlish et al.,
2015).

Since that time, another non-relativistic, action-based MOND theory has been
put forth. This quasilinear formulation of MOND (QUMOND, Milgrom, 2010) is
much more amenable to numerical simulations because the grid relaxation stage
is linear, like in Newtonian gravity. Despite using different field equations
to implement MOND consistently, QUMOND and AQUAL give rather similar results,
as demonstrated both numerically (Candlish, 2016) and analytically (Banik &
Zhao, 2018a).

QUMOND can be solved numerically using the publicly available $N$-body and
hydrodynamics code Phantom of RAMSES (Lüghausen et al., 2015), an adaptation
of the grid-based RAMSES algorithm widely used by astronomers due to its
adaptive mesh refinement feature (Teyssier, 2002).222The similar code RAyMOND
can also handle AQUAL (Candlish et al., 2015). It will eventually be made
public. As a result, QUMOND has become the main workhorse for simulations of
galaxy evolution and interactions (Thies et al., 2016; Renaud et al., 2016;
Thomas et al., 2018a, b; Bílek et al., 2018).

The question of disk stability in QUMOND remains to be addressed analytically
despite its importance for understanding the results of such simulations and
for establishing their initial conditions. Here, we study some aspects of disk
stability in QUMOND. We focus on deriving an analytic expression for the
QUMOND generalisation of the Toomre criterion (Toomre, 1964). This gives an
estimate of whether part of a disk is vulnerable to self-gravitating collapse.
In particular, the Toomre criterion gives a lower limit to the radial velocity
dispersion in a thin disk for it to remain stable against short-wavelength
axisymmetric perturbations.

After outlining the broader context in Section 1, we explain the steps in our
analytic derivation (Section 2). Our main result is a generalisation of the
Toomre condition to QUMOND (Equation 37). In Section 3, we apply this to
numerically determined RCs of a family of thin exponential disks with
different central surface density. Observational constraints on disk stability
are discussed in Section 4. We conclude in Section 5.

##  2 Stability analysis in QUMOND

QUMOND (Milgrom, 2010) is a modified gravity formulation of non-relativistic
MOND that is derivable from an action principle. It is also the non-
relativistic limit of a certain bimetric theory of relativistic MOND (Milgrom,
2009a). In QUMOND, the gravitational potential $\Phi$ is determined (given
suitable boundary conditions) by

$\displaystyle\nabla^{2}\Phi~{}=~{}\nabla\cdot\left[\nu\left(\frac{\left|\nabla\Phi_{N}\right|}{a_{{}_{0}}}\right)\nabla\Phi_{N}\right],$
(7)

where $\nu(y)$ is the so-called MOND interpolating function (in the $\nu$
version) between the Newtonian and modified regimes. $\Phi_{N}$ is an
auxiliary potential that solves the Poisson equation with the mass density
$\rho$ as its source. It is the Newtonian potential of the system, but this in
itself plays only an indirect role in the dynamics.

In this section, we determine the stability condition of a thin, flat
axisymmetric QUMOND disk to tightly wound Wentzel-Kramers-Brillouin (WKB)
axisymmetric perturbations. For this purpose, we consider perturbations
$\bar{\bm{g}}$ to the acceleration field
$-\nabla\Phi\equiv\bm{g}\equiv\bm{g}_{{}_{0}}+\bar{\bm{g}}$, where
$\bm{g}_{{}_{0}}\equiv-\nabla\Phi_{0}$ is the unperturbed or background field
(0 subscripts indicate unperturbed quantities while ${\,\bar{}\,}$
superscripts denote perturbations). In these terms, Equation 7 reads

$\displaystyle\nabla\cdot\left(\bm{g}_{{}_{0}}+\bar{\bm{g}}\right)~{}=~{}\nabla\cdot\left[\nu\left(\frac{\left|\bm{g}_{{}_{N0}}+\bar{\bm{g}}_{{}_{N}}\right|}{a_{{}_{0}}}\right)\left(\bm{g}_{{}_{N0}}+\bar{\bm{g}}_{{}_{N}}\right)\right].$
(8)

The N subscripts refer to Newtonian quantities, which must first be calculated
in order to solve Equation 8. Since perturbations in the auxiliary
intermediate Newtonian field drive those in the actual QUMOND field, we first
discuss how $\bm{g}_{{}_{N}}\equiv-\nabla\Phi_{{}_{N}}$ is affected by a
density perturbation.

###  2.1 Perturbation to the Newtonian potential

The Toomre criterion (which our work generalises) concerns a localised
(wavelength ${\ll r}$) axisymmetric perturbation to the Newtonian potential
$\Phi_{{}_{N}}$ (Toomre, 1964). This perturbation satisfies the Laplace
equation outside the disk because there is no matter there. In cylindrical
polar co-ordinates $\left(r,z\right)$ chosen so $z=0$ is the disk mid-plane,
this becomes

$\displaystyle\frac{\partial^{2}\bar{\Phi}_{{}_{N}}}{\partial
r^{2}}\,+\,\frac{\partial^{2}\bar{\Phi}_{{}_{N}}}{\partial z^{2}}~{}=~{}0\,.$
(9)

The WKB approximation made here greatly simplifies our analysis, which would
otherwise have to consider geometric corrections in addition to variation of
the background surface density within a single perturbation wavelength.
Nonetheless, use of the WKB assumption is an important shortcoming of our
work, which we discuss further in Section 3.4. Its main result is that our
analysis breaks down sufficiently close to the centre of a galaxy (Figure 5).

To analyse the stability of a disk to perturbations of different wavelengths,
we assume the surface density $\Sigma$ is perturbed with real radial wave-
vector $\alpha$ such that the density perturbation is the real part of

$\displaystyle\bar{\Sigma}~{}=~{}\widetilde{\Sigma}\,\mathrm{e}^{\mathrm{i}\alpha
r}.$ (10)

We expect that the resulting perturbation to $\Phi_{{}_{N}}$ takes the form

$\displaystyle\bar{\Phi}_{{}_{N}}~{}=~{}\widetilde{\Phi}_{{}_{N}}\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{r}}.$
(11)

where $\bm{k}\equiv\left(k_{r},k_{z}\right)$ is independent of the position
$\bm{r}$ relative to the disk centre.

Given the nature of the density perturbation $\bar{\Sigma}$ (Equation 10), we
assume $k_{r}=\alpha$. Outside the disk, Equation 9 implies that
${k_{r}}^{2}+{k_{z}}^{2}=0$, so $k_{z}$ must be imaginary. Choosing the sign
so that $\bar{\Phi}_{{}_{N}}$ decays away from the disk on both sides and
equating the discontinuity in $\frac{\partial\bar{\Phi}_{{}_{N}}}{\partial z}$
with $4\pi G\,\bar{\Sigma}$, we get that

$\displaystyle\bar{\Phi}_{{}_{N}}~{}=~{}-\frac{2\mathrm{\pi}G\widetilde{\Sigma}}{\left|\alpha\right|}\mathrm{e}^{\left(\mathrm{i}\alpha
r-\alpha\left|z\right|\right)}\,.$ (12)

In what follows, we consider the region $z>0$. Due to reflection symmetry
about the disk mid-plane, the results apply to both sides if we make the
identification $z\to\left|z\right|$. We also assume that $\alpha>0$ as local
disturbances with wave-vectors $\pm\alpha$ should behave in the same way given
the equality between the real parts of $\mathrm{e}^{\pm\mathrm{i}\alpha r}$.

###  2.2 Linearised QUMOND

For a WKB stability analysis, we consider small perturbations $\bar{\bm{g}}$
whose wavelength is much shorter than the scale over which the background
$\bm{g}_{{}_{0}}\left(\bm{r}\right)$ varies (typically the maximum of
$\left|\bm{r}\right|$ and the disk scale length). We use the governing
Equation 8 on the side $z>0$ in the region outside the disk. In this region,
$\bm{g}_{{}_{0}}$ varies only a little over many perturbation wavelengths, and
as usual in our Toomre-type analysis it can be assumed constant. As the
background fields satisfy Equation 7, we focus on the first order perturbative
terms in Equation 8.

$\displaystyle\nabla\cdot\bar{\bm{g}}~{}$ $\displaystyle=$
$\displaystyle~{}\nu_{{}_{0}}\left[\nabla\cdot\bar{\bm{g}}_{{}_{N}}+K_{0}(\widehat{\bm{w}}\cdot\nabla){\bar{g}}_{{}_{Nw}}\right],~{}\text{where}$
(13) $\displaystyle\nu_{{}_{0}}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\nu\left(\frac{g_{{}_{N0}}}{a_{{}_{0}}}\right)~{}~{}\text{and}$
(14) $\displaystyle K_{0}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\left.\frac{d\ln[\nu(y)]}{d\ln(y)}\right|_{y=g_{{}_{N0}}/a_{{}_{0}}}.$
(15)

Here, $g_{{}_{N0}}\equiv|\bm{g}_{{}_{N0}}|$ and $\widehat{\bm{w}}$ is a unit
vector in the direction of ${-\bm{g}_{{}_{0}}}$. This lets us define
${\bar{g}}_{{}_{Nw}}\equiv\bar{\bm{g}}_{{}_{N}}\cdot\widehat{\bm{w}}$ i.e.
${\bar{g}}_{{}_{Nw}}$ is the component of $\bar{\bm{g}}_{{}_{N}}$ in the
$\widehat{\bm{w}}$ direction. Equation 13 can be obtained by considering the
QUMOND dynamics of a system whose internal accelerations are much smaller than
the system’s acceleration in a constant ‘external’ field, represented here by
$\bm{g}_{{}_{0}}$ on one side of the disk. Such an external field-dominated
situation was treated in equation 67 of Milgrom (2010) and in equation 24 of
Banik & Zhao (2018a), with the latter work explaining the derivation.

###  2.3 Perturbation to the QUMOND potential

Substituting the Newtonian potential perturbation from Equation 12 into
Equation 13 gives the equation governing the QUMOND potential perturbation.

$\displaystyle\nabla^{2}\bar{\Phi}~{}$ $\displaystyle=$
$\displaystyle~{}Q\mathrm{e}^{\mathrm{i}\alpha r-\alpha z},$ (16)
$\displaystyle Q~{}$ $\displaystyle\equiv$
$\displaystyle~{}2\mathrm{\pi}G\nu_{{}_{0}}K_{0}\widetilde{\Sigma}\alpha\left(\cos
2\theta+\mathrm{i}\sin 2\theta\right),$ (17)

where $\theta$ is the angle $\widehat{\bm{w}}$ makes with the outwards radial
direction just above the disk plane.

As is well known, there are two parts to the general solution of an
inhomogeneous equation like Equation 16. The first is an arbitrary multiple of
the complementary function $\Phi_{c}$ that, in this case, corresponds to
solving the homogeneous Laplace equation with $\mathrm{e}^{\mathrm{i}\alpha
r}$ dependence at ${z=0}$. The second part of the solution is the particular
integral $\Phi_{p}$ of the inhomogeneous Equation 16. We already showed in
Section 2.1 that

$\displaystyle\Phi_{c}~{}=~{}A\mathrm{e}^{\mathrm{i}\alpha r-\alpha z},$ (18)

for any constant $A$.

The particular integral of Equation 16 must also have the same
$\mathrm{e}^{\mathrm{i}\alpha r}$ dependence at $z=0$. To exploit our
knowledge of $\Phi_{c}$, we try a solution of the form
$\Phi_{p}=f(z)\mathrm{e}^{\mathrm{i}\alpha r-\alpha z}$. Substituting this
into Equation 16, we obtain that $f^{\prime\prime}-2\alpha f^{\prime}=Q$,
where ′ indicates a radial derivative. Of the two resulting constants of
integration, one is absorbed by adjustments to $A$. The other is fixed by
requiring $\bar{\Phi}$ to decay away from the disk plane, implying that
$\Phi_{p}$ must do so. Thus, we obtain that

$\displaystyle\Phi_{p}~{}=~{}-\frac{Qz}{2\alpha}\mathrm{e}^{\left(\mathrm{i}\alpha
r-\alpha z\right)}\,.$ (19)

Noting that our derivation can be extended to both sides of the disk by
replacing $z\to\left|z\right|$, we see that the general solution is

$\displaystyle\bar{\Phi}=\mathrm{e}^{\left(\mathrm{i}\alpha
r-\alpha\left|z\right|\right)}\left[A-\mathrm{\pi}\left|z\right|\nu_{{}_{0}}K_{0}G\widetilde{\Sigma}\left(\cos
2\theta+\mathrm{i}\sin 2\theta\right)\right]\,.$ (20)

Similarly to the Newtonian case (Equation 12), $\bar{\Phi}$ decays
exponentially with $\left|z\right|$ over a perturbation wavelength, though the
decay is of the form $z\mathrm{e}^{-z}$ rather than $\mathrm{e}^{-z}$. A more
striking difference is that the QUMOND $\bar{\Phi}$ and $\bar{\Sigma}$ are not
in phase as $\sin 2\theta\neq 0$ in general. For a disk stability analysis,
only $\bar{\Phi}$ at ${z=0}$ is relevant and this does have the same phase as
$\bar{\Sigma}$. The phase difference becomes apparent only when ${z\neq 0}$.

In Newtonian gravity, such a phase offset is impossible because the problem is
symmetric with respect to ${r\to-r}$ at minima and maxima in $\bar{\Sigma}$.
This is not so in QUMOND because it depends non-linearly on the total
gravitational field, _including its background value_. This generally has some
radial component, thus picking out a preferred direction and breaking the
reflection symmetry. This unavoidable effect in MOND causes it to violate the
strong equivalence principle, a potentially testable prediction in the near
future (Pereira et al., 2016). It also underlies the so-called external field
effect in MOND (Milgrom, 1986), which is important for understanding the
internal dynamics of systems such as the Milky Way (Banik & Zhao, 2018b), its
satellite Crater 2 (McGaugh, 2016; Caldwell et al., 2017) and the more distant
NGC 1052-DF2 (Famaey et al., 2018; Kroupa et al., 2018).

###  2.4 The boundary condition

In this section, we fix the constant $A$ in Equation 20 using the boundary
condition on $\bar{\Phi}$ at the surface of the disk. For this purpose, we
apply Gauss’ theorem to the linearised QUMOND equation with constant
$\bm{g}_{{}_{0}}$ (Equation 13). This tells us that, at $z=0^{+}$, we must
have equality between the $z$ components of the vectors under the divergence
on both sides of this equation. Using z subscripts to indicate this component
of a vector, we get that

$\displaystyle{\bar{g}}_{{}_{z}}~{}=~{}\nu_{{}_{0}}\left({\bar{g}}_{{}_{Nz}}+K_{0}\,{\widehat{\bm{w}}}_{{}_{z}}{\bar{g}}_{{}_{Nw}}\right)\,.$
(21)

Upon differentiating the Newtonian potential from Equation 12, we get that

$\displaystyle{\bar{g}}_{{}_{Nw}}~{}=~{}2\mathrm{\pi}G\widetilde{\Sigma}\left(-\mathrm{i}\sin\theta+\cos\theta\right)\mathrm{e}^{\mathrm{i}\alpha
r}\,.$ (22)

As a result, the boundary condition becomes

$\displaystyle-A\alpha-\nu_{{}_{0}}\mathrm{\pi}GK_{0}\widetilde{\Sigma}\left(\cos^{2}\theta-\sin^{2}\theta+2\mathrm{i}\sin\theta\cos\theta\right)$
$\displaystyle=~{}2\mathrm{\pi}G\nu_{{}_{0}}\widetilde{\Sigma}\left[1+K_{0}\sin\theta\left(-\mathrm{i}\cos\theta+\sin\theta\right)\right]\,.$
(23)

This can be simplified to give

$\displaystyle
A~{}=~{}-\frac{\mathrm{\pi}G\nu_{{}_{0}}\widetilde{\Sigma}\left(2+K_{0}\right)}{\left|\alpha\right|}\,.$
(24)

In the Newtonian limit where $K_{0}=0$ and $\nu_{{}_{0}}=1$, this reduces to
$A=-\frac{2\mathrm{\pi}G\widetilde{\Sigma}}{\left|\alpha\right|}$, thereby
reproducing Equation 12 and equation 5.161 of (Binney & Tremaine, 2008).

###  2.5 The final result

The response of the disk to a density perturbation is governed by the
resulting potential perturbation within the disk plane. Restricting our more
general result for $\bar{\Phi}$ (Equation 20) to ${z=0}$ shows that disk
stability is governed by the parameter $A$. Comparing its value in Equation 24
with what it would be in Newtonian gravity, we get that the effective $G$
entering a local stability analysis is boosted by the factor

$\displaystyle\gamma~{}=~{}\nu_{{}_{0}}\left(1+\frac{K_{0}}{2}\right).$ (25)

Thus, the disk stability results in Binney & Tremaine (2008) can be obtained
by multiplying its equation 5.161 by $\gamma$ and following the change
through. This is because $G$ does not affect the restoring force arising from
velocity dispersion (analogous to pressure in a gas). Once the background
$\bm{g}_{{}_{0}}$ is fixed, so is the RC. In this case, $G$ also has no effect
on the restoring force from disk shear. Only the self-gravitating term (e.g.
Binney & Tremaine, 2008, equation 6.66) is affected, so $G$ only enters the
stability criterion inasmuch as it affects the relation between
$\widetilde{\Sigma}$ and $A$. Therefore, the Toomre stability condition for
Newtonian disks (Toomre, 1964) can be generalised to QUMOND by setting

$\displaystyle G~{}\to~{}G\nu_{{}_{0}}\left(1+\frac{K_{0}}{2}\right).$ (26)

Milgrom (1989) derived the corresponding result for AQUAL, whose governing
field equation is
$\nabla\cdot\left[\mu(|\bm{g}|/a_{{}_{0}})\bm{g}\right]=\nabla\cdot\bm{g}_{{}_{N}}$.
In this theory, the analogue of Equation 26 is

$\displaystyle G$ $\displaystyle\to$
$\displaystyle\frac{G}{\mu_{{}_{0}}\sqrt{1+L_{0}}},~{}\text{ where}$ (27)
$\displaystyle\mu_{{}_{0}}$ $\displaystyle\equiv$
$\displaystyle\mu\left(\frac{g_{{}_{0}}}{a_{{}_{0}}}\right)~{}\text{ and}$
(28) $\displaystyle L_{0}$ $\displaystyle\equiv$
$\displaystyle\left.\frac{d\ln\left[\mu\left(x\right)\right]}{d\ln\left(x\right)}\right|_{x=g_{{}_{0}}/a_{{}_{0}}}~{}~{}.$
(29)

The results in AQUAL and QUMOND are rather similar $-$ when we approximate
each theory by its algebraic analogue (exactly valid in spherical symmetry),
QUMOND reduces to
$\bm{g}=\nu\left(\frac{\left|\bm{g}_{{}_{N}}\right|}{a_{{}_{0}}}\right)\bm{g}_{{}_{N}}$
and AQUAL to
$\mu\left(\frac{\left|\bm{g}\right|}{a_{{}_{0}}}\right)\bm{g}=\bm{g}_{{}_{N}}$.
In this one-dimensional situation, the theories are equivalent provided
$\nu\left(\frac{\left|\bm{g}_{{}_{N}}\right|}{a_{{}_{0}}}\right)\mu\left(\frac{\left|\bm{g}\right|}{a_{{}_{0}}}\right)=1$.
Assuming this to be the case, Banik & Zhao (2018d) showed in their section 7.2
that this pair of $\bm{g}$ and $\bm{g}_{{}_{N}}$ satisfy

$\displaystyle\left(1+K_{0}\right)\left(1+L_{0}\right)~{}=~{}1.$ (30)

To first order in $K_{0}$ or $L_{0}$, this implies equality between the
factors of $\nu_{{}_{0}}\left(1+\frac{K_{0}}{2}\right)$ in Equation 26 and
$\frac{1}{\mu_{{}_{0}}\sqrt{1+L_{0}}}$ in Equation 27. The disk stability
condition is therefore very similar in both formulations of MOND. This
suggests that simulations of disk galaxies in the more computer-friendly
QUMOND should yield rather similar results to AQUAL (Candlish et al., 2015).

In general, Equation 30 is not exactly correct because the algebraic relations
are only approximations. However, we show in Section 3.5 that the QUMOND
version is very accurate beyond the central scale length of a thin exponential
disk galaxy. Similar results were previously obtained for AQUAL (Brada &
Milgrom, 1995).

##  3 Disk stability and surface density

Disk galaxies generally have an exponential surface density profile (Freeman,
1970). In this section, we consider infinitely thin exponential disks with a
range of parameters. Because MOND is an acceleration-dependent theory, the
only parameter we need to vary is the central surface density $\Sigma_{0}$. We
therefore consider a family of exponential disks with different $\Sigma_{0}$.
Our results hold for any disk scale length $r_{{}_{d}}$ as long as the total
mass is varied $\propto{r_{{}_{d}}}^{2}$. Otherwise, we must consider a
different member of our disk family with the appropriate $\Sigma_{0}$.

In this and subsequent sections, we consider only the unperturbed
gravitational field in disk galaxies. We therefore drop the convention that 0
subscripts refer to background quantities. Instead, all quantities are
understood to be background values.

###  3.1 Rotation curves

The first step to finding $\Phi$ is obtaining the Newtonian potential
$\Phi_{{}_{N}}$. This satisfies Equation 9 outside the disk as there is no
matter there. The presence of the disk imposes a boundary condition at all
points with $z=0$. Applying the Poisson equation to a ‘Gaussian pill-box’
around a small part of the disk shows that $\Phi_{{}_{N}}$ must satisfy

$\displaystyle\left.\frac{\partial\Phi_{{}_{N}}}{\partial
z}\right|_{z=0^{+}}~{}=~{}2\pi G\Sigma~{}=~{}2\pi
G\Sigma_{0}\rm{e}^{-\frac{r}{r_{{}_{d}}}}.$ (31)

We solve this using grid relaxation of Equation 9, focusing on the region
$z>0$ because of symmetry. To speed up the numerical convergence, we use
successive over-relaxation on a single grid with a customised radial
resolution. The details of our procedure are explained in appendix A of Banik
& Zhao (2018b), where we also mention our initial guess for $\Phi_{{}_{N}}$,
the boundary and stopping conditions and the choice of over-relaxation
parameter.

Once we have found $\bm{g}_{{}_{N}}$, we obtain $\nabla\cdot\bm{g}$ using
Equation 7. The exact result will depend on the assumed MOND interpolating
function, for which we use the ‘simple’ form (Milgrom, 1986; Famaey & Binney,
2005). Observational reasons for preferring this form were discussed in
section 7.1 of Banik & Zhao (2018d).

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

We obtain $\bm{g}$ from its divergence using direct summation done similarly
to section 2.2 of Banik & Zhao (2018b). This exploits the axisymmetric nature
of the problem using a ring library procedure. The situation is simpler here
because we do not consider any external field on the galaxy. This removes the
need to apply their equation 11, which is in any case incorrect (Banik & Zhao,
2018c).

Figure 1: Rotation curves of thin exponential disk galaxies in QUMOND with
different central surface density $\Sigma_{0}$, parameterised by $S$ in
Equation 33. Velocities are shown relative to the flatline level
$v_{{}_{f}}=\sqrt[4]{GMa_{{}_{0}}}$ for galaxies with scale length
$r_{{}_{d}}$ and total mass $M=2\pi{r_{{}_{d}}}^{2}\Sigma_{0}$. Our adopted
values for $S$ are listed in Table 1. The peak of each curve is indicated with
a pink star. The upper cyan curve has $S=0.06$ (similar to the Milky Way)
while the lower cyan curve has $S=0.15$ (similar to M31). To aid visibility,
curves with different $S$ are shown alternating between solid and dot-dashed
styles. The result for the deep-MOND limit ($S\to\infty$) is shown as a solid
black curve lying below the other curves. The two thin dashed blue curves show
the estimated $v_{c}\left(r\right)$ for $S=0.05$ and the deep-MOND limit using
the algebraic MOND relation (Equation 55). These agree very closely with the
corresponding numerical QUMOND curves.

We show our family of RCs in Figure 1, normalised according to the flatline
level $v_{{}_{f}}$ of each curve. The different RCs have different values of
the parameter $S$, which governs the importance of MOND to the galaxy. $S$ can
be thought of as the ratio between $\Sigma_{0}$ and the critical MOND surface
density $\Sigma_{M}$ (Milgrom, 2016).

$\displaystyle S~{}\equiv~{}\overbrace{\frac{a_{{}_{0}}}{2\pi
G}}^{\Sigma_{M}}\,\div\,\Sigma_{0}\,.$ (33)

Values of $S$
---
0.025 | 0.05
0.06 | 0.075
0.1 | 0.15
0.2 | 0.3
0.4 | 0.6
0.75 | 1
2 | 4
10 | $\infty$
Table 1: Values of the surface density parameter $S$ (Equation 33) in thin
exponential disk galaxies for which we show RCs and disk stability criteria in
our figures. $S\to\infty$ corresponds to the deep-MOND limit and $S\to 0$ to
the Newtonian limit.

In this way, we find the outwards radial gravity $g_{{}_{r}}$ at a finite
number of points along a radial transect within the disk plane. This allows us
to obtain the rotation speed

$\displaystyle v_{c}\left(r\right)~{}=~{}\sqrt{-rg_{{}_{r}}}\,.$ (34)

The curves shown in Figure 1 correspond to galaxies where $S$ takes the values
listed in Table 1. Galaxies with lower $S$ have a higher $v_{c}\left(r\right)$
that flattens out at larger $r$. We alternate the line styles between solid
and dashed to help identify them. The lowest black curve is the result for the
deep-MOND limit ($S\to\infty$), which we obtain by using
$\nu=\sqrt{\frac{a_{{}_{0}}}{\left|\bm{g}_{{}_{N}}\right|}}$ in Equation 7.

###  3.2 Minimum $\sigma_{r}$

The QUMOND RCs derived in Section 3.1 allow us to determine the minimum
$\sigma_{r}$ required for local stability of a stellar disk via application of
the Toomre stability criterion (Toomre, 1964), with $G$ adjusted according to
Equation 26.

$\displaystyle\sigma_{r}~{}$ $\displaystyle\geq$
$\displaystyle~{}\frac{3.36\,G\nu\Sigma\left(1+\frac{K_{0}}{2}\right)}{\Omega_{r}},~{}~{}\text{
where}$ (35) $\displaystyle{\Omega_{r}}^{2}~{}$ $\displaystyle=$
$\displaystyle~{}-\frac{3g_{{}_{r}}}{r}\,-\,\frac{\partial
g_{{}_{r}}}{\partial r}\,.$ (36)

Equation 35 can be written as a constraint on the so-called Toomre $Q_{*}$
parameter for stellar disks (Binney & Tremaine, 2008, equation 6.71).

$\displaystyle
Q_{*}~{}\equiv~{}\frac{\sigma_{r}\Omega_{r}}{3.36\,G\nu\Sigma\left(1+\frac{K_{0}}{2}\right)}~{}\geq~{}1$
(37)

In exactly the same way, we can generalise the corresponding result for
isothermal gas disks with sound speed $c_{s}$ (Binney & Tremaine, 2008,
equation 6.68).

$\displaystyle
Q_{gas}~{}\equiv~{}\frac{c_{s}\Omega_{r}}{\rm{\pi}\,G\nu\Sigma\left(1+\frac{K_{0}}{2}\right)}~{}\geq~{}1$
(38)

Our results in Figure 2 show that high surface density galaxies (with low $S$)
need rather high $\sigma_{r}$ in their central regions. Thus, MOND is unable
to stabilise high surface brightness (HSB) rotationally supported disks: doing
so would need such a large $\sigma_{r}$ that the disk would be dispersion-
supported. This is in line with the fact that a low $S$ implies MOND is not
very relevant to the galaxy, making it similar to a purely Newtonian
exponential disk lacking a DM halo. Such systems are dynamically unstable
(Hohl, 1971; Ostriker & Peebles, 1973).

Observationally, disk galaxies do indeed have an upper limit to their central
surface brightness. Due to the difficulty of detecting LSB galaxies, this was
first noticed in terms of all disk galaxies having only a very narrow range in
surface brightness (Freeman, 1970). It was later realised that the paucity of
LSB galaxies was likely a selection effect due to the brightness of the night
sky (Disney, 1976). However, higher surface brightness galaxies should be even
easier to detect against the sky, suggesting that the paucity of such galaxies
is a real aspect of our Universe (van der Kruit, 1987; McGaugh, 1996). The
maximum surface brightness is similar to that expected in MOND for a
reasonable mass to light ratio (Brada & Milgrom, 1999). More recent
observations confirm the presence of an upper limit (Fathi, 2010).

Figure 2: The minimum radial velocity dispersion $\sigma_{r}$ required for
local stability of our family of thin exponential QUMOND disks, using the
criterion in Equation 35. We use a similar style to Figure 1 and show results
for the same central surface densities (values of $S$ listed in Table 1).

In Figure 3, we show the minimum $\sigma_{r}$ profiles for $\Lambda$CDM-like
models in which the self-gravity of the disks are governed by Newtonian
gravity but their RCs follow the MOND predictions. This is because empirical
RCs follow MOND expectations very closely across a huge range of galaxy
surface brightness, size and mass (Lelli et al., 2017; Li et al., 2018). In a
$\Lambda$CDM context, this is due to properties of the DM halo. Although the
halo affects $v_{c}\left(r\right)$ and thus $\Omega_{r}$, it has no effect on
the perturbation to $\bm{g}_{{}_{N}}$ arising from a perturbation to $\Sigma$.
Consequently, we can use the classical Toomre condition (Toomre, 1964) to
analyse the stability of such disks, albeit with modified $\Omega_{r}$ for
observational consistency.

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

Figure 3: Similar to Figure 2, but for the equivalent $\Lambda$CDM disks.
These are assumed to have the same RC as the corresponding QUMOND model
(Figure 1). The stability condition follows from applying Equation 39 to these
RCs, implying no lower limit to $\sigma_{r}$ at very low surface density. Such
disks are essentially just test particles held together by the DM halo.

###  3.3 Minimum disk thickness

If we assume that the vertical velocity dispersion $\sigma_{z}=\sigma_{r}$ and
that the disk can be approximated locally as an isothermal slab, then we can
visualise our results as a minimum root mean square (rms) thickness to the
disk, $z_{{}_{rms}}$. A stratified (horizontally uniform) slab of isothermal
gas with sound speed $\sigma$ has a density $\rho$ which depends only on
height $z$ from the mid-plane at ${z\equiv 0}$, where the density is maximal.
Due to the symmetry of the slab, the Poisson equation becomes

$\displaystyle 4\mathrm{\pi}G\rho~{}=~{}-g^{\prime}\,.$ (40)

Here, we use a ′ to denote a derivative with respect to $z$ and $g$ to denote
$g_{{}_{z}}$ in this one-dimensional situation. Using this notation, the
equation of hydrostatic equilibrium becomes

$\displaystyle\sigma^{2}\rho^{\prime}~{}=~{}\rho g\,.$ (41)

The governing equations for the slab are thus of the form $\rho\propto
g^{\prime}$ and $\rho^{\prime}\propto\rho g$. These can be solved by
multiplying them together and cancelling the common factor $\rho$, yielding a
constraint of the form $\rho^{\prime}\propto\left(g^{2}\right)^{\prime}$.
Integrating this gives a relation between $\rho$ and $g$, with the additive
constant fixed by noting that $g\to\mp 2\pi G\Sigma$ and $\rho^{\prime}\to 0$
as $z\to\pm\infty$. Using Equation 40, we can express our result as a
differential equation for $g$ alone. This can be solved by separation of
variables, subject to the boundary condition ${g=0}$ at $z=0$ due to symmetry.
Given a solution for $g\left(z\right)$, it can be differentiated to obtain
$\rho$ (Equation 40). In this way, Spitzer (1942) found the density profile of
a Newtonian isothermal slab.

$\displaystyle\rho~{}$ $\displaystyle=$
$\displaystyle~{}\frac{\mathrm{\pi}G\Sigma^{2}}{2\sigma^{2}}\operatorname{sech}^{2}\left(\frac{z}{h}\right),$
(42) $\displaystyle h~{}$ $\displaystyle=$
$\displaystyle~{}\frac{\sigma^{2}}{\mathrm{\pi}G\Sigma}\,.$ (43)

The rms thickness $z_{{}_{rms}}$ of this $\operatorname{sech}^{2}$ profile is
$\frac{\mathrm{\pi}h}{\sqrt{12}}$ (Appendix A). As the vertical gravity just
outside the disk plane would be enhanced by a factor of $\nu$ in QUMOND, we
assume our Newtonian $z_{{}_{rms}}$ can be generalised to QUMOND if we
multiply the local value of $G$ by $\nu$.

$\displaystyle
z_{{}_{rms}}~{}=~{}\frac{{\sigma_{r}}^{2}}{G\nu\Sigma\sqrt{12}}~{}.$ (44)

Our constraint on $z_{{}_{rms}}$ is only meant to be illustrative as our
assumptions may not hold, in particular that $\sigma_{z}=\sigma_{r}$. Even so,
it provides a useful guide to how dynamically hot the disk must be at
different radii. Our results are shown in the top panel of Figure 4. The thick
pink curve near the top is the result for a bare Newtonian exponential disk,
whose minimum $\sigma_{r}$ follows from Equation 39 applied to the Newtonian
RC.

An interesting aspect of our results is what they reveal about disk stability
in the deep-MOND limit ${S\to\infty}$ (solid black curve below the other
curves). Even in this limit, a galaxy can only gain a limited amount of extra
stability compared to the purely Newtonian case. To understand this, we
consider how changes in $\bm{g}$ and $\bm{g}_{{}_{N}}$ relate to changes in
the surface density $\Sigma$. In Newtonian gravity, we get that

$\displaystyle\frac{\Delta\left|\bm{g}_{{}_{N}}\right|}{\left|\bm{g}_{{}_{N}}\right|}~{}\approx~{}\frac{\Delta\Sigma}{\Sigma}\,.$
(45)

Figure 4: Our results on the minimum $\sigma_{r}$ required for disk stability
(Equation 35), shown here as an equivalent rms thickness obtained using
Equation 44. The thick pink curve is the purely Newtonian result for a galaxy
lacking DM ($S\to 0$) and thus with a Keplerian RC at large radii. _Top_ :
Results for QUMOND. _Bottom_ : The analogous Newtonian results, found assuming
$\nu=1$ and $K_{0}=0$. Except for the thick pink curve, we assume a DM halo is
present such that the RCs are given by QUMOND, making them the same for both
panels. We assume the halo does not respond to the disk.

In QUMOND, this is still true but $\bm{g}\neq\bm{g}_{{}_{N}}$ as the fields
are related by Equation 8. Thus, we expect that

$\displaystyle\frac{\Delta\left|\bm{g}\right|}{\left|\bm{g}\right|}~{}$
$\displaystyle\approx$
$\displaystyle~{}\frac{\Delta\left|\bm{g}_{{}_{N}}\right|}{\left|\bm{g}_{{}_{N}}\right|}\,+\,\frac{\Delta\nu}{\nu}$
(46) $\displaystyle=$
$\displaystyle~{}\frac{\left(1+K_{0}\right)\Delta\left|\bm{g}_{{}_{N}}\right|}{\left|\bm{g}_{{}_{N}}\right|}\,.$
(47)

Because ${K_{0}<0}$, the MOND gravitational field is less perturbed than the
Newtonian one by a change to $\Sigma$. However, this extra stability is
limited because $K_{0}>-\frac{1}{2}$, making the factor of
$\left(1+K_{0}\right)$ at least $\frac{1}{2}$.

Similar results would be obtained in AQUAL, whose governing equation
$\nabla\cdot\left(\mu\bm{g}\right)=-4\mathrm{\pi}G\rho$ implies that

$\displaystyle\frac{\left(1+L_{0}\right)\Delta\left|\bm{g}\right|}{\left|\bm{g}\right|}~{}$
$\displaystyle\approx$
$\displaystyle~{}\frac{\Delta\left|\bm{g}_{{}_{N}}\right|}{\left|\bm{g}_{{}_{N}}\right|}\,.$
(48)

In this case, $L_{0}<1$, limiting the factor of $\left(1+L_{0}\right)$ to at
most 2. Roughly speaking, this means that self-gravitating disks in AQUAL or
QUMOND can never be more than twice as stable as a Newtonian disk. Either
modification to Newtonian gravity endows disks with a similar amount of extra
stability because $\left(1+L_{0}\right)\left(1+K_{0}\right)\approx 1$
(Equation 30).

In the bottom panel of Figure 4, we show how our results on $z_{{}_{rms}}$
would change for $\Lambda$CDM-like galaxy models whose stability is governed
by Equation 39 but with RCs that follow QUMOND predictions. For very high
values of $S$, the gravitational field must then be dominated by the DM halo.
Consequently, the disk can be arbitrarily cold dynamically and yet remain
stable. As just discussed, this is not the case in either version of MOND that
we consider.

###  3.4 The critical wavelength

Our analysis assumes that the radius $r$ greatly exceeds the wavelength
$\lambda_{crit}$ most unstable to a perturbation. In this section, we explore
the validity of this WKB approximation.

According to equation 22 of Toomre (1964),

$\displaystyle\lambda_{crit}~{}=~{}0.55\times\frac{4\pi^{2}G\Sigma}{{\Omega_{r}}^{2}}$
(49)

The factor of 0.55 was derived numerically in their section 5c under the
assumption that $\sigma_{r}$ satisfies Equation 39. Our work shows that the
appropriate MOND generalisation of Equation 49 is

$\displaystyle\lambda_{crit}~{}=~{}\frac{2.2\pi^{2}G\Sigma\nu\left(1+\frac{K_{0}}{2}\right)}{{\Omega_{r}}^{2}}$
(50)

At large $r$, we expect the RC to flatline such that $\Omega_{r}\propto 1/r$
so that $\lambda_{crit}\propto r^{2}\mathrm{e}^{-r}$. This proves that the WKB
approximation should be very accurate when ${r\gg r_{{}_{d}}}$ but is likely
to break down when ${r\ll r_{{}_{d}}}$.

Toomre (1964) used numerical experiments to show that the WKB approximation
works quite well even when $\lambda_{crit}\approx r$ (see their section 4b). A
likely explanation is that $\bar{\bm{g}}$ at any point mostly arises from
material within $\frac{1}{4}$ of a perturbation wavelength (Safronov, 1960).
This is due to the steep inverse square law of Newtonian gravity.
Perturbations in MOND also follow an inverse square law in both the AQUAL and
QUMOND formulations, though there is an additional angular dependence which is
absent in Newtonian gravity (Banik & Zhao, 2018a). Thus, we assume that our
WKB approximation should work well when $\lambda_{crit}<r$.

We use Figure 5 to show the results of Equations 49 and 50, each time showing
a line of equality for the reasons just discussed. The bottom panel shows
$\lambda_{crit}$ for $\Lambda$CDM-like models, where the RC is
indistinguishable from MOND. For comparison, we also use a thick pink curve to
show the result for a bare Newtonian disk without any DM halo. In Newtonian
models with a halo, the boost to the rotation curve increases ${\Omega_{r}}$,
thus reducing $\lambda_{crit}$. The reduction is more significant for a galaxy
with lower $\Sigma_{0}$ (higher $S$) because such galaxies need a more
substantial DM halo in order to explain the observed properties of LSBs.
Consequently, the overall appearance of Figure 5 is very similar to Figure 4,
our estimate of the minimum $z_{{}_{rms}}$ required by the disk to ensure its
stability.

This can be understood more quantitatively by substituting Equation 35 into
Equation 44 and comparing the result to Equation 50. In both cases, the result
is $\propto\frac{G\nu\Sigma}{{\Omega_{r}}^{2}}$. However, the constant of
proportionality is $0.15\left(1+\frac{K_{0}}{2}\right)\times$ smaller for
$z_{{}_{rms}}$ than for $\lambda_{crit}$. As a result,
$\lambda_{crit}=8.9\,z_{{}_{rms}}$ in the deep-MOND limit.

Figure 5: Our results on the wavelength most unstable to self-gravitating
collapse, calculated using Equation 50 for MOND models (_top_) and Equation 49
for Newtonian $\Lambda$CDM-like models (_bottom_). The black line of equality
is shown in both panels.

In the MOND case, the extra factor of $\nu\left(1+\frac{K_{0}}{2}\right)$
entering into a stability analysis partially counters the boost to
${\Omega_{r}}$. Nonetheless, our results in the top panel of Figure 5 indicate
that $\lambda_{crit}$ is still smaller than for a bare Newtonian disk. To
understand this, suppose that the circular orbit frequency is ${\Omega_{c}}$.
It can be straightforwardly shown that

$\displaystyle{\Omega_{r}}^{2}~{}$ $\displaystyle=$
$\displaystyle~{}\left(n+3\right){\Omega_{c}}^{2}$ (51) $\displaystyle n~{}$
$\displaystyle\equiv$
$\displaystyle~{}\frac{\partial\ln\left|g_{r}\right|}{\partial\ln r}\,.$ (52)

Here, $n$ is the logarithmic radial derivative of $\left|g_{r}\right|$, the
magnitude of the radial gravity. Combining this relation with the fact that
MOND approximately boosts $g_{r}$ by a factor of $\nu$ (Section 3.5), we can
write Equation 50 in terms of the Newtonian angular frequency
${\Omega_{c,N}}$.

$\displaystyle\lambda_{crit}~{}=~{}\frac{2.2\pi^{2}G\Sigma\nu\left(1+\frac{K_{0}}{2}\right)}{\left(n+3\right)\nu{\Omega_{c,N}}^{2}}$
(53)

The factors of $\nu$ cancel between the numerator and denominator. This is
because $\nu$ enhances both the global RC and the local response of $\bm{g}$
to a density perturbation.333A similar cancellation is also evident in
Equations 47 and 48. As a result, $\lambda_{crit}$ does not differ too greatly
from the Newtonian result in the absence of DM.

Even so, some differences do exist because we are left with an extra factor of
$\left(1+\frac{K_{0}}{2}\right)<1$ in the numerator of Equation 53. In the
deep-MOND regime, this factor is $\frac{3}{4}$. Regardless of the central
$\Sigma_{0}$, this regime is always reached in the outskirts of a galaxy.

Some role is also played by the factor of ${\left(n+3\right)}$. By definition,
${n=-2}$ in the outskirts of a bare Newtonian disk. In MOND, Equation 1
implies that the analogous result is ${n=-1}$. This difference reduces
$\lambda_{crit}$ by another factor of 2 compared to the Newtonian case.

Therefore, our results indicate that the WKB approximation should be even more
accurate for MOND than for bare Newtonian disks of the sort analysed by Toomre
(1964). Comparing our $\lambda_{crit}$ curves with the lines of equality in
Figure 5, it is clear that the least stable wavelength is much smaller than
the radius beyond the central ${\approx 2r_{{}_{d}}}$. Within this region, we
expect our analysis to be less accurate. The exact details depend on $S$: the
WKB approximation remains valid down to lower $r$ for galaxies with larger
$S$.

This effect is much stronger for $\Lambda$CDM-like models (bottom panel of
Figure 5). Thus, the $\Lambda$CDM version of our analysis should be valid
almost everywhere within a LSB galaxy. This is fortunate given the importance
of LSBs in distinguishing between $\Lambda$CDM and MOND (Section 4).

###  3.5 Comparison with algebraic MOND relations

We can compare our numerically determined QUMOND RCs for thin exponential
disks against two algebraic MOND relations (ALMs) between the Newtonian and
MOND accelerations. One such ALM relates the MOND acceleration
$g^{{}^{A}}_{{}_{r}}$ to the Newtonian $g_{{}_{Nr}}$ in the disk mid-plane.

$\displaystyle
g^{{}^{A}}_{{}_{r}}~{}=~{}\nu\left(\frac{g_{{}_{Nr}}}{a_{{}_{0}}}\right)g_{{}_{Nr}}\,.$
(54)

This is the original formulation of MOND (Milgrom, 1983a, equation 2). When
applied to RC analyses, it implies a unique relation between the accelerations
predicted by Newtonian gravity and the factor by which observed accelerations
exceed this prediction. This ‘mass discrepancy-acceleration relation’ (MDAR)
has been used in most subsequent MOND analyses of RCs, for example in the
recent detailed analysis of Li et al. (2018).444They use the term ‘radial
acceleration relation’ (RAR) as the underlying cause may not be missing mass.
In modified-inertia interpretations of MOND, this equation is exactly correct
for the mid-plane accelerations and can thus be used for making exact RC
predictions (Milgrom, 1994, 2011).

For a modified gravity interpretation of MOND, a more appropriate ALM would
equate the argument of $\nu$ with the total Newtonian acceleration just
outside the disk, not the mid-plane one where ${g_{{}_{Nz}}=0}$ (Brada &
Milgrom, 1995).

$\displaystyle g^{{}^{A}}_{{}_{r}}~{}$ $\displaystyle=$
$\displaystyle~{}\nu\left(a_{{}_{0}}^{-1}\sqrt{{g_{{}_{Nr}}}^{2}+{g_{{}_{Nz}}}^{2}}\right)g_{{}_{Nr}},~{}\text{
where}$ (55) $\displaystyle g_{{}_{Nz}}~{}$ $\displaystyle=$
$\displaystyle~{}\mp 2\pi G\Sigma\,.$ (56)

Whichever specific version of MOND one uses, both forms of the ALM are
equivalent and exact in cases of spherical symmetry, though they differ in
their application to disk galaxies. For example, Brada & Milgrom (1995) showed
in their section 3 that the ALM of Equation 55 coincides exactly with AQUAL
RCs for a Kuzmin disk, but the ALM of Equation 54 does not.555For an isolated
Kuzmin disk, the gravitational field in any theory is that of a point mass
offset from the galactic centre. Thus, Equation 55 holds exactly everywhere
outside the disk. Moreover, $g_{{}_{r}}$ is continuous vertically across the
disk.

In Figure 1, we showed that RCs obtained with Equation 55 are very similar to
those based on QUMOND in the two cases considered (dashed dark blue curves),
even for points quite close to the disk centre. This is similar to the results
obtained by Angus et al. (2012) and Jones-Smith et al. (2018).

Figure 6: The fractional error $g_{{}_{r}}\div g^{{}^{A}}_{{}_{r}}-1$ in the
radial gravity $g_{{}_{r}}$ that arises from using the algebraic approximation
to QUMOND (Equation 55). Results are shown down to $0.32\,r_{{}_{d}}$ as a
purely exponential profile is unlikely to remain accurate at arbitrarily low
radii. A very high surface brightness galaxy ($S\to 0$) would be purely
Newtonian, making the ALM exactly correct. Such a system would appear on this
graph as a flat line at 0.

We now consider in more detail the fractional difference in $g_{{}_{r}}$
between Equation 55 and QUMOND, following on from previous calculations for
AQUAL (Milgrom, 1986; Brada & Milgrom, 1995). For this purpose, we use Figure
6 to show $g_{{}_{r}}\div g^{{}^{A}}_{{}_{r}}-1$. Evidently, this ALM differs
very little from QUMOND beyond the central $0.5\,r_{{}_{d}}$, even in the
deep-MOND limit. In the Newtonian limit, MOND has no effect on the dynamics,
making the ALM in either form an exact representation of MOND.

The ALM of Equation 54 differs more significantly from QUMOND (Figure 7).
Consequently, it is important to use Equation 55 rather than Equation 54 if
one is trying to approximate QUMOND with a simple algebraic relation. Even
Equation 55 is not exactly equivalent to QUMOND, but our results suggest that
it is rather close in situations with a high degree of symmetry.

Figure 7: Similar to Figure 6, but using a different version of the ALM that
neglects the vertical gravity when calculating $\nu$ (Equation 54). The
resulting estimates of the radial gravity now differ much more from our
numerical QUMOND calculations.

Given the increasing accuracy of observations, they may be able to distinguish
between the ALMs presented here. A recent analysis comparing these ALMs with
observations favoured Equation 55, suggesting that MOND arises from a
modification of gravity (Frandsen & Petersen, 2018). Some works even calculate
RCs using a rigorous solution of QUMOND obtained with a Poisson solver (Angus
et al., 2012, 2015). If such predictions more closely agree with observations
than a simple application of Equation 54, then a modified gravity
interpretation of MOND would be favoured. Of crucial importance to such a test
would be the region ${r\approx\left(1-2\right)r_{{}_{d}}}$, though accurate
observations at even smaller radii would also be very helpful (Figure 7).

##  4 Observational context

In the $\Lambda$CDM picture, the observed internal kinematics of LSBs imply
that they must be surrounded by dominant DM halos (e.g. McGaugh & de Blok,
1998). Such inert halos would lend strong stabilising support to LSB disks,
allowing them to remain dynamically stable with a very low velocity dispersion
(Figure 3). Observed LSBs have rather higher velocity dispersions and thus
appear to be dynamically overheated (Saburova, 2011). If so, it would be
difficult for them to sustain spiral density waves, the leading explanation
for observed spiral features in HSBs (Lin & Shu, 1964). Interestingly, LSBs
also have spiral features (McGaugh et al., 1995). It has been argued that this
and other features of LSBs suggest that their gravitating mass mostly resides
in their disk, contradicting $\Lambda$CDM expectations (McGaugh & de Blok,
1998, section 3.3).

Assuming the density wave theory for spiral structures, counting the number of
spiral arms gives an idea of the critical wavelength most unstable to
amplification by disk self-gravity. Indeed, D’Onghia (2015) performed analytic
calculations for the number of spiral arms in galaxies observed as part of the
DiskMass survey (Bershady et al., 2010). She found good agreement with
observations if she made reasonable assumptions on the mass to light ratio.
The DiskMass survey ‘selects against LSB disks’, making the work of D’Onghia
(2015) an important check on the validity of the approach in a regime where
$\Lambda$CDM and MOND predictions do not greatly differ.

The theory should also apply to LSBs, which provide a good opportunity to make
a priori predictions because such galaxies were generally not known about in
the 1960s. Our results in Section 3.4 show that, in a $\Lambda$CDM context,
the critical wavelength for LSB galaxies is expected to be much shorter than
the radius almost everywhere (bottom panel of Figure 5). Thus, the WKB
approximation should work particularly well in Newtonian LSB disks with
massive DM halos.

Applying a similar spiral-counting technique to LSBs in a Newtonian context,
Fuchs (2003) found that their disks would be much too stable to allow the
formation of their observed spiral arms if their stellar masses are similar to
those suggested by stellar population synthesis models (e.g. Bell & de Jong,
2001). The spiral structure could be explained in $\Lambda$CDM only if much of
the mass needed to explain their elevated RCs resides not in a stabilising DM
halo but rather within the disk itself. This would make the disk very massive,
sometimes requiring a mass to light ratio ${>10\times}$ the Solar value in the
$R$-band.

This can be understood by considering the terms in Equation 39. An elevated RC
implies a high epicyclic frequency ${\Omega_{r}}$, making the observed
$\sigma_{r}$ much higher than the minimum required for stability if we assume
the disk has a conventional mass to light ratio. Physically, this arises
because density perturbations in such a system would wind up quickly. This
makes the disk very stable, making it difficult to form spiral arms. In
Newtonian gravity, their existence implies a much higher $\Sigma$.

An unusually massive disk is also required by the Newtonian analysis of Peters
& Kuzio de Naray (2018) to explain the pattern speeds of bars in LSBs, which
are faster than expected in 3 of the 4 galaxies they considered. Similar
results were obtained by Algorry et al. (2017) upon comparing a larger sample
of observed galaxies (Corsini, 2011; Aguerri et al., 2015) with results from
the EAGLE hydrodynamical simulations (Schaye et al., 2015; Crain et al.,
2015). Higher numerical resolution is required to be more certain of this
result.

To some extent, bars and spiral features in galaxies can be triggered by
interactions with satellites (Hu & Sijacki, 2018). However, without disk self-
gravity, any spirals formed in this way would rapidly wind up and decay due to
differential rotation of the disk (Fall & Lynden-Bell, 1981, page 111). Even
in a galaxy like M31, the simulations of Dubinski et al. (2008) indicate that
interactions with a realistic satellite population only cause mild disk
heating in excess of that which arises in a control simulation without
satellites. Thus, evidence has been mounting over several decades that the
gravity in a LSB generally comes from its disk. This contradicts the
$\Lambda$CDM expectation that it should mostly come from its near-spherical
halo of DM given the large acceleration discrepancy at all radii in LSBs.

It has been claimed that the degree of dynamical stability observed in some
galaxies appears more consistent with $\Lambda$CDM expectations than with MOND
(Sánchez-Salcedo et al., 2016). Their analysis found that NGC 6503 should have
a bar in MOND but currently did not. In fact, this galaxy does have a faint
end-on bar (Kuzio de Naray et al., 2012). Given that bar strengths are
expected to change with time, we might merely be observing its bar at a time
when it is weak. Moreover, the analysis of Sánchez-Salcedo et al. (2016) did
not get a good MOND fit to the RC of this galaxy by assuming a distance of 5.2
Mpc and allowing a 15% uncertainty (see their section 5.1). Recently, Li et
al. (2018) showed in their figure A1 that a good MOND fit can be obtained for
a distance of 6.5 Mpc, outside this range but rather consistent with the 6.25
Mpc measured by Tully et al. (2013). Given the importance of the RC to the
stability of a galaxy, there is clearly some doubt regarding the tension
claimed by Sánchez-Salcedo et al. (2016) between the MOND-predicted and
observed properties of NGC 6503. At a larger heliocentric distance, the same
rotation velocity implies a lower acceleration at the analogous position in
the galaxy e.g. at its half-light radius. In MOND, the lack of a DM halo means
this is only possible if the disk has a lower surface density. This would take
it deeper into the MOND regime, endowing it with more stability (larger $S$ in
Figure 2).

##  5 Conclusions

We have generalised the Toomre disk stability condition (Toomre, 1964) to the
quasilinear formulation of Modified Newtonian Dynamics (QUMOND, Milgrom,
2010). Our main result is Equation 37. We use this to estimate the minimum
radial velocity dispersion $\sigma_{r}$ required by thin exponential disk
galaxies to avoid local self-gravitating collapse. In Newtonian gravity, all
such galaxies are identical up to scaling. This is no longer true in MOND as
it depends non-linearly on the typical acceleration, which is fully determined
by the central surface density. Thus, we set up a one-parameter family of
exponential disks and numerically determine their rotation curves (Figure 1)
and minimum $\sigma_{r}$ profiles (Figure 2). We also consider the stability
of the analogous galaxies in $\Lambda$CDM, which we assume have a DM halo that
causes their RC to match QUMOND predictions based on the baryons alone.

Throughout this work, we make use of the WKB approximation, namely that
perturbations are much smaller than the galactocentric radius. In Section 3.4,
we test the validity of this approximation. Our results show that it works
quite well beyond the central ${\approx 2r_{{}_{d}}}$ of MOND exponential
disks. It is expected to work particularly well for LSB disks in $\Lambda$CDM,
where the WKB approximation should be very accurate almost everywhere (bottom
panel of Figure 5).

The very central regions of galaxies may also be affected by non-axisymmetric
instabilities such as bars (e.g. Sellwood, 1981). Moreover, our results
pertain only to thin disks and thus cannot be applied to pressure-supported
galaxies or regions thereof. In the real Universe, even disk galaxies often
have a centrally concentrated bulge. We expect that our analysis applies
without modification to the regions outside the bulge, once the RC is
calculated appropriately and used to determine $\Omega_{r}$ in Equation 35.
The bulge would provide an additional source of stability in the regions
outside it by enhancing the RC666and thus $\Omega_{r}$ but not the disk
surface density. In this sense, a bulge would have a somewhat similar effect
to a DM halo, though an important difference is that bulges are often directly
observed while DM halos remain speculative (Liu et al., 2017).

Our analytic results provide a stability criterion but shed no light on how
exactly instability would develop and whether it could be saturated by non-
linear effects. Such questions need to be addressed with numerical
simulations. These indicate that rotationally supported self-gravitating disks
are unstable in Newtonian gravity (Hohl, 1971). They are generally thought to
be stabilised by a massive surrounding DM halo (Ostriker & Peebles, 1973). In
principle, this halo can grant an unlimited amount of stability to the disk
depending on their relative masses.

The situation is very different in MOND, where the modification to gravity can
only endow disks with a limited amount of extra stability. This remains true
for galaxies with arbitrarily low surface density, even though MOND has a very
large effect on the dynamics of such systems. Equation 47 gives a rough
understanding of why this is the case.

Our numerical calculations allow us to check the closeness of the action-based
QUMOND to the often-used algebraic MOND expression $\bm{g}=\nu\bm{g}_{{}_{N}}$
(Section 3.5). Our results indicate that beyond the central ${0.5\,r_{d}}$,
the radial force in the disk mid-plane differs by ${<3\%}$, even if the galaxy
has a very low surface density. However, this is only true if $\nu$ in the ALM
expression (Equation 55) is based on both the radial and vertical components
of $\bm{g}_{{}_{N}}$, with the latter assumed to be $2\pi G\Sigma$ as
appropriate for a point just outside the disk. If this component of
$\bm{g}_{{}_{N}}$ is neglected, then similarly good agreement between the ALM
and QUMOND is only attained beyond $\approx 2.5\,r_{d}$ (Figure 7).
Nonetheless, this form of the ALM (Equation 54) is correct in some modified
inertia interpretations of MOND (Milgrom, 1994, 2011). This may provide a way
to distinguish whether MOND is best understood as a modification of gravity or
of inertia. A recent analysis favoured the former (Frandsen & Petersen, 2018).

Our analytic disk stability condition for QUMOND (Equation 35) should prove
useful when setting up stable disk galaxies in the efficient $N$-body codes
Phantom of RAMSES (Lüghausen et al., 2015) and RAyMOND (Candlish et al., 2015)
that implement this theory. We are currently attempting to use the former to
simulate a past flyby interaction between the Milky Way and Andromeda
galaxies, thus extending our work in Banik et al. (2018) by performing
$N$-body simulations similar to those of Bílek et al. (2018). We hope to
clarify if this flyby scenario can form structures similar to those observed
in the Local Group.

##  6 Acknowledgements

IB is supported by an Alexander von Humboldt research fellowship. The authors
wish to thank the referee Agris Kalnajs for helping to clarify our assumptions
and their range of validity. The algorithms were set up using
matlab${}^{\text{\textregistered}}$.

## References

  * Aguerri et al. (2015) Aguerri J. A. L., et al., 2015, A&A, 576, A102
  * Algorry et al. (2017) Algorry D. G., et al., 2017, MNRAS, 469, 1054
  * Angus et al. (2012) Angus G. W., van der Heyden K. J., Famaey B., Gentile G., McGaugh S. S., de Blok W. J. G., 2012, MNRAS, 421, 2598
  * Angus et al. (2015) Angus G. W., Gentile G., Swaters R., Famaey B., Diaferio A., McGaugh S. S., Heyden K. J. v. d., 2015, MNRAS, 451, 3551
  * Apostol (1983) Apostol T. M., 1983, The Mathematical Intelligencer, 5, 59
  * Banik & Zhao (2018a) Banik I., Zhao H., 2018a, SciFed Journal of Astrophysics, 1, 1000008
  * Banik & Zhao (2018b) Banik I., Zhao H., 2018b, MNRAS, 473, 419
  * Banik & Zhao (2018c) Banik I., Zhao H., 2018c, MNRAS, 478, 3878
  * Banik & Zhao (2018d) Banik I., Zhao H., 2018d, MNRAS, 480, 2660
  * Banik et al. (2018) Banik I., O’Ryan D., Zhao H., 2018, MNRAS, 477, 4768
  * Begeman et al. (1991) Begeman K. G., Broeils A. H., Sanders R. H., 1991, MNRAS, 249, 523
  * Bekenstein & Milgrom (1984) Bekenstein J., Milgrom M., 1984, ApJ, 286, 7
  * Bell & de Jong (2001) Bell E. F., de Jong R. S., 2001, ApJ, 550, 212
  * Bershady et al. (2010) Bershady M. A., Verheijen M. A. W., Swaters R. A., Andersen D. R., Westfall K. B., Martinsson T., 2010, ApJ, 716, 198
  * Bílek et al. (2018) Bílek M., Thies I., Kroupa P., Famaey B., 2018, A&A, 614, A59
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Brada & Milgrom (1995) Brada R., Milgrom M., 1995, MNRAS, 276, 453
  * Brada & Milgrom (1999) Brada R., Milgrom M., 1999, ApJ, 519, 590
  * Caldwell et al. (2017) Caldwell N., et al., 2017, ApJ, 839, 20
  * Candlish (2016) Candlish G. N., 2016, MNRAS, 460, 2571
  * Candlish et al. (2015) Candlish G. N., Smith R., Fellhauer M., 2015, MNRAS, 446, 1060
  * Combes (2016) Combes F., 2016, in Laurikainen E., Peletier R., Gadotti D., eds, Astrophysics and Space Science Library Vol. 418, Galactic Bulges. p. 413 (arXiv:1501.03603), doi:10.1007/978-3-319-19378-6˙15
  * Corsini (2011) Corsini E. M., 2011, Memorie della Societa Astronomica Italiana Supplementi, 18, 23
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * D’Onghia (2015) D’Onghia E., 2015, ApJL, 808, L8
  * Desmond (2017a) Desmond H., 2017a, MNRAS, 464, 4160
  * Desmond (2017b) Desmond H., 2017b, MNRAS, 472, L35
  * Disney (1976) Disney M. J., 1976, Nature, 263, 573
  * Dubinski et al. (2008) Dubinski J., Gauthier J.-R., Widrow L., Nickerson S., 2008, in Funes J. G., Corsini E. M., eds, Astronomical Society of the Pacific Conference Series Vol. 396, Formation and Evolution of Galaxy Disks. p. 321 (arXiv:0802.3997)
  * Fall & Lynden-Bell (1981) Fall S. M., Lynden-Bell D., 1981, The structure and evolution of normal galaxies. NATO Advanced Study Institute held at the Institute of Astronomy and Clare College, Cambridge, UK, 3 - 15 August 1980.. Cambridge University Press
  * Famaey & Binney (2005) Famaey B., Binney J., 2005, MNRAS, 363, 603
  * Famaey & McGaugh (2012) Famaey B., McGaugh S. S., 2012, Living Reviews in Relativity, 15, 10
  * Famaey et al. (2018) Famaey B., McGaugh S., Milgrom M., 2018, MNRAS, 480, 473
  * Fathi (2010) Fathi K., 2010, ApJL, 722, L120
  * Frandsen & Petersen (2018) Frandsen M. T., Petersen J., 2018, preprint, Arxiv (arXiv:1805.10706)
  * Freeman (1970) Freeman K. C., 1970, ApJ, 160, 811
  * Fuchs (2003) Fuchs B., 2003, ApSS, 284, 719
  * Gradshteyn et al. (2015) Gradshteyn I. S., Ryzhik I. M., Zwillinger D., Moll V., 2015, Table of Integrals, Series, and Products, eighth edition edn. Academic Press, Boston
  * Hohl (1971) Hohl F., 1971, ApJ, 168, 343
  * Hu & Sijacki (2018) Hu S., Sijacki D., 2018, MNRAS, 478, 1576
  * Jones-Smith et al. (2018) Jones-Smith K., Abraham R., Kell L., Mathur H., 2018, New Journal of Physics, 20, 063042
  * Kahn & Woltjer (1959) Kahn F. D., Woltjer L., 1959, ApJ, 130, 705
  * Kroupa et al. (2018) Kroupa P., et al., 2018, Nature, 561, E4
  * Kuzio de Naray et al. (2012) Kuzio de Naray R., Arsenault C. A., Spekkens K., Sellwood J. A., McDonald M., Simon J. D., Teuben P., 2012, MNRAS, 427, 2523
  * Lelli et al. (2017) Lelli F., McGaugh S. S., Schombert J. M., Pawlowski M. S., 2017, ApJ, 836, 152
  * Li et al. (2018) Li P., Lelli F., McGaugh S., Schombert J., 2018, A&A, 615, A3
  * Lin & Shu (1964) Lin C. C., Shu F. H., 1964, ApJ, 140, 646
  * Liu et al. (2017) Liu J., Chen X., Ji X., 2017, Nature Physics, 13, 212
  * Londrillo & Nipoti (2009) Londrillo P., Nipoti C., 2009, Memorie della Societa Astronomica Italiana Supplementi, 13, 89
  * Lüghausen et al. (2015) Lüghausen F., Famaey B., Kroupa P., 2015, Canadian Journal of Physics, 93, 232
  * McGaugh (1996) McGaugh S. S., 1996, MNRAS, 280, 337
  * McGaugh (2011) McGaugh S. S., 2011, Physical Review Letters, 106, 121303
  * McGaugh (2016) McGaugh S. S., 2016, ApJL, 832, L8
  * McGaugh & Milgrom (2013) McGaugh S., Milgrom M., 2013, ApJ, 775, 139
  * McGaugh & Wolf (2010) McGaugh S. S., Wolf J., 2010, ApJ, 722, 248
  * McGaugh & de Blok (1998) McGaugh S. S., de Blok W. J. G., 1998, ApJ, 499, 66
  * McGaugh et al. (1995) McGaugh S. S., Schombert J. M., Bothun G. D., 1995, AJ, 109, 2019
  * Milgrom (1983a) Milgrom M., 1983a, ApJ, 270, 365
  * Milgrom (1983b) Milgrom M., 1983b, ApJ, 270, 371
  * Milgrom (1986) Milgrom M., 1986, ApJ, 302, 617
  * Milgrom (1988) Milgrom M., 1988, ApJ, 333, 689
  * Milgrom (1989) Milgrom M., 1989, ApJ, 338, 121
  * Milgrom (1994) Milgrom M., 1994, Annals of Physics, 229, 384
  * Milgrom (1999) Milgrom M., 1999, Phys. Lett. A, 253, 273
  * Milgrom (2009a) Milgrom M., 2009a, Physical Review D, 80, 123536
  * Milgrom (2009b) Milgrom M., 2009b, ApJ, 698, 1630
  * Milgrom (2010) Milgrom M., 2010, MNRAS, 403, 886
  * Milgrom (2011) Milgrom M., 2011, Acta Physica Polonica B, 42, 2175
  * Milgrom (2012) Milgrom M., 2012, Physical Review Letters, 109, 131101
  * Milgrom (2014) Milgrom M., 2014, Scholarpedia, 9
  * Milgrom (2016) Milgrom M., 2016, Physical Review Letters, 117, 141101
  * Ostriker & Peebles (1973) Ostriker J. P., Peebles P. J. E., 1973, ApJ, 186, 467
  * Pazy (2013) Pazy E., 2013, Phys. Rev. D, 87, 084063
  * Pereira et al. (2016) Pereira J. P., Overduin J. M., Poyneer A. J., 2016, Physical Review Letters, 117, 071103
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, ApJ, 517, 565
  * Peters (1981) Peters P. C., 1981, American Journal of Physics, 49, 564
  * Peters & Kuzio de Naray (2018) Peters W., Kuzio de Naray R., 2018, MNRAS, 476, 2938
  * Renaud et al. (2016) Renaud F., Famaey B., Kroupa P., 2016, MNRAS, 463, 3637
  * Riess et al. (1998) Riess A. G., et al., 1998, AJ, 116, 1009
  * Roberts & Whitehurst (1975) Roberts M. S., Whitehurst R. N., 1975, ApJ, 201, 327
  * Rogstad & Shostak (1972) Rogstad D. H., Shostak G. S., 1972, ApJ, 176, 315
  * Rubin & Ford (1970) Rubin V. C., Ford Jr. W. K., 1970, ApJ, 159, 379
  * Saburova (2011) Saburova A. S., 2011, Astronomy Reports, 55, 409
  * Safronov (1960) Safronov V. S., 1960, Annales d’Astrophysique, 23, 979
  * Salucci & Turini (2017) Salucci P., Turini N., 2017, preprint, Arxiv (arXiv:1707.01059)
  * Sánchez-Salcedo et al. (2016) Sánchez-Salcedo F. J., Martínez-Gómez E., Aguirre-Torres V. M., Hernández-Toledo H. M., 2016, MNRAS, 462, 3918
  * Sanders (1996) Sanders R. H., 1996, ApJ, 473, 117
  * Sanders & Noordermeer (2007) Sanders R. H., Noordermeer E., 2007, MNRAS, 379, 702
  * Sanders & Verheijen (1998) Sanders R. H., Verheijen M. A. W., 1998, ApJ, 503, 97
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Sellwood (1981) Sellwood J. A., 1981, A&A, 99, 362
  * Serra et al. (2012) Serra P., et al., 2012, MNRAS, 422, 1835
  * Serra et al. (2016) Serra P., Oosterloo T., Cappellari M., den Heijer M., Józsa G. I. G., 2016, MNRAS, 460, 1382
  * Smolin (2017) Smolin L., 2017, Physical Review D, 96, 083523
  * Spitzer (1942) Spitzer Jr. L., 1942, ApJ, 95, 329
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Thies et al. (2016) Thies I., Kroupa P., Famaey B., 2016, preprint, Arxiv (arXiv:1606.04942)
  * Thomas et al. (2018a) Thomas G. F., Famaey B., Ibata R., Renaud F., Martin N. F., Kroupa P., 2018a, A&A, 609, A44
  * Thomas et al. (2018b) Thomas G. F., Famaey B., Ibata R., Renaud F., Martin N. F., Kroupa P., 2018b, A&A, 609, A44
  * Tiret & Combes (2007) Tiret O., Combes F., 2007, A&A, 464, 517
  * Tiret & Combes (2008) Tiret O., Combes F., 2008, A&A, 483, 719
  * Toomre (1964) Toomre A., 1964, ApJ, 139, 1217
  * Tully et al. (2013) Tully R. B., et al., 2013, AJ, 146, 86
  * Verlinde (2016) Verlinde E. P., 2016, SciPost Physics, 2, 16
  * Zhao et al. (2013) Zhao H., Famaey B., Lüghausen F., Kroupa P., 2013, A&A, 557, L3
  * de Blok & McGaugh (1997) de Blok W. J. G., McGaugh S. S., 1997, MNRAS, 290, 533
  * den Heijer et al. (2015) den Heijer M., et al., 2015, A&A, 581, A98
  * van der Kruit (1987) van der Kruit P. C., 1987, A&A, 173, 59

##  Appendix A Root mean square thickness of a sech2 profile

To allow easier comparison between different density profiles, a standard of
comparison is required. This is most easily done using the widely used root
mean square (rms) measure of thickness. The relation between rms and
$\operatorname{sech}^{2}$ scale heights can be found by solving

$\displaystyle I~{}$ $\displaystyle\equiv$
$\displaystyle~{}\int_{0}^{\infty}x^{2}\operatorname{sech}^{2}x\,dx$ (57)
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{4x^{2}\mathrm{e}^{2x}}{\left(\mathrm{e}^{x}+\mathrm{e}^{-x}\right)^{2}\mathrm{e}^{2x}}\,dx$
(58) $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{\infty}\frac{x^{2}\mathrm{e}^{x}}{\left(\mathrm{e}^{x}+1\right)^{2}}\,dx\,.$
(59)

On the last line, we substitute for ${2x}$. Some further progress can be made
using integration by parts.

$\displaystyle
I~{}=~{}\frac{1}{2}\left(\left.-\,\frac{x^{2}}{\mathrm{e}^{x}+1}\right|_{0}^{\infty}\,+\,2\int_{0}^{\infty}\frac{x}{\mathrm{e}^{x}+1}\,dx\right)\,.$
(60)

The first term does not contribute to $I$ because it is 0 both when $x=0$ and
as ${x\to\infty}$. The second term can be expressed as an infinite sum.

$\displaystyle I~{}$ $\displaystyle=$
$\displaystyle~{}\int_{0}^{\infty}\frac{x\mathrm{e}^{-x}}{1+\mathrm{e}^{-x}}\,dx$
(61) $\displaystyle=$
$\displaystyle~{}-\sum_{n=1}^{\infty}\int_{0}^{\infty}x\mathrm{e}^{-nx}\left(-1\right)^{n}\,dx\,.$
(62)

Each term in this sum is of the form

$\displaystyle\int_{0}^{\infty}x\mathrm{e}^{-nx}\,dx~{}=~{}\frac{1}{n^{2}}\,.$
(63)

This allows us to write $I$ in terms of the Riemann zeta function of 2, which
can be shown to have a value of $\frac{\mathrm{\pi}^{2}}{6}$ using only basic
integrals and infinite sums similar to those just used (Apostol, 1983). In our
case, the contributions from even values of $n$ must be subtracted rather than
added, reducing the final result by a factor of
$\left(1-\frac{2}{2^{2}}\right)$. Therefore,

$\displaystyle I~{}=~{}\frac{\mathrm{\pi}^{2}}{12}\,.$ (64)

This reproduces the result stated in section 3.536 of Gradshteyn et al.
(2015). Given that $\int_{0}^{\infty}\operatorname{sech}^{2}x\,dx=1$, a
density profile with
$\rho\propto\operatorname{sech}^{2}\left(\frac{z}{h}\right)$ must have a rms
thickness of $\frac{\mathrm{\pi}h}{\sqrt{12}}$, which is only 9.3% less than
$h$.