# Weak-Coupling Theory of Neutron Scattering as a Probe of Altermagnetism

Thomas A. Maier  Computational Sciences and Engineering Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37831, USA  Satoshi Okamoto
Materials Science and Technology Division, Oak Ridge National Laboratory, Oak
Ridge, Tennessee 37831, USA

###### Abstract

Inelastic neutron scattering provides a powerful probe of the magnetic
excitations of quantum magnets. Altermagnets have recently emerged as a new
class of magnets with vanishing net magnetization characteristic of
antiferromagnets and with a spin-split electronic structure typical of
ferromagnets. Here we introduce a minimal Hubbard model with two-sublattice
orthorhombic anisotropy as a framework to study altermagnetism. Using
unrestricted Hartree-Fock calculations, we find an altermagnetic state for
this model that evolves from a metallic state to an insulating state with
increasing Hubbard-$U$ Coulomb repulsion. We then examine the inelastic
neutron scattering response in these states using random-phase approximation
calculations of the dynamic spin susceptibility $\chi^{\prime\prime}({\bf
q},\omega)$. We find that the magnetic excitation spectrum depends on its
chirality for ${\bf q}$ along certain directions in reciprocal space, an
observation that may be used in inelastic neutron scattering experiments as a
probe of altermagnetism.

##  I Introduction

Altermagnetism has emerged as a new class of magnetism [1, 2, 3, 4]. As
conventional antiferromagnets, altermagnets have an alternating spin
polarization in real space. In contrast to antiferromagnets, however, the
magnetic ordering in altermagnets breaks the $\cal PT$ (parity times time-
reversal) symmetry. As a consequence, the spin splitting also appears in
reciprocal space in the electronic band structure. A large number of candidate
materials have already been identified to host altermagnetism, including RuO2
[5, 6, 7], FeSb2 [8], and other materials [1, 2]. When rotational ${\cal
C}_{4}$ symmetry is preserved, such as in RuO2, the band-dependent spin
polarization has $d_{x^{2}-y^{2}}$-symmetry, and the analogy with
unconventional $d$-wave superconductivity as observed in the cuprates has been
pointed out [2]. Moreover, based on their unique electronic structure, several
phenomena have been proposed that are especially useful for spintronics
applications with itinerant electrons [1, 2, 9, 10].

While the unique electronic structure of altermagnets has received most
attention, especially in the context of spintronics applications, the momentum
dependent spin splitting and spin dynamics have just now started to become the
focus of research. For example, Ghorashi et al. proposed a heterostructure
consisting of an altermagnet and a superconductor as a platform to create
Majorana zero modes with net zero magnetization [11]. Steward et al. used a
phenomenological model for the spin dynamics to propose a way to assess the
spin excitations in altermagnets directly from the phonon spectrum utilizing
the coupling between spin excitations and lattices degree of freedom [12].
Moreover, Šmejkal et al. theoretically examined the spin excitations in RuO2
[13] based on strong coupling linear spin-wave theory, and found an
alternating chirality splitting in the magnon bands different from the
chirality-degenerate magnon bands in conventional antiferromagnets. Many
candidate altermagnets, however, are weakly correlated itinerant systems, such
as RuO2 for example, and especially the materials that are more favorable for
practical applications in spintronics. In that case, a more natural starting
point is a weak coupling theory based on an itinerant picture.

Here we use weak coupling theory to examine the spin dynamics in altermagnets
and explore how inelastic neutron scattering in the magnetically ordered state
can be used to distinguish between a conventional antiferromagntic and an
altermagnetic state. To this end, we introduce a minimal microscopic Hubbard
model with an on-site Coulomb $U$ repulsion for an altermagnet, and a
corresponding model for an (anisotropic) antiferromagnet. We first study the
phase diagram of these models for varying $U$ using Hartree Fock (HF) theory.
Next, we proceed to explore the structure of the inelastic neutron scattering
in the magnetic state for both altermagnetic and antiferromagnetic models
using a weak-coupling random phase approximation.

##  II Model

To model the altermagnetic state, we start with a two-dimensional (2D) Hubbard
model on a square lattice with the total number lattice sites $L$. Its
Hamiltonian

$\displaystyle
H=-\sum_{ij,\sigma}(t_{ij}+\mu\delta_{ij})c^{\dagger}_{i\sigma}c^{\phantom{\dagger}}_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}\,,$
(1)

is formulated in terms of the usual electron creation
($c^{\dagger}_{i\sigma}$) and annihilation ($c_{i\sigma}$) operators for sites
$i$ and spin $\sigma$, their corresponding number operators
$n_{i\sigma}=c^{\dagger}_{i\sigma}c^{\phantom{\dagger}}_{i\sigma}$, the
single-particle hopping amplitude $t_{ij}$ between sites $i$ and $j$, and the
on-site Coulomb repulsion $U$. $\mu$ is the chemical potential. As shown in
Fig. 1, panels (a) and (b), the 2D square lattice is divided into A- and
B-sublattices to allow for different intra-sublattice (A-A and B-B) hopping
terms. We include three terms for the hopping amplitude $t_{ij}$: A nearest
neighbor (inter-sublattice A-B) term $t_{0}$ and third-nearest-neighbor
(intra-sublattice A-A and B-B) terms $t_{1}$ and $t_{2}$, which break the
sublattice ${\cal C}_{4}$ symmetry. The minimal model we propose to describe
the altermagnetic state is illustrated in Fig. 1 (a). On sublattice A, it has
intra-sublattice hopping $t_{1}$ along $x$ and $t_{2}$ along $y$. On the
B-sublattice, $t_{1}$ and $t_{2}$ are switched. Fourier-transforming the first
(non-interacting) term of the Hamiltonian in Eq. (1) to reciprocal space gives

$H_{0}=\sum_{{\bf k},\sigma}\psi^{\dagger}_{{\bf k}\sigma}{\bf h}({\bf
k})\psi^{\phantom{\dagger}}_{{\bf k}\sigma}\,.$ (2)

Here, ${\bf k}$ is the crystal momenta of a sublattice, $\psi^{\dagger}_{{\bf
k}\sigma}=(c^{\dagger}_{A{\bf k}\sigma},c^{\dagger}_{B{\bf k}\sigma})$ and

$\displaystyle{\bf h}({\bf k})$ $\displaystyle=$
$\displaystyle\frac{1}{2}[\epsilon_{AA}({\bf k})+\epsilon_{BB}({\bf
k})]\tau_{0}+\epsilon_{AB}({\bf k})\tau_{1}$ (3) $\displaystyle+$
$\displaystyle\frac{1}{2}[\epsilon_{AA}({\bf k})-\epsilon_{BB}({\bf
k})]\tau_{3}-\mu\tau_{0}$

with

$\displaystyle\epsilon_{AA}({\bf k})$ $\displaystyle=$
$\displaystyle-2t_{1}\cos 2k_{x}-2t_{2}\cos 2k_{y},$ (4)
$\displaystyle\epsilon_{BB}({\bf k})$ $\displaystyle=$
$\displaystyle-2t_{2}\cos 2k_{x}-2t_{1}\cos 2k_{y},$ (5)
$\displaystyle\epsilon_{AB}({\bf k})$ $\displaystyle=$
$\displaystyle-2t_{0}(\cos k_{x}+\cos k_{y})\,.$ (6)

and $\tau_{i}$ the Pauli matrices in sublattice space. As we will see, the
bandstructure of this model leads to an altermagnetic state that breaks $\cal
PT$ symmetry.

Figure 1: Illustration of the models considered in the calculations: (a)
altermagnet when sublattices A and B have opposite spin polarizations and (b)
antiferromagnet which has an orthorhombic anisotropy.

For comparison, we will also consider an analogous model with translationally
invariant hopping amplitudes that supports a conventional antiferromagnetic
state. As illustrated in Fig. 1 (b), this model has the same $t_{1}$ and
$t_{2}$ hopping terms on both A- and B-sublattices so that $\epsilon_{BB}({\bf
k})=\epsilon_{AA}({\bf k})$. In this case, we will see that its magnetic
bandstructure is two-fold (spin) degenerate and therefore does not violate the
$\cal PT$ symmetry.

In the following we take $t_{0}=1$ as the unit of energy and set
$t_{1}=0.4t_{0}$ and $t_{2}=0.2t_{0}$. The chemical potential $\mu$ is
adjusted so that the model is half-filled with an average site filling
$\langle n\rangle$ = 1.

##  III Hartree-Fock phase diagram

We start by analyzing the models shown in Fig. 1 with Hamiltonian in Eq. (1)
using an unrestricted Hartree-Fock approximation. This approximation allows
one to investigate possible symmetry-broken phases of interacting models and
excitations from such phases.

To proceed, we introduce a mean-field decoupling in the interaction term of
Eq. (1) as

$\displaystyle Un_{i\uparrow}n_{i\downarrow}$ $\displaystyle\approx$
$\displaystyle U\bigl{\\{}n_{i\uparrow}\langle n_{i\downarrow}\rangle+\langle
n_{i\uparrow}\rangle n_{i\downarrow}-\langle n_{i\uparrow}\rangle\langle
n_{i\downarrow}\rangle\bigr{\\}}$ (7) $\displaystyle=$
$\displaystyle-\sum_{\sigma}\frac{1}{2}U(M_{i}\sigma-N_{i})n_{i\sigma},$

where $M_{i}\equiv\langle n_{i\uparrow}-n_{i\downarrow}\rangle$ and
$N_{i}\equiv\langle n_{i\uparrow}+n_{i\downarrow}\rangle$ are the spin
antisymmetric and spin symmetric combinations of $\langle n_{i\sigma}\rangle$,
respectively. In the second line of Eq. (7), the constant term $U\langle
n_{i\uparrow}\rangle\langle n_{i\downarrow}\rangle$ is neglected.

We explore a self-consistent solution, in which the electron density is
uniform while the spin polarization alternates between sublattices A and B,
such that $N_{A}=N_{B}=N$ and $M_{A}=-M_{B}=M$. This is achieved by
diagonalizing the $4\times 4$ mean-field Hamiltonian

$H_{0}=\sum_{{\bf k}}\psi^{\dagger}_{{\bf k}}[{\bf h}({\bf
k})\otimes\sigma_{0}-h\tau_{3}\otimes\sigma_{3}]\psi^{\phantom{\dagger}}_{{\bf
k}},$ (8)

and self-consistently computing the mean field order parameters $\langle
n_{i\sigma}\rangle$. Here, $\psi^{\dagger}_{{\bf k}}=(c^{\dagger}_{A{\bf
k}\uparrow},c^{\dagger}_{B{\bf k}\uparrow},c^{\dagger}_{A{\bf
k}\downarrow},c^{\dagger}_{B{\bf k}\downarrow})$, and and $\sigma_{i}$ are the
Pauli matrices acting on the spin space. $h=\frac{1}{2}UM$ is the staggered
Zeeman field, and $\frac{1}{2}UN$ is absorbed into the chemical potential
$\mu$. The self-consistency condition is given by

$\langle n_{\ell}\rangle=\frac{1}{L}\sum_{{\bf k},\nu}a_{\nu}^{\ell*}({\bf
k})a_{\nu}^{\ell}({\bf k})f(\xi_{\nu{\bf k}}-\mu),$ (9)

where $\ell=(l,\sigma)$ labels sublattice, $l=A,B$, and spin,
$\sigma=\uparrow,\downarrow$, $\xi_{\nu{\bf k}}$ are the $\bf k$-dependent
energy eigenvalues of band $\nu$ with the corresponding eigenvectors given by
$a^{\ell}_{\nu}({\bf k})=\langle\mu{\bf k}|\ell{\bf k}\rangle$, and
$f(\xi_{\nu{\bf k}}-\mu)$ is the Fermi distribution function.

Figure 2: Hartree-Fock results for altermagnetic (a-c) and antiferromagnetic
(d-f) models. (a, d) show the magnetization $M$ for the A-sublattice and the
corresponding Zeeman field $h$. (b, e) and (c, f) show the bandstructures in
the metallic states with $U=3.6t_{0}$ and the insulating states with
$U=4.0t_{0}$, respectively. Red and blue bands in (b, c) indicate the spin
$\uparrow$ and $\downarrow$ components, respectively, showing
$d_{x^{2}-y^{2}}$-type ${\bf k}$-dependent spin splitting, while gray bands in
(e, f) indicate spin degeneracy. The insets of (a, d) show the Fermi surfaces
in the metallic states with $U=3.6t_{0}$, with gray broken lines indicating
the reduced Brillouin zone of the staggered magnetic state. The Fermi surfaces
consist of hole pockets at $(\pm\pi/2,\pm\pi/2)$ and electron pockets at
$(\pm\pi,0)$ and $(0,\pm\pi)$.

Fig. 2 (a) shows the magnetization $M=\langle
n_{\uparrow}-n_{\downarrow}\rangle$ for the A-sublattice (the magnetization on
the B-sublattice is $-M$) as a function of the Coulomb $U$ repulsion. For the
altermagnetic model, as shown in panel (a), we find a magnetic state with
staggered spin polarization ($M$ on sublattice A, $-M$ on sublattice B) above
a threshold of $U\gtrsim 3.5t_{0}$. Right above this critical $U$, we find a
metallic state for a small range of $U$ as shown by the orange region, in
which the sublattice magnetization $M$ and the corresponding Zeeman mean field
$h$ rapidly rise with increasing $U$. For larger $U$, the bandstructure
becomes gapped and the increase of $M$ and $h$ slows down significantly in the
green shaded insulating region. The corresponding bandstructure in the
altermagnetic state is illustrated in panel (b) for the metallic system at
$U=3.6t_{0}$, and in panel (c) for an insulating system at $U=4t_{0}$. In both
cases, the bandstructure has four different bands resulting from the
sublattices and the spins, i.e. unlike in the antiferromagnetic state, the
spin degeneracy is lifted, and the splitting between spin up (red) and down
(blue) depends on the crystal momentum ${\bf k}$. As can be seen from the
color shading in the figure, the splitting has a $d_{x^{2}-y^{2}}$ structure
in momentum space. For the smaller $U$ metallic state in panel (b), the bands
cross the Fermi level resulting in four electron pockets at the corners
[$(\pi,0)$, $(0,\pi)$, etc.] of the reduced Brillouin zone, and four hole
pockets at the zone edges, as shown in the inset of Fig. 2 (a). One sees that
${\cal C}_{4}$ rotational symmetry of the bandstructure and Fermi surface
requires an additional spin-flip $\sigma\rightarrow\bar{\sigma}$ operation,
due to the sign reversal of the $h\tau_{3}\otimes\sigma_{3}$ term in Eq. (8)
under the ${\cal C}_{4}$ operation.

In the ordered state, ${\cal PT}$ symmetry is broken, since for a general
momentum ${\bf k}$, the degeneracy between ${\bf k}\uparrow$ and ${\bf
k}\downarrow$ states is lifted. This can be unterstood from the form of the
Hamiltonian in Eq. (8), since, under the ${\cal PT}$ operation, the
$\frac{1}{2}[\epsilon_{AA}({\bf k})-\epsilon_{BB}({\bf k})]\tau_{3}$ term in
${\bf h}({\bf k})$ transforms into its negative value. For the altermagnetic
model, this term is finite except along the diagonal $|k_{x}|=|k_{y}|$ and
zone boundary $|k_{x}\pm k_{y}|=\pi$ due to the difference in $t_{1}$ and
$t_{2}$ hoppings between the two sublattices.

The phase diagram shown in Fig. 2 (d) for the corresponding translationally
invariant model [Fig. 1 (b)] with $\epsilon_{BB}=\epsilon_{AA}$ is similar,
exhibiting an antiferromagnetic phase for $U\gtrsim 3.45t_{0}$ that evolves
from a metallic state in the orange region to an insulating state in the green
region. In contrast to the altermagnetic state, the corresponding
bandstructure of this model shown in panels (e) and (f) is spin degenerate and
lacks ${\cal C}_{4}$ symmetry since $t_{1}\neq t_{2}$. This can be seen from
the Fermi surface for the $U=3.6t_{0}$ metallic state shown in the inset of
panel (d). In this case, ${\cal PT}$ symmetry is preserved, because the
$\frac{1}{2}[\epsilon_{AA}({\bf k})-\epsilon_{BB}({\bf k})]\tau_{3}$ term in
${\bf h}({\bf k})$ vanishes for the antiferromagnetic model.

##  IV Inelastic neutron scattering

The dynamic spin susceptibility that determines the neutron scattering
intensity is given by

$\displaystyle\chi_{ij}({\bf q},i\omega_{m})$ $\displaystyle=$
$\displaystyle\sum_{l,l^{\prime}}\int^{1/T}_{0}d\tau
e^{i\omega_{m}\tau}\langle{\cal T}_{\tau}S_{i}^{l}({\bf
q},\tau)S_{j}^{l^{\prime}}({\bf-q},0)\rangle\,,$

for transferred momentum ${\bf q}$ and bosonic Matsubara frequencies
$\omega_{m}=2m\pi T$ with temperature $T$. We will consider both scattering
with positive chirality given by $\chi_{+-}$ and scattering with negative
chirality given by $\chi_{-+}$. These susceptibilities are calculated with the
spin-flip operators $S_{+}^{l}({\bf q})=\sum_{\bf
k}\psi^{\dagger}_{l\uparrow}({\bf
k+q})\psi^{\phantom{\dagger}}_{l\downarrow}({\bf k})$ and $S_{-}^{l}({\bf
q})=\sum_{\bf k}\psi^{\dagger}_{l\downarrow}({\bf
k+q})\psi^{\phantom{\dagger}}_{l\uparrow}({\bf k})$. In the multi-orbital RPA
approximation [14, 15], $\chi_{ij}$ is calculated from an orbital dependent
spin-susceptibility tensor

$\displaystyle\chi^{\ell_{1}\ell_{2}}_{\ell_{3}\ell_{4}}({\bf
q},i\omega_{m})=\left\\{\chi_{0}({\bf q},i\omega_{m})[1-{\cal U}\chi_{0}({\bf
q},i\omega_{m})]^{-1}\right\\}^{\ell_{1}\ell_{2}}_{\ell_{3}\ell_{4}}\,.$

Here, ${\cal U}$ is a 16$\times$16 matrix containing the on-site Coulomb
interaction $U$ with matrix elements ${\cal
U}^{\ell_{1}\ell_{2}}_{\ell_{3}\ell_{4}}=U\delta_{l_{1}l_{2}}\delta_{l_{1}l_{3}}\delta_{l_{1}l_{4}}[\delta_{\sigma_{1}\bar{\sigma_{2}}}\delta_{\sigma_{1}\sigma_{3}}\delta_{\sigma_{1}\bar{\sigma_{4}}}-\delta_{\sigma_{1}\sigma_{2}}\delta_{\sigma_{1}\bar{\sigma_{3}}}\delta_{\sigma_{1}\bar{\sigma_{4}}}]$,
and the matrix elements of the bare susceptibility $\chi_{0}$ are given by

$\displaystyle(\chi_{0})^{\ell_{1}\ell_{2}}_{\ell_{3}\ell_{4}}({\bf
q},i\omega_{m})=$ $\displaystyle-$ $\displaystyle\frac{T}{L}\sum_{{\bf
k},\varepsilon_{n}}G_{0}^{\ell_{3}\ell_{1}}({\bf
k+q},i\varepsilon_{n}+i\omega_{m})$ (10) $\displaystyle\times$ $\displaystyle
G_{0}^{\ell_{2}\ell_{4}}({\bf k},i\varepsilon_{n})\,,$

with $G_{0}^{\ell\ell^{\prime}}({\bf k},i\varepsilon_{n})$ the bare Green’s
function

$\displaystyle G_{0}^{\ell\ell^{\prime}}({\bf k},i\varepsilon_{n})$
$\displaystyle=$ $\displaystyle\sum_{\nu}\frac{a^{\ell}_{\nu}({\bf
k})a^{\ell^{\prime}*}_{\nu}({\bf k})}{i\varepsilon_{n}-\xi_{\nu{\bf k}}+\mu}$
(11)

and $\varepsilon_{n}=(2n+1)\pi T$ the fermionic Matsubara frequencies.
Carrying out the usual analytic continuation of Matsubara frequencies
$\omega_{m}$ to the real frequency $\omega$ axis, the physical spin
susceptibility $\chi_{ij}({\bf q},\omega)$ is then obtained as

$\displaystyle\chi_{ij}({\bf q},\omega)=$
$\displaystyle\sum_{l,l^{\prime}}\sum_{\alpha,\beta,\gamma,\delta}\sigma^{i}_{\alpha\beta}\sigma^{j}_{\gamma\delta}$
$\displaystyle\times$
$\displaystyle\chi^{l\alpha,l\beta}_{l^{\prime}\gamma,l^{\prime}\delta}({\bf
q},i\omega_{m}\rightarrow\omega+i\eta)\,.$

Here, the operators $\sigma^{+}=\frac{1}{2}(\sigma_{1}+i\sigma_{2})$ and
$\sigma^{-}=\frac{1}{2}(\sigma_{1}-i\sigma_{2})$ with the Pauli matrices
$\sigma_{i}$, and $\eta$ is a positive infinitesimal.

Figure 3: RPA spin susceptibility for the insulating state at $U=4t_{0}$: (a)
$\chi^{\prime\prime}_{+-}({\bf q},\omega)$ and (b)
$\chi^{\prime\prime}_{-+}({\bf q},\omega)$, and (c) their average, for the
altermagnetic model shown in Fig. 1(a), and (d) $\chi^{\prime\prime}_{+-}({\bf
q},\omega)=\chi^{\prime\prime}_{-+}({\bf q},\omega)$ for the orthorhombic
antiferromagnetic model shown in Fig. 1(b).

In Fig. 3 we show the imaginary part of the dynamic spin susceptibilities
$\chi^{\prime\prime}_{+-}({\bf q},\omega)$ [panel(a)] and
$\chi^{\prime\prime}_{-+}({\bf q},\omega)$ [panel (b)] and their average
[panel (c)] for the altermagnetic state calculated in the insulating region
for $U=4t_{0}$ along high-symmetry directions in the Brillouin zone. As a
comparison, panel (d) shows the spectrum obtained for the antiferromagnet with
orthorhombic anisotropy. For all cases, the spin-wave dispersion has minima at
${\bf q}=0$, $(\pi,\pi)$ and $(\pi/2,\pi/2)$ with the largest intensity at
$(\pi,\pi)$ because of the staggered magnetic ordering. A closer inspection of
the detailed dispersion near ${\bf q}=0$, however, reveals important
differences: Due to the orthorhombic anisotropy in the hopping amplitudes, the
dispersion along $q_{x}$ from ${\bf q}=0$ to $(\pi,0)$ is different from the
dispersion along $q_{y}$ from ${\bf q}=0$ to $(0,\pi)$. For the altermagnetic
model, this difference is reversed between the neutron scattering response
with positive chirality described by $\chi_{+-}$ shown in panel (a) and the
response with negative chirality shown by $\chi_{-+}$ in panel (b). This
difference is similar to that observed previously in calculations for an
effective Heisenberg model of RuO2 in Ref. [13]. As a result of the chirality
dependent spin-wave dispersion, the average of $\chi_{+-}$ and $\chi_{-+}$
shown in panel (c) displays two branches dispersing out of ${\bf q}=0$. For
the antiferromagnetic case shown in panel (d) obtained for the model with
orthorhombic anisotropy, a similar asymmetry between $q_{x}$ and $q_{y}$ is
seen. For this case, however, the spectrum does not depend on chirality, i.e.
$\chi_{+-}=\chi_{-+}$, and the total scattering intensity shows only a single
branch.

Figure 4: RPA spin susceptibility for the metallic state at $U=3.6t_{0}$: (a)
$\chi^{\prime\prime}_{+-}({\bf q},\omega)$ and (b)
$\chi^{\prime\prime}_{-+}({\bf q},\omega)$, and (c) their average, for the
altermagnetic model shown in Fig. 1(a), and (d) $\chi^{\prime\prime}_{+-}({\bf
q},\omega)=\chi^{\prime\prime}_{-+}({\bf q},\omega)$ for the orthorhombic
antiferromagnetic model shown in Fig. 1(b).

Fig. 4 shows the corresponding spectrum obtained for the metallic state at
$U=3.6t_{0}$. In this case, the overall intensity drops to about half of that
for the insulating state and the increased Landau damping due to the metallic
character is obvious. Similar to the insulating case (Fig. 3), however, a
clear asymmetry in the low energy spin excitations is visible near ${\bf
q}=0$. For the altermagnetic case [panels (a) and (b)], this asymmetry is
again reversed between different chiralities, and the average of $\chi_{+-}$
and $\chi_{-+}$ has two distinct branches for energies up to $\sim t_{0}$
[panel (c)].

In the RPA approach, the chirality dependence of the response can be traced to
an asymmetry in the propagation between spin $\uparrow$ and spin $\downarrow$
particles. As shown in the insets of Figs. 3(a) and 4(a), the neutron
scattering process with positive chirality involves a spin $S=1$ particle-hole
excitation where momentum ${\bf q}$ is transferred to a spin $\uparrow$
electron. Conversely, as shown in the insets of Figs. 3(b) and 4(b), a
scattering process with negative chirality involves a spin $S=-1$ particle-
hole excitation where momentum ${\bf q}$ is transferred to a spin $\downarrow$
electron. Since spin $\uparrow$ and spin $\downarrow$ electrons propagate
mainly on different sublattices, the difference in the magnetic excitations
described by $\chi_{+-}$ and $\chi_{-+}$ can be traced to the differences in
the electronic structure between the two sublattices.

## Conclusion

We have introduced a microscopic two-sublattice Hubbard model with
alternating, ${\cal C}_{4}$ symmetry breaking orthorhombic anisotropy in the
hopping parameters, which provides a minimal framework to study
altermagnetism. Using an unrestricted Hartree-Fock approximation, we have
found an altermagnetic state with staggered magnetization and spin-split
bandstructure, which evolves from a metallic state at small $U$ to an
insulating state at larger $U$. We have used this model and an RPA
approximation of the dynamic spin susceptibility to study the inelastic
neutron scattering response in the altermagnetic state, and compared this
response with that of a corresponding conventional antiferomagnetic model with
orthorhombic anisotropy, for both the metallic and insulating states. As one
would expect, we have found that the magnetic spectrum described by
$\chi^{\prime\prime}({\bf q},\omega)$ breaks ${\cal C}_{4}$ symmetry for both
these models, and in the case of the altermagnetic model, depends also on its
chirality, described by either $\chi^{\prime\prime}_{+-}$ or
$\chi^{\prime\prime}_{-+}$, for ${\bf q}$ along certain directions in the
Brillouin zone. For the conventional antiferromagnetic state, the magnetic
excitation spectrum is chirality-degenerate. The observation of chirality-
split magnetic excitations in inelastic neutron scattering experiments may
therefore be used as evidence for altermagnetism.

###### Acknowledgements.

We acknowledge useful discussions with Andrew Christianson, Hu Miao, and
Pyeongjae Park. This work was supported by the U.S. Department of Energy,
Office of Science, Basic Energy Sciences, Materials Sciences and Engineering
Division. This manuscript has been authored by UT-Battelle, LLC, under
Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The United
States Government retains and the publisher, by accepting the article for
publication, acknowledges that the United States Government retains a
nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for the
United States Government purposes. The Department of Energy will provide
public access to these results of federally sponsored research in accordance
with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-
access-plan).

## References

  * [1] L. Šmejkal, J. Sinova, and T. Jungwirth, Phys. Rev. X 12, 031042 (2022).
  * [2] L. Šmejkal, J. Sinova, and T. Jungwirth, Phys. Rev. X 12, 040501 (2022).
  * [3] S. Hayami, Y. Yanagi, and H. Kusunose, J. Phys. Soc. Jpn. 88, 123702 (2019).
  * [4] S. Hayami, Y. Yanagi, and H. Kusunose, Phys. Rev. B 102, 144441 (2020).
  * [5] T. Berlijn, P. C. Snijders, O. Delaire, H.-D. Zhou, T. A. Maier, H.-B. Cao, S.-X. Chi, M. Matsuda, Y. Wang, M. R. Koehler, P. R. C. Kent, and H. H. Weitering, Phys. Rev. Lett. 118, 077201 (2017).
  * [6] K.-H. Ahn, A. Hariki, K.-W. Lee, and J. Kuneš, Phys. Rev. B 99, 184432 (2019).
  * [7] L. Šmejkal, R. González-Hernandez, T. Jungwirth, and J. Sinova, Sci. Adv. 6, eaaz8809 (2020).
  * [8] I. Mazin, Phys. Rev. X 12 040002 (2022).
  * [9] A. Bose, N. J. Schreiber, R. Jain, D.-F. Shao, H. P. Nair, J. Sun, X. S. Zhang, D. A. Muller, E. Y. Tsymbal, D. G. Schlom, and D. C. Ralph, Nat. Electron 5, 267 (2022)
  * [10] Z. Feng, X. Zhou, L. Šmejkal, L. Wu, Z. Zhu, H. Guo, R. González-Hernández, X. Wang, H. Yan, P. Qin, X. Zhang, H. Wu, H. Chen, Z. Meng, L. Liu, Z. Xia, J. Sinova, T. Jungwirth, and Z. Liu, Nat. Electron 5, 735 (2022).
  * [11] S. A. A. Ghorashi, T. L. Hughes, J. Cano, preprint arXiv:2306.09413 (2023).
  * [12] C. R. W. Steward, R. M. Fernandes, J. Schmalian, preprint arxiv:2307.01855 (2023).
  * [13] L. Šmejkal, A. Marmodoro, K.-H. Ahn, R. González-Hernández, I. Turek, S. Mankovsky, H. Ebert, S. W. D’Souza, O. Šipr, J. Sinova, and T. Jungwirth, preprint arXiv:2211.13806 (2022).
  * [14] K. Kubo, Phys.Rev.B 75, 224509 (2007).
  * [15] S. Graser, T. A. Maier, P. J. Hirschfeld, and D. J. Scalapino, New J. Phys. 11, 25016 (2009).


# The underlying radial acceleration relation

Harry Desmond1
1Institute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama
Building, Portsmouth, PO1 3FX, UK

E-mail<EMAIL_ADDRESS>
###### Abstract

The radial acceleration relation (RAR) of late-type galaxies relates their
dynamical acceleration, $g_{\text{obs}}$, to that sourced by baryons alone,
$g_{\text{bar}}$, across their rotation curves. Literature fits to the RAR
have fixed the galaxy parameters on which the relation depends—distance,
inclination, luminosity and mass-to-light ratios—to their maximum a priori
values with an uncorrelated Gaussian contribution to the uncertainties on
$g_{\text{bar}}$ and $g_{\text{obs}}$. In reality these are free parameters of
the fit, contributing systematic rather than statistical error. Assuming a
range of possible functional forms for the relation with or without intrinsic
scatter (motivated by Modified Newtonian Dynamics with or without the external
field effect), I use Hamiltonian Monte Carlo to perform the full joint
inference of RAR and galaxy parameters for the SPARC dataset. This reveals the
intrinsic RAR underlying that observed. I find an acceleration scale
$a_{0}=(1.19\pm 0.04\,\text{(stat)}\pm 0.09\,\text{(sys)})\>\times\>10^{-10}$
m s-2, an intrinsic scatter $\sigma_{\text{int}}=(0.034\pm
0.01\,\text{(stat)}\pm 0.01\,\text{(sys)})$ dex (assuming the SPARC error
model is reliable) and weak evidence for the external field effect. I make
summary statistics of all my analyses publicly available for future SPARC
studies or applications of a calibrated RAR, for example redshift-independent
distance measurement.

###### keywords:

galaxies: formation – galaxies: fundamental parameters – galaxies: kinematics
and dynamics – galaxies: statistics – dark matter

††pubyear: 2023††pagerange: The underlying radial acceleration
relation–LABEL:LastPage

##  1 Introduction

Galaxies are observed to follow several tight and regular scaling relations
between their internal motions and morphology. The classical correlations are
the Tully–Fisher relation between rotation velocity and mass or luminosity in
late-type galaxies (e.g. Tully & Fisher 1977; McGaugh et al. 2000; Pizagno et
al. 2007) and the Fundamental Plane relating luminosity, size and velocity
dispersion (e.g. Djorgovski & Davis 1987; Dressler et al. 1987; Cappellari et
al. 2013, including its projection onto the mass–velocity plane, the
Faber–Jackson relation; Faber & Jackson 1976) in early types. A more recent
addition is the mass discrepancy–acceleration or radial acceleration relation
(Sanders, 1990; McGaugh, 2004; Lelli et al., 2017), relating the local total
acceleration, $g_{\text{obs}}$, to that sourced by baryons, $g_{\text{bar}}$,
across rotation curves. This subsumes the Tully–Fisher relation and provides
more detailed information about the gravitational potential in late-type
galaxies.

These relations provide the key evidence concerning the _mass discrepancy
problem_ in galaxies, namely that the motions of stars and gas imply far
higher dynamical than baryonic masses in a Newtonian analysis. In the
prevailing $\Lambda$ Cold Dark Matter ($\Lambda$CDM) cosmology this difference
is assumed to be made up by dark matter, leading to attempts to explain the
relations through the modelling of galaxy formation, the galaxy–halo
connection and halo mass distributions (e.g. Gnedin et al. 2007; Blanton et
al. 2008; Desmond & Wechsler 2015, 2017; Di Cintio & Lelli 2016; Ludlow et al.
2017; Navarro et al. 2017; Keller & Wadsley 2017; Desmond 2017; Tenneti et al.
2018; Paranjape & Sheth 2021). This explanation is made difficult by the fact
that galaxy formation in $\Lambda$CDM proceeds in a highly stochastic and
complicated manner, at odds with the simplicity of galaxy phenomenology at the
present day.

An alternative hypothesis is that Newtonian gravity breaks down at the galaxy
scale. Surprisingly, galaxy dynamics can be explained well by a model in which
$g_{\text{obs}}=g_{\text{bar}}$ for $g_{\text{bar}}\gg a_{0}$ and
$g_{\text{obs}}\propto g_{\text{bar}}^{1/2}$ as $g_{\text{bar}}\ll a_{0}$,
where $a_{0}\approx 10^{-10}$ m s-2 is a new fundamental constant. This
naturally leads to the observed simplicity in the aforementioned scaling
relations. Supplemented by an “interpolating function” that connects the
Newtonian and modified gravity regimes, this theory is known as Modified
Newtonian Dynamics (MOND; Milgrom 1983a, c, b) and has achieved some success
at explaining and even predicting galaxy behaviour (e.g. Famaey & McGaugh
2012; McGaugh & Milgrom 2013; Chae et al. 2020b). MOND has been incorporated
into a range of nonrelativistic and relativistic theories over the past four
decades, as reviewed most recently in Banik & Zhao (2022).

The RAR is MOND written in terms of observables for late-type galaxies. This
at once gives the relation central importance in the missing mass debate and
makes it the most sensitive probe of gravitational parameters within the MOND
paradigm. These include the acceleration constant $a_{0}$ marking the onset of
modified dynamics, the intrinsic scatter $\sigma_{\text{int}}$ and possibly a
parameter $e_{\text{N}}\equiv g_{\text{ext}}/a_{0}$ describing the influence
of mass surrounding the galaxy (external field effect, EFE; Milgrom 1983a),
where $g_{\text{ext}}$ is the strength of the gravitational field in which the
galaxy is embedded. $\sigma_{\text{int}}$ bears on the question of whether the
RAR manifests law-like gravitational behaviour as posited by MOND, and also
determines the precision with which the relation may be used to calibrate
galaxy properties such as distance (analogously to the Tully–Fisher relation).
$e_{\text{N}}$ addresses the key question within the MOND paradigm of the
extent to which—and manner in which—modified gravity or inertia violates the
strong equivalence principle. Previous fits have found $a_{0}\approx 1.2\times
10^{-10}$ m/s2, $e_{\text{N}}\approx 0.003$ and a small intrinsic scatter
$\sigma_{\text{int}}<0.1$ dex (Lelli et al., 2017; Li et al., 2018; Chae et
al., 2021, 2022).

$g_{\text{bar}}$ and $g_{\text{obs}}$ depend on a number of properties of the
galaxies, most importantly their distance $D$, inclination $i$, luminosity $L$
and mass-to-light ratios $\Upsilon$ of their various components. These are
nuisance parameters when determining the properties of the RAR, although of
course of interest in their own right. Past RAR studies have either fixed
these to their maximum a priori values given other measurements and then
propagated their uncertainties into $g_{\text{bar}}$ and $g_{\text{obs}}$ as
if they were random and uncorrelated (Lelli et al., 2017), or varied both the
nuisance _and_ RAR parameters galaxy-by-galaxy, effectively assuming a
different RAR for each galaxy (Li et al., 2018; Chae et al., 2020b, 2021,
2022). Assuming an underlying universal form for the RAR, a superior inference
constrains _global_ RAR parameters along with the local galaxy properties. The
main advantage of this is that it propagates the prior distributions of the
galaxy parameters as systematic rather than statistical uncertainties, thus
capturing the correlations across rotation curves that fluctuations in these
parameters induce. For example, a higher (lower) $\Upsilon$ than expected in a
particular galaxy causes a higher (lower) $g_{\text{bar}}$ across its rotation
curve, yet modelling it as a statistical uncertainty implicitly assumes that a
fluctuation in $\Upsilon$ could scatter $g_{\text{bar}}$ up at one point and
down at the next. The full inference also captures the degeneracies between
the RAR and galaxy parameters, which have a non-trivial impact on the relation
through the shape of the galaxy priors. This is the analysis I perform here.

Although conceptually simple, the full inference is technically challenging
because it implies a vastly higher-dimensional parameter space than the
simplified versions. The analysis of Lelli et al. (2017) has two parameters
($a_{0}$, $\sigma_{\text{int}}$), while that of Li et al. (2018) has four
($D$, $i$, $\Upsilon_{\text{disk}}$, $\Upsilon_{\text{bulge}}$) repeated
$N=147$ times for $N$ SPARC galaxies. (Chae et al. 2020b additionally sample
$\Upsilon_{\text{gas}}$ and $e_{\text{N}}$.) Li et al. and Chae et al. cannot
accommodate parameters that couple the galaxies, so fix $a_{0}=1.2\times
10^{-10}$ m s-2 a priori and can at best reconstruct $\sigma_{\text{int}}$
post-hoc from the distribution of residuals, thus neglecting its degeneracy
with the other variables. The full inference has up to $6N+n+2$ parameters
($a_{0}$, $\sigma_{\text{int}}$, $N\times e_{\text{N}}$, $N\times D$, $N\times
i$, $N\times L_{3.6}$, $N\times\Upsilon_{\text{disk}}$,
$N\times\Upsilon_{\text{gas}}$, $n\times\Upsilon_{\text{bulge}}$) where $n=31$
is the number of galaxies with bulges. 915 parameters is beyond many sampling
methods but routine for Hamiltonian Monte Carlo. This will enable a robust
determination of the RAR parameters for arbitrary priors and assumptions about
the underlying functional form. $g_{\text{bar}}$ and $g_{\text{obs}}$
transformed according to the best-fit galaxy parameter values (Fig. 1) reveals
the RAR that underlies the sampling distributions of those parameters.

The structure of this paper is as follows. In Sec. 2 I describe the SPARC data
and selection criteria I employ. Sec. 3 gives the methodology, including the
likelihood model, priors, treatment of the galaxy parameters and details of
the sampler. The results are presented in Sec. 4. Sec. 5 discusses the broader
ramifications of the study, remaining systematic uncertainties and useful
further work, while Sec. 6 concludes. Throughout, $\log$ has base 10 and
accelerations are given in $10^{-10}$ m s-2 unless otherwise stated.

##  2 Observational Data

I analyse the SPARC sample (Lelli et al.,
2016),111http://astroweb.cwru.edu/SPARC/ comprising 175 rotation curves from
the literature with Spitzer photometry at 3.6$\mu$m. I apply the quality cuts
recommended by Lelli et al. (2017), removing galaxies with quality flag 3
(indicating large asymmetries, non-circular motions and/or offsets between
stellar and Hi distributions) or maximum a priori $i<30\deg$, and points for
which the quoted fractional uncertainty on the observed rotation velocity is
greater than 10 per cent. This leaves $2696$ points from $147$ galaxies, of
which all have mass in a stellar disk but only $31$ have mass in a central
bulge.

Distances are determined by a variety of methods with a corresponding range of
uncertainties (Lelli et al., 2016), while the inclinations are estimated from
tilted-ring fits to the velocity fields. I use these as Gaussian priors in the
inference. The total luminosity at 3.6$\>\mu$m, $L_{3.6}$, is well-measured
but its uncertainty is quoted so I include it as a Gaussian prior for
completeness and to eliminate statistical uncertainty in the independent
($g_{\text{bar}}$) direction which complicates the inference (see Sec. 3.2).
The disk and bulge mass-to-light ratios, $\Upsilon_{\text{disk}}$ and
$\Upsilon_{\text{bulge}}$, are believed to be $\sim 0.5$ and $\sim 0.7$
respectively, with a $\sim 25$ per cent uncertainty (Meidt et al., 2014;
McGaugh & Schombert, 2014; Lelli et al., 2016). I use these as lognormal
priors, which are marginally favoured over Gaussian given the way the
parameters are determined (S. McGaugh and F. Lelli, priv. comm.).

With Hi mass measured, a correction factor must be applied to calculate the
total gas mass and hence the gas contribution to $g_{\text{bar}}$. The
fiducial SPARC analysis uses a conversion factor of 1.33 (accounting for
primordial helium), but a more accurate determination includes a scaling of
the hydrogen fraction with the stellar mass $M_{*}$ of the galaxy (McGaugh et
al., 2020):

$M_{\text{gas}}=X^{-1}\>M_{\text{HI}}$ (1)

where

$X(M_{*})=0.75-38.2\>\left(M_{*}/(1.5\times 10^{24}M_{\odot})\right)^{0.22}.$
(2)

As the Hi mass has already been scaled by 1.33 in SPARC, I define

$\bar{\Upsilon}_{\text{gas}}(M_{*})=1/(1.33\>X(M_{*}))$ (3)

where overbar denotes maximum a priori value. ($M_{*}$ must be determined
after sampling $L_{3.6}$, $\Upsilon_{\text{disk}}$ and
$\Upsilon_{\text{bulge}}$.) This scales $\Upsilon_{\text{gas}}$ relative to
the value assumed in SPARC when calculating $V_{\text{gas}}$, as
$\Upsilon_{\text{disk}}$ and $\Upsilon_{\text{bulge}}$ do for
$V_{\text{disk}}$ and $V_{\text{bul}}$. The results are not significantly
altered compared to $\bar{\Upsilon}_{\text{gas}}=1$. $\Upsilon_{\text{gas}}$
is given a lognormal prior with 10 per cent width (Lelli et al., 2016).222To
convert between normal and lognormal distributions I use the full equations
relating their means and standard deviation
$\displaystyle\mu=\exp(\tilde{\mu}+\tilde{\sigma}^{2}/2),$
$\displaystyle\sigma^{2}=(\exp(\tilde{\sigma}^{2})-1)\exp(2\tilde{\mu}+\tilde{\sigma}^{2}),$
with inverse  $\displaystyle\tilde{\mu}=\ln(\mu)-\ln(1+\sigma^{2}/\mu^{2})/2,$
(5) $\displaystyle\tilde{\sigma}^{2}=\ln(1+\sigma^{2}/\mu^{2}),$ where a tilde
indicates the lognormal. The uncertainties are sufficiently small in most
cases for this not to differ appreciably from the more common first-order
approximation.

##  3 Method

Figure 1: The underlying RAR of the SPARC sample is obtained by transforming
$g_{\text{bar}}$ and $g_{\text{obs}}$ according to the best-fit galaxy
parameters, in this case those at the median of the posterior for the
inference with intrinsic scatter but without EFE. The Simple IF fit with
$a_{0}$ in its $2\sigma$ allowed range is overplotted in red, and the median
errorbar size, deriving solely from the statistical uncertainty in
$V_{\text{obs}}$, is shown as a magenta bar in the lower right.

###  3.1 Modelling the RAR

I fit two functions to the RAR. The first is the “Simple interpolating
function (IF)” (Famaey & Binney, 2005):

$g_{\text{obs}}^{\text{pred}}=g_{\text{bar}}/2+\sqrt{g_{\text{bar}}^{2}/4+g_{\text{bar}}\>a_{0}}.$
(6)

Although in tension with Solar System measurements this function is highly
successful for galaxy dynamics (Famaey & McGaugh, 2012), and may readily be
tweaked to circumvent local constraints without appreciably altering its
larger-scale behaviour. One such modification is the “RAR IF” of Lelli et al.
(2017), which I have checked yields almost identical results to the Simple IF.

The reason I use the Simple IF is that the second function I consider is
designed to reduce to it in the zero-external-field limit. This is the EFE
formula for the nonrelativistic AQUAL (Bekenstein & Milgrom, 1984) theory of
MOND designed in Chae & Milgrom (2022):

$\displaystyle
g_{\text{obs}}^{\text{pred}}=g_{\text{bar}}\left(\frac{1}{2}+\left(\frac{1}{4}+\left(\left(\frac{g_{\text{bar}}}{a_{0}}\right)^{2}+(1.1e_{\text{N}})^{2}\right)^{-\frac{1}{2}}\right)^{\frac{1}{2}}\right)\times$
$\displaystyle\Biggl{(}1+\tanh\left(\frac{1.1e_{\text{N}}}{g_{\text{bar}}/a_{0}}\right)^{1.2}\times\left(-\frac{1}{3}\right)\times$
(7)
$\displaystyle\frac{\left(\left(\left(\frac{g_{\text{bar}}}{a_{0}}\right)^{2}+(1.1e_{\text{N}})^{2}\right)^{-\frac{1}{2}}\right)\left(\frac{1}{4}+\left(\left(\frac{g_{\text{bar}}}{a_{0}}\right)^{2}+(1.1e_{\text{N}})^{2}\right)^{-\frac{1}{2}}\right)^{-\frac{1}{2}}}{1+\left(\frac{1}{2}+2\left(\left(\frac{g_{\text{bar}}}{a_{0}}\right)^{2}+(1.1e_{\text{N}})^{2}\right)^{-\frac{1}{2}}\right)^{\frac{1}{2}}}\Biggr{)},$

where $e_{\text{N}}\equiv g_{\text{ext}}/a_{0}$ describes the strength of the
external field at the galaxy in question. The EFE arises in most formulations
of MOND due to its nonlinearity: the strong equivalence principle is violated
because the acceleration of a system as a whole cannot be transformed away in
calculation of its internal motions. This implies that otherwise identical
galaxies in different gravitational environments have different kinematics. A
stronger external field pushes the system towards the Newtonian regime by
reducing the gravitational boost of MOND, causing a downturn in the RAR at low
$g_{\text{bar}}$ where $g_{\text{ext}}$ can be a non-negligible fraction of
$g_{\text{obs}}$. While several fitting formulae for the EFE exist (e.g. Banik
& Zhao 2015; Haghi et al. 2019a; Zonoozi et al. 2021), Eq. 3.1 is the most
sophisticated in allowing for variable disk thickness and scale length—and the
orientation of the field relative to the disk axis through azimuthal
averaging—and has been shown to yield good agreement with the SPARC data (Chae
et al., 2022; Chae, 2022a).

I consider both the case of $e_{\text{N}}$ as a global parameter describing
the average external field over the sample, and as a parameter varying galaxy-
by-galaxy to describe their separate local environments. In the former case I
use a uniform prior sufficiently broad to enclose the full posterior; in the
latter, where there is insufficient information in the data for a meaningful
constraint on $e_{\text{N}}$, I impose a prior based on the environmental
field estimates of the SPARC galaxies from Desmond et al. (2018); Chae et al.
(2021). These are determined entirely independently of the SPARC data by
summing contributions to the gravitational field from the baryonic masses of
surrounding objects, including a sophisticated treatment of survey
incompleteness and other missing mass.

As my fiducial analysis I use the results assuming that missing baryons are
strongly clustered around visible objects (“maximum clustering”) because this
is expected in MOND and was shown in Chae et al. (2021, 2022) to give good
agreement with the SPARC rotation curves. I also consider an “average
clustering” model that assumes a prior distribution midway between the “max
clustering” and “no clustering” (missing baryons uncorrelated with visible
objects) results, with a width given by half the difference between the two.
This systematic uncertainty is larger than the statistical uncertainty in
either clustering case separately. The most precise calculation of Chae et al.
(2021) uses data from the Sloan Digital Sky Survey and hence is only valid
within the footprint of that survey, which includes 90 galaxies in my sample.
For the remaining 57 I take $\bar{e}_{\text{N}}$ to be the median
$e_{\text{N}}$ over all SPARC galaxies (in the corresponding clustering
model), with an uncertainty twice the median uncertainty for all SPARC
galaxies. This corresponds to a conservatively wide prior for galaxies without
object-specific prior information, while still leveraging information on the
$e_{\text{N}}$ distribution across the population. Combined with the no-EFE
case (Eq. 6), these EFE models ought roughly to span the space of possible EFE
behaviour and hence indicate the level of systematic uncertainty that the
unknown EFE behaviour induces.

Table 1 summarises the free parameters of the inference and their priors.

###  3.2 Inference procedure

The parameters inferred in the fiducial model are
$a_{0},\sigma_{\text{int}},147\times e_{\text{N}},147\times D,147\times
i,147\times
L_{3.6},147\times\Upsilon_{\text{disk}},147\times\Upsilon_{\text{gas}}$ and
$31\times\Upsilon_{\text{bulge}}$. At any point in parameter space I calculate
$g_{\text{bar}}$ and $g_{\text{obs}}$ as

$g_{\text{bar}}=\frac{\left(\Upsilon_{\text{gas}}V_{\text{gas}}^{2}+L_{36}/\bar{L}_{36}\>(\Upsilon_{\text{disk}}V_{\text{disk}}^{2}+\Upsilon_{\text{bulge}}V_{\text{bul}}^{2})\right)}{r},$
(8)
$g_{\text{obs}}=\frac{V_{\text{obs}}^{2}}{r}\>\frac{\sin(\bar{i})^{2}}{\sin(i)^{2}}\>\frac{\bar{D}}{D},$
(9)

where $V_{\text{gas}}$, $V_{\text{disk}}$, $V_{\text{bul}}$, $V_{\text{obs}}$
and $r$ are the velocities generated by the gas, disk and bulge, observed
velocity and galactocentric radius as quoted in the SPARC database (i.e.
assuming $D=\bar{D}$, $i=\bar{i}$, $L_{3.6}=\bar{L}_{3.6}$ and all
$\Upsilon=1$) respectively.

The only remaining uncertainty to treat as statistical is the contribution of
$\delta V_{\text{obs}}$ to $g_{\text{obs}}$.333This is a combination of a
formal error from the entire disk fit and a contribution from the difference
between the velocities of the approaching and receding sides of the disk (see
Lelli et al. 2016, eq. 1). These noise terms are not uncorrelated Gaussian
random variables, so a further improvement to the method would be to either
sample them or model their covariance structure (see Sec. 5). I assume this is
lognormal, so that

$\delta\log(g_{\text{obs}})=\sqrt{\log(1+(\delta
g_{\text{obs}}/g_{\text{obs}})^{2})}/\ln(10)$ (10)

where

$\delta g_{\text{obs}}/g_{\text{obs}}=2\>\delta
V_{\text{obs}}/V_{\text{obs}}.$ (11)

I then use either Eq. 6 or 3.1 to calculate the predicted $g_{\text{obs}}$ at
each $g_{\text{bar}}$. $\sigma_{\text{int}}$ simply adds in quadrature with
$\delta\log(g_{\text{obs}})$, so the likelihood is:

$\displaystyle\mathcal{L}(d|\vec{p})=$
$\displaystyle\prod_{j}\frac{1}{\sqrt{2\pi\sigma_{\text{tot}}^{2}}}\times$
(12)
$\displaystyle\exp\\{-(\log(g_{\text{obs,j}})-\log(g_{\text{obs,j}}^{\text{pred}}))^{2}/(2\>\sigma_{\text{tot,j}}^{2})\\}$

where $\vec{p}$ is the parameter vector and

$\sigma_{\text{tot,j}}^{2}\equiv\delta\log(g_{\text{obs,j}})^{2}+\sigma_{\text{int}}^{2}.$
(13)

$j$ runs over the $2696$ data points.

This fiducial analysis assumes no statistical uncertainty on the velocities
sourced by the gas, disk and bulge. However the calculation of
$V_{\text{disk}}$, $V_{\text{bul}}$ and $V_{\text{gas}}$ in Lelli et al.
(2016) made assumptions about the 3D geometry of these baryons, particularly
in the thickness of the disk components. Variation may be expected to alter
the baryon velocities at the $\sim$10-15 per cent level (F. Lelli, priv.
comm.). I therefore also consider models in which these velocities are each
assigned 10 per cent uncorrelated Gaussian uncertainties,444These assumptions
are unlikely to hold in detail because, as with the other nuisance parameters,
variation in disk thickness, disk flaring or the oblateness of bulges will
cause correlated deviations across the rotation curves. The effects will be
larger at smaller $r$ where higher order multipoles of the potential are more
important for the velocity field. My leading-order assumption is meant merely
to assess the characteristic impact of uncertainties of this magnitude on the
RAR parameters, especially $\sigma_{\text{int}}$. which are propagated
according to

$\displaystyle\delta V_{\text{bar}}^{2}=$
$\displaystyle(0.2\>\Upsilon_{\text{gas}}\>V_{\text{gas}}|V_{\text{gas}}|)^{2}+(0.2\>\Upsilon_{\text{disk}}\>V_{\text{disk}}^{2})^{2}$
(14) $\displaystyle+(0.2\>\Upsilon_{\text{bulge}}\>V_{\text{bul}}^{2})^{2}$
$\delta\log(g_{\text{bar}})=\sqrt{\log(1+(\delta
V_{\text{bar}}^{2}/V_{\text{bar}}^{2})^{2})}/\ln(10).$ (15)

The application of these equations will be indicated by “boosted
uncertainties”.

In this case, the presence of uncertainties in the $x$ direction of the RAR
plane introduces latent variables describing the true position of each point
on the $x$-axis in the Bayesian hierarchical model. This makes the likelihood
function for the parameters of interest alone ambiguous. Two approaches to
remove the latent nuisance parameters without sampling them are to marginalise
over them with a uniform prior, or to maximise the likelihood with respect to
each of them (as a function of the other parameters in the inference) to
produce a profile likelihood for the other parameters. These result in
different maximum-likelihood points and parameter constraints. Tests on mock
data (in agreement with literature results; Berger et al. 1999; Hadzhiyska et
al. 2023) show that the marginalised likelihood recovers the correct intrinsic
scatter and weakly biased shape parameters (e.g. $a_{0}$ and $e_{\text{N}}$)
while the profile likelihood recovers unbiased shape parameters but can bias
$\sigma_{\text{int}}$ significantly low. As I am mainly interested in whether
the boosted uncertainties allow for an intrinsic-scatter-free RAR I opt for
the marginalised likelihood, which replaces Eq. 13 by

$\sigma_{\text{tot,j}}^{2}=\delta\log(g_{\text{obs,j}})^{2}+\sigma_{\text{int}}^{2}+\frac{\text{d}\log
g_{\text{obs}}^{\text{pred}}}{\text{d}\log
g_{\text{bar}}}\bigg{|}_{g_{\text{bar,j}}}^{2}\delta\log(g_{\text{bar,j}})^{2}$
(16)

in Eq. 12 (for the derivation see e.g. Desmond et al. 2023 sec. 3.2).

I perform the inference using the No U-Turns Sampler (NUTS; Hoffman & Gelman
2011) method of Hamiltonian Monte Carlo (HMC), as implemented in numpyro (Phan
et al., 2019; Bingham et al., 2019). I initialise the sampler to the median of
20,000 points randomly drawn from the prior, which I find to yield good
convergence behaviour. For each inference I concatenate 28 separate chains run
in parallel, manually tuning the number of warmup and sampling steps to ensure
that burn-in is complete and that there are enough effective samples for the
Gelman-Rubin statistic (Gelman & Rubin, 1992) to satisfy $|r-1|<0.001$. This
requires $\sim$1000 warmup steps, $\sim$4000 sampling steps and takes $\sim$1
hour per model to run.

Parameter | Definition | Prior
---|---|---
$a_{0}$ | Acceleration constant ($10^{-10}$ ms-2) | $\mathcal{U}(0.1,5)$
$\sigma_{\text{int}}$ | Intrinsic scatter in $g_{\text{obs}}$ (dex) | $\mathcal{U}(0,1)$
$e_{\text{N}}$ | External field strength relative to $a_{0}$ |  $\mathcal{U}(0,0.5)$ or Lognormal($\log(\bar{e}_{\text{N}}),\delta\log(e_{\text{N}})$)1
$D$ | Distance | $\mathcal{N}\text{t}(\bar{D},\delta D;0,\text{None})$
$i$ | Inclination | $\mathcal{N}\text{t}(\bar{i},\delta i;0,180\deg)$
$L_{3.6}$ | Luminosity at 3.6$\mu$m | $\mathcal{N}(\bar{L}_{3.6},\delta L_{3.6})$
$\Upsilon_{\text{disk}}$ |  $M$/$L$ of disk | Lognormal(-0.72346, 0.24622)2
$\Upsilon_{\text{bulge}}$ |  $M$/$L$ of bulge | Lognormal(-0.38699, 0.24622)2
$\Upsilon_{\text{gas}}$ |  $M$/$L$ of gas | Lognormal($\log(\bar{\Upsilon}_{\text{gas}}(M_{*}))-0.5\log(1.01)$, 0.0997513)2

1 See Sec. 3.1 2 See footnote 2
---

Table 1: The free parameters of the inference and their priors. $\mathcal{N}$t
denotes a truncated normal with lower and upper limits given by the final two
arguments.

###  3.3 Validation with mock data

My analysis uses uniform priors on $a_{0}$, $\sigma_{\text{int}}$ and global
$e_{\text{N}}$, which are not reparametrisation invariant and cannot, pace
popular opinion, be considered uninformative. Such priors are prone to
contributing volume effects to the posterior, which can lead to significant
biases when applied to parameters to which the likelihood is relatively
insensitive (e.g. Hadzhiyska et al. 2023). In addition, the uncertainties and
finite sample size lead to scatter in the maximum-likelihood parameters around
the population values (sample variance). To assess the impact of these effects
I analyse mock data generated by the following procedure:

  1. (1)

Randomly sample the galaxy parameters from their prior distributions

  2. (2)

Rescale the “observed” $g_{\text{bar}}$ according to Eq. 8 to calculate
$g_{\text{bar}}^{\text{true}}$

  3. (3)

Use either the no-EFE, global-EFE or max-clustering-EFE model, with some true
$a_{0}$, $e_{\text{N}}$, $\sigma_{\text{int}}$, to calculate
$g_{\text{obs}}^{\text{true}}$ from the $g_{\text{bar}}^{\text{true}}$

  4. (4)

Transform to observed $g_{\text{obs}}$ through Eq. 9, and hence to
$V_{\text{obs}}$ assuming the same $r$ values in the mock data as in the SPARC
data

  5. (5)

Replacing the SPARC $V_{\text{obs}}$ by these values, calculate the maximum-
likelihood values of all parameters and run the inference to compute their
posteriors

  6. (6)

Repeat twice with different random seeds.

Choosing $a_{0}=1.2$, $\sigma_{\text{int}}=0.05$ dex and $e_{\text{N}}=0.01$
in the case with global EFE, the results are shown in Table 2. I find neither
the maximum-likelihood parameters nor their posteriors to be significantly
different to their true values, showing the above effects not to be important
for the data and models under consideration. The $e_{\text{N}}$ values and
galaxy parameters (not shown) are similarly unbiased. This gives confidence to
proceed with the analysis of the real data.

| $a_{0}\>(1.2)$ | $\sigma_{\text{int}}\>(0.05)$
---|---|---
Model | ML | numpyro | ML | numpyro
No EFE | 1.21 | $1.21\pm 0.04$ | 0.051 | $0.051\pm 0.001$
| 1.19 | $1.18\pm 0.04$ | 0.049 | $0.049\pm 0.001$
| 1.21 | $1.15\pm 0.04$ | 0.050 | $0.050\pm 0.001$
Global EFE | 1.20 | $1.29\pm 0.05$ | 0.050 | $0.050\pm 0.001$
| 1.20 | $1.17\pm 0.04$ | 0.050 | $0.050\pm 0.001$
| 1.22 | $1.18\pm 0.04$ | 0.050 | $0.049\pm 0.001$
Max-clust EFE | 1.21 | $1.18\pm 0.04$ | 0.050 | $0.051\pm 0.001$
| 1.18 | $1.14\pm 0.04$ | 0.051 | $0.051\pm 0.001$
| 1.20 | $1.24\pm 0.04$ | 0.049 | $0.050\pm 0.001$
Table 2: RAR parameters inferred from mock data generated according to the
prescription of Sec. 3.3. “ML” stands for maximum likelihood and bracketed
numbers after the parameters show the generating values. The HMC uncertainties
are $1\sigma$.

##  4 Results

Table 3 shows the median and $1\sigma$ uncertainty of the RAR parameters for
each of the models considered, along with their maximum log-likelihood
($\ln(\hat{L})$), maximum log-posterior ($\ln(\hat{P})$) and Bayesian
information criterion (BIC) relative to the first model. The BIC should be
taken as a very rough estimator only for the Bayesian evidence, both because
the number of data points does not greatly exceed the number of parameters and
because the parameter priors are not necessarily slowly varying at the maximum
a-posteriori point. For example, replacing $2\ln(\hat{L})$ by $2\ln(\hat{P})$
in the BIC formula changes $\Delta$BIC to 1.44 for the “No scatter, global
EFE” model, turning “decisive” evidence on the Jeffrey’s scale against the
inclusion of $e_{\text{N}}$ to “barely worth mentioning.” Regardless of this,
a clear result from the goodness-of-fit statistics is that either non-zero
intrinsic scatter or boosted uncertainties is strongly preferred, mainly
through a large increase in the likelihood. The fiducial model, which I
considered a priori to be most likely, is “Scatter, max-clustering EFE”.

Model | $\mathbf{a_{0}}$ | $\mathbf{e_{\text{N}}}$ | $\mathbf{\sigma_{\text{int}}}$ | $\Delta\ln(\hat{L})$ | $\Delta\ln(\hat{P})$ |  $\Delta$BIC
---|---|---|---|---|---|---
No scatter, no EFE | $1.134^{+0.028}_{-0.027}$ | — | — | 0 | 0 | 0
No scatter, global EFE | $1.138^{+0.028}_{-0.027}$ | $0.0016^{+0.0005}_{-0.0005}$ | — | $-3.44$ | 3.24 | 14.8
No scatter, max-clustering EFE | $1.309^{+0.039}_{-0.037}$ | $0.0050^{+0.0203}_{-0.0033}$ | — | 292 | 55.8 | 577
No scatter, avg-clustering EFE | $1.307^{+0.038}_{-0.036}$ | $0.0022^{+0.0231}_{-0.0017}$ | — | 275 | $-19.2$ | 611
Scatter, no EFE | $1.070^{+0.032}_{-0.031}$ | — | $0.035^{+0.001}_{-0.001}$ | 2020 | 2381 | $-4032$
As above, boosted uncertainties | $1.121^{+0.035}_{-0.034}$ | — | $<0.0020$ | 1936 | 2354 | $-3864$
Scatter, global EFE | $1.077^{+0.033}_{-0.032}$ | $0.0017^{+0.0009}_{-0.001}$ | $0.035^{+0.001}_{-0.001}$ | 2019 | 2375 | $-4023$
Scatter, max-clustering EFE | $1.236^{+0.043}_{-0.041}$ | $0.0049^{+0.0125}_{-0.003}$ | $0.033^{+0.001}_{-0.001}$ | 2052 | 2220 | $-2935$
As above, boosted uncertainties | $1.275^{+0.047}_{-0.044}$ | $0.0047^{+0.0103}_{-0.0029}$ | $<0.0019$ | 1948 | 2182 | $-2727$
Scatter, avg-clustering EFE | $1.272^{+0.047}_{-0.045}$ | $0.0020^{+0.0134}_{-0.0016}$ | $0.033^{+0.001}_{-0.001}$ | 2074 | 2161 | $-2978$
Table 3: Constraints on RAR parameters and goodness-of-fit statistics for the
models considered. For the maximum-clustering and average-clustering EFE
models the quoted $e_{\text{N}}$ constraints describe the stacked posteriors
over all galaxies. The final three columns are the maximum log-likelihood,
maximum log-posterior and Bayesian information criterion relative to the first
model.

The value of $a_{0}$ is slightly reduced if the RAR is assumed to possess
intrinsic scatter, and more significantly increased by including the EFE with
external field strength priors from the baryonic large-scale structure. This
is because these priors somewhat increase $e_{\text{N}}$ relative to
likelihood alone and there is a positive degeneracy between $e_{\text{N}}$ and
$a_{0}$. Differences in $a_{0}$ between the models are up to a few times
larger than their statistical uncertainties, and a naive averaging over all
the models implies $a_{0}=1.19\pm 0.04\,\text{(stat)}\pm 0.09\,\text{(sys)}$.

The intrinsic scatter is similar in all models in which it is included, except
in the “boosted uncertainties” case. The results suggest that
$\sigma_{\text{int}}=0.034\pm 0.001\,\text{(stat)}\pm 0.001\,\text{(sys)}$
dex, with the significant caveat that this assumes the SPARC error model is
reliable. A 10 per cent uncertainty on $V_{\text{disk}}$, $V_{\text{bul}}$ and
$V_{\text{disk}}$—as may be expected from deviations from the assumed 3D
baryon geometry—is sufficient to set the preferred $\sigma_{\text{int}}$ to 0.
The intrinsic scatter is not driven by the most egregious outliers: removing
the 11 points with $g_{\text{bar}}>g_{\text{obs}}$ and $g_{\text{obs}}<2$
after transformation for the model without EFE (see Fig. 1) reduces it only to
0.031 dex.

There is weak evidence for $\langle e_{\text{N}}\rangle>0$ when it is given a
uniform prior, as evidenced both by the inferred $e_{\text{N}}$ being
consistent with 0 within $\sim$3$\sigma$ and by the small gain in maximum
posterior value across the chain ($\Delta\ln(\hat{P})$). Adding this parameter
is not favoured by the BIC. The model perhaps most similar to the MOND
expectation with global $e_{\text{N}}$ (not shown in Table 3) has
$\sigma_{\text{int}}=0$ and boosted uncertainties: in this case $\langle
e_{\text{N}}\rangle=0.0024\pm 0.0009$, again a 3$\sigma$ “detection” but a
slightly larger value. This model has $\Delta$BIC$=-3860$, very similar to the
case with $\sigma_{\text{int}}$ and boosted uncertainties but without EFE,
showing that again the addition of a global $e_{\text{N}}$ is not favoured.

The average prior values of $e_{\text{N}}$ over all the galaxies for the
maximum-clustering and average-clustering cases are 0.0050 and 0.0018
respectively, while the posteriors average to 0.0050 and 0.0022 without
$\sigma_{\text{int}}$ and 0.0049 and 0.0020 with. This shows that when
modelled galaxy-by-galaxy the data does not disfavour significant values of
$e_{\text{N}}$ in agreement with the large-scale structure expectation, and in
fact in both cases the maximum-likelihood value is increased over the case of
no or global $e_{\text{N}}$, significantly without intrinsic scatter and
moderately with it. Prior evidence is however required to support the presence
of the EFE in a model comparison sense as the BIC strongly disfavours the
addition of 147 $e_{\text{N}}$ parameters. While the average-clustering EFE
model gives a slightly higher $\hat{L}$ when including $\sigma_{\text{int}}$,
it gives a lower $\hat{L}$ without it and a lower $\hat{P}$ in both cases.
This constitutes weak evidence in favour of the maximum-clustering prior, but
the data is far from sufficient to distinguish robustly between them.

Fig. 1 uses Eqs. 8 and 9 to transform $g_{\text{bar}}$ and $g_{\text{obs}}$
according to the parameters at the median of the posterior for the inference
with $\sigma_{\text{int}}$ but without EFE (the model preferred by the BIC),
including the $2\sigma$ model prediction. This illustrates the extreme
tightness of the underlying (posterior) RAR, including relative to the
traditional prior RAR (e.g. Desmond et al. 2023, fig. 1). The few outlying
points are also outliers of the prior RAR, and are brought slightly closer to
the line by the transformation.

Fig. 2 shows an excerpt from the corner plot of the inference that also
includes the EFE with a global $e_{\text{N}}$. The posteriors on three
parameters for galaxy 1 (D512-2) are compared to their maximum prior values
(blue lines), with which they agree well. $a_{0}$, $e_{\text{N}}$ and
$\sigma_{\text{int}}$ are not strongly degenerate with other parameters, while
the full parameter space exhibits the degeneracies expected from Eqs. 8 and 9.
The median distance to the SPARC galaxies is reduced by 1.7 per cent going
from the inference with $\sigma_{\text{int}}$ but without EFE to that with
local $e_{\text{N}}$ and maximum-clustering prior, although within the latter
inference $D$ is positively correlated with $e_{\text{N}}$ on average within
the chain. This is because at higher $e_{\text{N}}$,
$g_{\text{obs}}^{\text{pred}}$ is lower at low $g_{\text{bar}}$ and hence
closer to $g_{\text{obs}}$ at higher $D$ (or $i$).

To show more generally the differences between the prior and posterior values
of the galaxy parameters, Fig. 3 shows the distribution of normalised
residuals for the model with intrinsic scatter and galaxy-by-galaxy
$e_{\text{N}}$ with maximum-clustering prior. For galaxy parameter $X$, the
normalised residual is defined as

$(\text{med}(X_{\text{post}})-\bar{X})/\text{std}(X_{\text{post}})$ (17)

where subscript “post” and “prior” denote the posterior and prior
distributions, and “med” and “std” stand for median and standard deviation.
Larger values of the residual indicate a posterior significantly shifted from
the prior due to the influence of the likelihood. The width of the
distributions therefore reflect the sensitivity of the data to the parameters:
those with distributions sharply peaked at 0 such as $L_{3.6}$ and
$\Upsilon_{\text{gas}}$ are relatively unimportant so that
$\text{med}({X}_{\text{post}})\approx\bar{X}$, while those with very broad
distributions have likelihood often peaked far from the prior centres. This is
especially pronounced for $\Upsilon_{\text{bulge}}$ and
$\Upsilon_{\text{disk}}$ which determine $g_{\text{bar}}$ in the galaxies’
central regions. $\Upsilon_{\text{bulge}}$ in particular has a slight negative
offset to reduce the number of inner points for which
$g_{\text{bar}}>g_{\text{obs}}$. The black dashed lines are what one would
expect if $X_{\text{post}}$ scatters around $\bar{X}$ with a Gaussian
distribution of width given by std($X_{\text{post}}$).

Finally, I show in Fig. 4 smoothed distributions of the fractional
uncertainties in the galaxy-specific parameters in the prior and posterior for
the inference with intrinsic scatter and local $e_{\text{N}}$ with maximum-
clustering prior. In most cases the RAR constraint has increased the precision
with which the parameters are known. This is especially marked for the
distance, where galaxies in the second mode of the prior distribution (those
at $\lesssim$60 Mpc with only redshift distances; Lelli et al. 2016) are
brought into a single posterior mode at $\sim$10 per cent uncertainty. This
illustrates the utility of the RAR as a redshift-independent distance
measurement method. There is analogous behaviour for $e_{\text{N}}$, where the
higher prior mode corresponds to galaxies outside the SDSS footprint which are
assigned higher uncertainties (see Sec. 3.1); these become better known on
applying the RAR constraint. The uncertainties on $\Upsilon_{\text{disk}}$ and
$\Upsilon_{\text{bulge}}$ also fall markedly, partly to reduce the number of
inner points with posterior probability at
$g_{\text{bar}}>g_{\text{obs}}$.555An alternative explanation for apparent
$g_{\text{bar}}>g_{\text{obs}}$ at low $r$ is the presence of bars or
asymmetries in the inner regions of the disk. In some cases it may be
preferable to excise such data (e.g. Katz et al. 2014) to prevent it from
biasing $\Upsilon_{\text{disk}}$ low, which has a knock-on effect across the
rotation curve. The analogues of Figs. 1–4 for the other models are
qualitatively similar, in line with the variation in their results shown in
Table 3.

Supplementary tables available at https://zenodo.org/record/7752545 (Desmond,
2023) contain the mean, median and 1, 2 and 3$\sigma$ constraints on the
parameters for each of the models considered. The format is illustrated in
Table 4, which shows the first five and last two rows for the model with
$\sigma_{\text{int}}$ and global $e_{\text{N}}$. An example use of these is to
resolve distance variations among galaxies in the Ursa Major cluster, all of
which have $\bar{D}=18$ Mpc. A reference table contains the galaxy names at
each index as well as the bulge luminosity and means and standard deviations
of their Gaussian priors on $D$, $i$ and $L_{3.6}$, transcribed from the SPARC
database. Plots of the prior and posterior constraints on $D$, $i$, $L_{3.6}$,
$e_{\text{N}}$, $\Upsilon_{\text{disk}}$, $\Upsilon_{\text{bulge}}$ and
$\Upsilon_{\text{gas}}$ for all galaxies under each model are also included.

Parameter | 0.135% | 2.275% | 15.87% | 50% | mean | 84.13% | 97.725% | 99.865%
---|---|---|---|---|---|---|---|---
$a_{0}$ | 0.9847 | 1.0145 | 1.0455 | 1.0773 | 1.078 | 1.1102 | 1.145 | 1.1792
$\sigma_{\text{int}}$ | 0.0321 | 0.033 | 0.0339 | 0.0349 | 0.0349 | 0.0358 | 0.0368 | 0.0379
$e_{\text{N}}$ | 0.0 | 0.0001 | 0.0007 | 0.0017 | 0.0017 | 0.0026 | 0.0035 | 0.0044
Dist[0] | 2.0319 | 2.2868 | 2.5538 | 2.8242 | 2.8252 | 3.0978 | 3.3688 | 3.6426
Dist[1] | 6.6289 | 8.4419 | 10.8253 | 13.8181 | 14.0917 | 17.3551 | 21.3172 | 25.6289
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮
ML_gas[145] | 0.6955 | 0.7658 | 0.8427 | 0.9271 | 0.9308 | 1.0185 | 1.1198 | 1.2295
ML_gas[146] | 0.7595 | 0.832 | 0.9139 | 1.0046 | 1.009 | 1.1038 | 1.2145 | 1.3344
Table 4: Excerpt of posterior summaries for the model with intrinsic scatter
and global $e_{\text{N}}$. The full tables for all models are available
online. Figure 2: Corner plot of selected parameters from the inference
including intrinsic scatter and global $e_{\text{N}}$. For the distance (in
Mpc), inclination (in degrees) and disk mass-to-light ratio of the first
galaxy (D512-2), the maximum a priori values are shown by the blue lines.
Figure 3: The distribution of normalised residuals (Eq. 17) of the galaxy-
specific parameters in the inference including intrinsic scatter and local
$e_{\text{N}}$ with maximum-clustering prior. The mirroring of some of the
curves about the $x$-axis is for visual clarity only. A standard normal
distribution—corresponding to $X_{\text{post}}$ scattering around $\bar{X}$ as
a Gaussian with width given by std($X_{\text{post}}$)—is shown in dashed
black. Figure 4: Distribution of prior and posterior fractional uncertainties
on galaxy-specific parameters for the inference with scatter and local
$e_{\text{N}}$ with maximum-clustering prior. In the lower panels the vertical
dashed lines are the fractional prior uncertainties (25 and 10 per cent),
which are the same for all galaxies.

##  5 Discussion

While the value of $a_{0}$ that I infer ($1.1\lesssim a_{0}\lesssim 1.3$) is
in full agreement with literature results, the RAR intrinsic scatter
$\sigma_{\text{int}}\approx 0.034$ dex is significantly smaller. Lelli et al.
(2017) quote a total scatter of 0.13 dex and argue that most of this comes
from observational uncertainties; using their model and Eqs. 12 and 16 I find
$\sigma_{\text{int}}=0.082\pm 0.003$ dex. Li et al. (2018) quote an intrinsic
scatter of 0.057 dex from residuals around their best-fit relation with fixed
$a_{0}=1.2$. A reduction in intrinsic scatter when marginalising over galaxy
parameters is not guaranteed: although the denominator of Eq. 12 favours
smaller $\sigma_{\text{int}}$, this is offset by lower prior probabilities of
the galaxy parameters at values that bring the points closer to the
theoretical line to reduce the exponent in Eq. 12. The statistical
uncertainties on $g_{\text{bar}}$ and $g_{\text{obs}}$ also fall greatly (the
former to 0) when the galaxy parameters are inferred, so a lower total scatter
of the points around the line need not translate into lower
$\sigma_{\text{int}}$. That the preferred intrinsic scatter in the full
analysis is extremely small hints towards the RAR being at base a practically
monotonic correlation.

The question of the fundamentality of the RAR has important ramifications for
the mass discrepancy problem on galaxy scales. Lelli et al. (2017) claimed the
relation to be law-like, a result supported here. It will be challenging for
simulations or semi-analytic models in $\Lambda$CDM to reproduce the tightness
of the underlying RAR given that even the scatter of the prior RAR, $\sim$0.1
dex, is non-trivial (Di Cintio & Lelli, 2016; Desmond, 2017; Keller & Wadsley,
2017; Ludlow et al., 2017). It is important to bear in mind however that SPARC
is only a small fraction of the data pertinent to the RAR: one may also use
ultra-diffuse galaxies (e.g. Freundlich et al. 2022), local dwarf spheroidals
(e.g. McGaugh & Wolf 2010; McGaugh & Milgrom 2013), low-acceleration regions
including the outer Milky Way (Oman et al., 2020), stacked weak lensing
(Brouwer et al., 2021) and groups and clusters of galaxies (Chae et al., 2019;
Chae et al., 2020a; Chan & Del Popolo, 2020; Tian et al., 2020; Pradyumna &
Desai, 2021; Gopika & Desai, 2021). Some of this data appears to deviate from
the MOND expectation. SPARC is however one of the few datasets with
uncertainties under sufficient control for the present analysis to be feasible
and meaningful. Future work may extend it to new regimes.

Naively one expects $\sigma_{\text{int}}=0$ in MOND, but in fact this is only
true in modified inertia formulations in the limit of perfectly circular
orbits. Even in such models a scatter is introduced by deviations from
circularity because the dynamics of an object depends on its entire past
trajectory (Milgrom, 2011, 2022). This may already be sufficient to account
for the 0.034 dex (8 per cent) scatter present in the underlying RAR.
Additional scatter is present in modified gravity formulations where the
algebraic MOND relation holds only in spherical symmetry (Famaey & McGaugh,
2012), a condition clearly violated in disk galaxies. It would be interesting
to quantify the scatter introduced by these effects in SPARC-like galaxies,
further testing the MOND paradigm and providing a novel way to distinguish
between the modified gravity and modified inertia interpretations (Petersen &
Lelli, 2020; Chae, 2022b).

The presence or absence of the EFE is controversial within the MOND
literature, with some studies claiming strong evidence for it (McGaugh &
Milgrom, 2013; Haghi et al., 2019b; Chae et al., 2020b) and others strong
evidence against (Hernandez et al., 2019; Freundlich et al., 2022). My work
does not resolve this issue: there is weak evidence for a positive average
external field strength across the sample, while inferring it galaxy-by-galaxy
with a prior from independent measurements of environment improves the
likelihood but is not favoured by the Bayesian information criterion. An
important caveat is that the fitting formula I use (Eq. 3.1) was only designed
for the outer regions of rotation curves (Chae & Milgrom, 2022)—while I am
applying it to them in their entirety—and is only valid within the AQUAL
model. Different EFE formulae may give different constraints on
$e_{\text{N}}$, and hence $a_{0}$ with which it is degenerate, but would not
be expected to alter $\sigma_{\text{int}}$ appreciably. The environmental
priors are also highly uncertain due to the possibility of clustered unseen
baryonic mass. Further work on modelling the EFE and constraining the external
field is therefore required to reach a definitive conclusion concerning the
existence of the EFE in galaxy dynamics.

Besides calibrating the RAR my work provides strong constraints on the
properties of the SPARC galaxies under the assumption that the underlying RAR
is as I model it. Summaries of these constraints are made public in online
tables to facilitate future studies using SPARC. One such application is to
use the RAR as a redshift-independent distance probe. The current study
calibrates the relation using a sample with informative distance priors. The
distance to any galaxy with a (partially) resolved rotation curve may be
inferred by fitting it to the calibrated RAR, including marginalisation over
the other relevant properties of the galaxy but not necessarily the parameters
of the RAR itself. This is analogous to the well-established Tully–Fisher and
Fundamental Plane methods, but achieves higher precision ($\sim$10 per cent
uncertainty on $D$ rather than 20-25 per cent; Tully et al. 2023) at the cost
of requiring resolved kinematics. My analysis also supplies enhanced kinematic
inclinations, as well as constraints on mass-to-light ratios which may be
correlated with other galaxy properties to advance understanding of stellar
populations and galaxies’ gas content.

Outliers of the underlying RAR may either be individual rotation curve points
with large residuals from the best-fit line (e.g. in Fig. 1) or entire
galaxies with parameters strongly shifted from their prior centres to achieve
a good fit (readily visible in the supplementary figures). Studying these on a
case-by-case basis may help identify peculiar galactic features and signpost
the need for more sophisticated modelling. For example, NGC 2915 seems to
require very high $\Upsilon_{\text{disk}}\approx 1.2$ in all models. This is a
starburst dwarf galaxy with significant radial motion in the inner regions,
the modelling of which affects the entire rotation curve. It has a highly
complex structure and is likely not in dynamical equilibrium towards the
centre, while further out there is a strong warp (Meurer et al., 1996; Elson
et al., 2010, 2011a, 2011b; Tang et al., 2022). The ability to correlate such
properties with the results of this analysis across the sample would lend
weight to the interpretation of the underlying RAR as fundamental.

Provided the fitting functions are good, the constraints on galaxy parameters
do not assume MOND any more than modelling the Tully–Fisher relation as a
power-law. Both may be thought of as empirical descriptions of the data
without consideration of their theoretical significance. Recently, however,
Desmond et al. (2023) have challenged the optimality of MOND functions (those
with Newtonian and deep-MOND or EFE-driven regimes) for fitting the RAR,
finding that the majority of functions that most efficiently compress the
SPARC data tend to constant $g_{\text{obs}}$ at low $g_{\text{bar}}$. These
functions and their parameters have no known theoretical significance, but
would lead to different “underlying” relations and hence galaxy parameter
constraints. This systematic uncertainty in e.g. distance measurements could
be assessed by repeating the present inference with one or more of these
functions. It is however not known how the results of Desmond et al. would be
affected by marginalising over the galaxy parameters separately for each
function as done here, which would alter the function ranking, or by extending
Exhaustive Symbolic Regression to higher complexity.

I see three technical ways in which this inference could be improved. The
first is to separate $\delta V_{\text{obs}}$ into a truly statistical and
systematic part, including its full covariance across the galaxies’ rotation
curves. The average $\delta\log(V_{\text{obs}})$ across the sample is 0.03
dex, so more accurate modelling of this has the potential to alter the best-
fit value of $\sigma_{\text{int}}$, which is similar. In particular, an
important contribution to $\delta V_{\text{obs}}$ is from the difference in
velocity between the approaching and receding sides of the disk, which is
likely strongly correlated over $r$. This may itself be correlated with
inclination, which in detail may vary across the disk. The second is to
characterise better the molecular gas and 3D baryon geometry. The former could
be done by estimating (or sampling) $M_{H_{2}}$ from the $M_{*}-M_{H_{2}}$
relation of a population of similar galaxies and adding a contribution to
$V_{\text{bar}}$ from a thin $H_{2}$ disk, and the latter by solving Poisson’s
equation for different assumptions about disk thickness, bulge oblateness and
potential asymmetries. For the SPARC galaxies these uncertainties are however
highly subdominant to those that I model explicitly. The third is to use a
Jeffrey’s rather than uniform prior for $a_{0}$, $e_{\text{N}}$ and
$\sigma_{\text{int}}$, thus eliminating potential volume effects. The mock
tests of Sec. 3.3 show these not to bias the results significantly. A
Jeffrey’s prior would however enable the inference of $e_{\text{N}}$ galaxy-
by-galaxy without importing large-scale structure information; when I tried
this on mock data using a uniform prior I found $a_{0}$ to be biased high,
presumably due to the allowed and poorly-constrained volume towards high
$e_{\text{N}}$. Another extension would be to try a more complex intrinsic
scatter model; Li et al. (2018) for example find a superposition of two
Gaussians to fit the residuals better than one, perhaps reflecting the
separate formal error and kinematic asymmetry contributions to
$V_{\text{obs}}$.

##  6 Conclusion

I have uncovered underlying RARs in the SPARC data by fitting the parameters
of the Simple interpolating function, with and without intrinsic scatter and
the external field effect, simultaneously with all relevant galaxy properties.
The preferred intrinsic scatter is very small, $0.034\pm 0.002$ dex, and
additional plausible uncertainties are capable of reducing this to 0. The
acceleration constant is in the range $1.1\lesssim a_{0}\lesssim 1.3$, in good
agreement with literature results. I find weak evidence for the external field
effect using an average external field strength over the full sample with a
uniform prior. Allowing the field strength to vary galaxy-by-galaxy with a
prior from large-scale structure observations improves the overall likelihood
of the data but is not favoured by the Bayesian information criterion. My
results suggest near-monotinicity and a high degree of regularity in the the
RAR, providing a fresh challenge to galaxy formation models. The constraints I
produce on the SPARC galaxies’ parameters—distance, inclination, luminosity
and disk, bulge and gas mass-to-light ratios—are the most precise to date
(although subject to systematic error if the underlying RAR is not as I model
it). I publicly release summaries of all posteriors for analyses that may
benefit from this information.

##  7 Data availability

The mean, median and 1, 2 and 3$\sigma$ confidence intervals of the parameters
for all models are available at https://zenodo.org/record/7752545. The
remaining data generated here, including the full HMC chains to explore
degeneracies, will be made available on reasonable request.

## Acknowledgements

I thank Deaglan Bartlett, Kyu-Hyun Chae, Federico Lelli, Stacy McGaugh, James
Prideaux-Ghee, Richard Stiskalek and Tariq Yasin for useful discussions.

I am supported by a Royal Society University Research Fellowship (grant no.
211046). This project has received funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No 693024).

For the purpose of open access, I have applied a Creative Commons Attribution
(CC BY) licence to any Author Accepted Manuscript version arising.

## References

  * Banik & Zhao (2015) Banik I., Zhao H., 2015, arXiv e-prints, p. arXiv:1509.08457
  * Banik & Zhao (2022) Banik I., Zhao H., 2022, Symmetry, 14, 1331
  * Bekenstein & Milgrom (1984) Bekenstein J., Milgrom M., 1984, ApJ, 286, 7
  * Berger et al. (1999) Berger J. O., Liseo B., Wolpert R. L., 1999, Statistical Science, 14, 1
  * Bingham et al. (2019) Bingham E., et al., 2019, J. Mach. Learn. Res., 20, 28:1
  * Blanton et al. (2008) Blanton M. R., Geha M., West A. A., 2008, ApJ, 682, 861
  * Brouwer et al. (2021) Brouwer M. M., et al., 2021, A&A, 650, A113
  * Cappellari et al. (2013) Cappellari M., et al., 2013, MNRAS, 432, 1709
  * Chae (2022a) Chae K.-H., 2022a, arXiv e-prints, p. arXiv:2207.11069
  * Chae (2022b) Chae K.-H., 2022b, Astrophys. J., 941, 55
  * Chae & Milgrom (2022) Chae K.-H., Milgrom M., 2022, ApJ, 928, 24
  * Chae et al. (2019) Chae K.-H., Bernardi M., Sheth R. K., Gong I.-T., 2019, ApJ, 877, 18
  * Chae et al. (2020a) Chae K.-H., Bernardi M., Domínguez Sánchez H., Sheth R. K., 2020a, ApJ, 903, L31
  * Chae et al. (2020b) Chae K.-H., Lelli F., Desmond H., McGaugh S. S., Li P., Schombert J. M., 2020b, ApJ, 904, 51
  * Chae et al. (2021) Chae K.-H., Desmond H., Lelli F., McGaugh S. S., Schombert J. M., 2021, ApJ, 921, 104
  * Chae et al. (2022) Chae K.-H., Lelli F., Desmond H., McGaugh S. S., Schombert J. M., 2022, arXiv e-prints, p. arXiv:2209.07357
  * Chan & Del Popolo (2020) Chan M. H., Del Popolo A., 2020, MNRAS, 492, 5865
  * Desmond (2017) Desmond H., 2017, MNRAS, 464, 4160
  * Desmond (2023) Desmond H., 2023, doi:10.5281/zenodo.7752545
  * Desmond & Wechsler (2015) Desmond H., Wechsler R. H., 2015, MNRAS, 454, 322
  * Desmond & Wechsler (2017) Desmond H., Wechsler R. H., 2017, MNRAS, 465, 820
  * Desmond et al. (2018) Desmond H., Ferreira P. G., Lavaux G., Jasche J., 2018, MNRAS, 474, 3152
  * Desmond et al. (2023) Desmond H., Bartlett D. J., Ferreira P. G., 2023, arXiv e-prints, p. arXiv:2301.04368
  * Di Cintio & Lelli (2016) Di Cintio A., Lelli F., 2016, MNRAS, 456, L127
  * Djorgovski & Davis (1987) Djorgovski S., Davis M., 1987, ApJ, 313, 59
  * Dressler et al. (1987) Dressler A., Lynden-Bell D., Burstein D., Davies R. L., Faber S. M., Terlevich R., Wegner G., 1987, ApJ, 313, 42
  * Elson et al. (2010) Elson E. C., de Blok W. J. G., Kraan-Korteweg R. C., 2010, Monthly Notices of the Royal Astronomical Society, 404, 2061
  * Elson et al. (2011a) Elson E. C., de Blok W. J. G., Kraan-Korteweg R. C., 2011a, Monthly Notices of the Royal Astronomical Society, 411, 200
  * Elson et al. (2011b) Elson E. C., de Blok W. J. G., Kraan-Korteweg R. C., 2011b, Monthly Notices of the Royal Astronomical Society, 415, 323
  * Faber & Jackson (1976) Faber S. M., Jackson R. E., 1976, ApJ, 204, 668
  * 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
  * Freundlich et al. (2022) Freundlich J., Famaey B., Oria P.-A., Bílek M., Müller O., Ibata R., 2022, A&A, 658, A26
  * Gelman & Rubin (1992) Gelman A., Rubin D. B., 1992, Statistical Science, 7, 457
  * Gnedin et al. (2007) Gnedin O. Y., Weinberg D. H., Pizagno J., Prada F., Rix H.-W., 2007, ApJ, 671, 1115
  * Gopika & Desai (2021) Gopika K., Desai S., 2021, Physics of the Dark Universe, 33, 100874
  * Hadzhiyska et al. (2023) Hadzhiyska B., Wolz K., Azzoni S., Alonso D., García-García C., Ruiz-Zapatero J., Slosar A., 2023, arXiv e-prints, p. arXiv:2301.11895
  * Haghi et al. (2019b) Haghi H., et al., 2019b, MNRAS, 487, 2441
  * Haghi et al. (2019a) Haghi H., et al., 2019a, Monthly Notices of the Royal Astronomical Society, 487, 2441
  * Hernandez et al. (2019) Hernandez X., Cortés R. A. M., Allen C., Scarpa R., 2019, International Journal of Modern Physics D, 28, 1950101
  * Hoffman & Gelman (2011) Hoffman M. D., Gelman A., 2011, arXiv e-prints, p. arXiv:1111.4246
  * Katz et al. (2014) Katz H., McGaugh S. S., Sellwood J. A., de Blok W. J. G., 2014, MNRAS, 439, 1897
  * Keller & Wadsley (2017) Keller B. W., Wadsley J. W., 2017, ApJ, 835, L17
  * Lelli et al. (2016) Lelli F., McGaugh S. S., Schombert J. M., 2016, AJ, 152, 157
  * 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
  * Ludlow et al. (2017) Ludlow A. D., et al., 2017, Phys. Rev. Lett., 118, 161103
  * McGaugh (2004) McGaugh S. S., 2004, ApJ, 609, 652
  * McGaugh & Milgrom (2013) McGaugh S., Milgrom M., 2013, ApJ, 775, 139
  * McGaugh & Schombert (2014) McGaugh S. S., Schombert J. M., 2014, AJ, 148, 77
  * McGaugh & Wolf (2010) McGaugh S. S., Wolf J., 2010, ApJ, 722, 248
  * McGaugh et al. (2000) McGaugh S. S., Schombert J. M., Bothun G. D., de Blok W. J. G., 2000, ApJ, 533, L99
  * McGaugh et al. (2020) McGaugh S. S., Lelli F., Schombert J. M., 2020, Research Notes of the AAS, 4, 45
  * Meidt et al. (2014) Meidt S. E., et al., 2014, ApJ, 788, 144
  * Meurer et al. (1996) Meurer G. R., Carignan C., Beaulieu S., Freeman K. C., 1996, Astron. J., 111, 1551
  * Milgrom (1983a) Milgrom M., 1983a, ApJ, 270, 365
  * Milgrom (1983b) Milgrom M., 1983b, ApJ, 270, 371
  * Milgrom (1983c) Milgrom M., 1983c, ApJ, 270, 384
  * Milgrom (2011) Milgrom M., 2011, arXiv e-prints, p. arXiv:1111.1611
  * Milgrom (2022) Milgrom M., 2022, Phys. Rev. D, 106, 064060
  * Navarro et al. (2017) Navarro J. F., Benítez-Llambay A., Fattahi A., Frenk C. S., Ludlow A. D., Oman K. A., Schaller M., Theuns T., 2017, MNRAS, 471, 1841
  * Oman et al. (2020) Oman K. A., Brouwer M. M., Ludlow A. D., Navarro J. F., 2020, arXiv e-prints, p. arXiv:2006.06700
  * Paranjape & Sheth (2021) Paranjape A., Sheth R. K., 2021, MNRAS, 507, 632
  * Petersen & Lelli (2020) Petersen Lelli 2020, A&A, 636, A56
  * Phan et al. (2019) Phan D., Pradhan N., Jankowiak M., 2019, arXiv preprint arXiv:1912.11554
  * Pizagno et al. (2007) Pizagno J., et al., 2007, AJ, 134, 945
  * Pradyumna & Desai (2021) Pradyumna S., Desai S., 2021, Physics of the Dark Universe, 33, 100854
  * Sanders (1990) Sanders R. H., 1990, A&ARv, 2, 1
  * Tang et al. (2022) Tang et al., 2022, A&A, 668, A179
  * Tenneti et al. (2018) Tenneti A., Mao Y.-Y., Croft R. A. C., Di Matteo T., Kosowsky A., Zago F., Zentner A. R., 2018, MNRAS, 474, 3125
  * Tian et al. (2020) Tian Y., Umetsu K., Ko C.-M., Donahue M., Chiu I. N., 2020, ApJ, 896, 70
  * Tully & Fisher (1977) Tully R. B., Fisher J. R., 1977, A&A, 54, 661
  * Tully et al. (2023) Tully R. B., et al., 2023, ApJ, 944, 94
  * Zonoozi et al. (2021) Zonoozi A. H., Lieberz P., Banik I., Haghi H., Kroupa P., 2021, MNRAS, 506, 5468