11institutetext: Dipartimento di Fisica, Università di Torino, Via P. Giuria
1, I-10125 Torino, Italy 22institutetext: Istituto Nazionale di Fisica
Nucleare (INFN), Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy
33institutetext: I. Physikalisches Institut, Universität zu Köln , Zülpicher
Str. 77, D-50937 Köln, Germany

# Probing modified Newtonian dynamics
with hypervelocity stars

Sankha Subhra Chakrabarty<EMAIL_ADDRESS>Luisa Ostorero
1122 Arianna Gallo 1122 Stefano Ebagezio 3311 Antonaldo Diaferio 1122

(Received date; accepted date)

We show that measuring the velocity components of hypervelocity stars (HVSs)
can discriminate between modified Newtonian dynamics (MOND) and Newtonian
gravity. Hypervelocity stars are ejected from the Galactic center on radial
trajectories with a null tangential velocity component in the reference frame
of the Galaxy. They acquire tangential components due to the nonspherical
components of the Galactic gravitational potential. Axisymmetric potentials
only affect the latitudinal components, $v_{\theta}$, and non-null azimuthal
components, $v_{\phi}$, originate from non-axisymmetric matter distributions.
For HVSs with sufficiently high ejection speed, the azimuthal velocity
components are proportionate to the deviation of the gravitational potential
from axial symmetry. The ejection velocity threshold is $\sim$ 750 km s-1 for
4 $M_{\sun}$ stars and increases with decreasing HVS mass. We determine the
upper limit of $v_{\phi}$ as a function of the galactocentric distance for
these high-speed HVSs if MOND, in its quasi-linear formulation QUMOND, is the
correct theory of gravity and either the triaxial Galactic bulge or a
nonspherical hot gaseous halo is the primary source of the azimuthal
component, $v_{\phi}$. In Newtonian gravity, the HVSs within 60 kpc of the
Galactic center may easily have $v_{\phi}$ values higher than the QUMOND upper
limit if the dark matter halo is triaxial or if the dark matter halo and the
baryonic components are axisymmetric but their two axes of symmetry are
misaligned. Therefore, even a limited sample of high-speed HVSs could in
principle allow us to distinguish between the QUMOND scenario and the dark
matter model. This test is currently limited by (i) the lack of a proper
procedure to assess whether a star originates from the Galactic center and
thus is indeed an HVS in the model one wishes to constrain; and (ii) the large
uncertainties on the galactocentric azimuthal velocity components, which
should be reduced by at least a factor of $\sim 10$ to make this test
conclusive. A proper procedure to assess the HVS nature of the observed stars
and astrometric measurements with microarcsecond precision would make this
test feasible.

###### Key Words.:

Gravitations – Cosmology: dark matter – Galaxy: general – Galaxy: structure –
Galaxy: kinematics and dynamics

##  1 Introduction

The mass discrepancy in the Universe originates from a set of independent
observations on galactic and cosmological scales. On the scale of galaxies,
nearly flat rotation curves of disk galaxies at large radii (Bosma, 1978;
Rubin et al., 1982, 1985; Persic et al., 1996; McGaugh et al., 2000; Sofue &
Rubin, 2001; Martinsson et al., 2013; Bhattacharjee et al., 2014), the
stability of dynamically cold stellar disks (Hohl, 1971; Ostriker & Peebles,
1973; Fall & Efstathiou, 1980; Raha et al., 1991; Sellwood, 2014; Sellwood et
al., 2019), the dynamics of the outer regions of ellipticals (de Zeeuw &
Franx, 1991; Franx et al., 1991; Cappellari et al., 2006; Pota et al., 2013;
Cappellari, 2016; Pulsoni et al., 2018), and the large mass-to-light ratios of
dwarf galaxies (Kormendy, 1987; Irwin & Hatzidimitriou, 1995; Mateo, 1998; van
den Bosch & Swaters, 2001; Walker et al., 2009; Łokas, 2009; Simon, 2019)
cannot be explained in the standard theory of gravity without assuming the
existence of a large amount of dark matter (Blumenthal et al., 1984; White &
Frenk, 1991; Wechsler & Tinker, 2018; Zavala & Frenk, 2019). On the scale of
galaxy clusters (Biviano, 2000; Voit, 2005; Diaferio et al., 2008; Walker et
al., 2019), the dynamics of the member galaxies (Kneib et al., 1996; Giodini
et al., 2009; Geller et al., 2013; Rines et al., 2013; Sohn et al., 2017; Tian
et al., 2020), the X-ray emission of the intracluster gas (Sarazin, 1986;
Rosati et al., 2002; Böhringer & Werner, 2010; Sohn et al., 2019; Clerc et
al., 2020), and gravitational lensing (Tyson et al., 1990; Postman et al.,
2012; Sereno et al., 2018; Meneghetti et al., 2020; Umetsu, 2020) imply mass-
to-light ratios of order $\sim 100-400\leavevmode\nobreak\ M_{\sun}/L_{\sun}$
(Girardi et al., 2002; Rines et al., 2004; Proctor et al., 2015). On
cosmological scales, the structure formation from the nearly homogeneous
matter distribution implied by the small temperature anisotropies of the
cosmic microwave background (CMB) requires a stronger gravitational pull than
provided by the baryonic matter alone (Silk, 1967; Davis et al., 1985; Spergel
et al., 2007; Planck Collaboration et al., 2020).

The most popular and widely investigated solution to the problem of the mass
discrepancy is to assume the presence of cold dark matter (CDM) that is non-
baryonic and interacts with the baryons only via gravity (Peebles, 1982; Bond
et al., 1982; Blumenthal et al., 1982, 1984; Spergel et al., 2007; Frenk &
White, 2012; Strigari, 2013; Planck Collaboration et al., 2020). However, to
date, none of the elementary particles suggested as candidates of dark matter
has been detected (Tanabashi et al., 2018). The allowed windows for the
parameters associated with various dark matter particles are also shrinking
due to constraints from various terrestrial experiments and from astrophysical
observations (Bertone et al., 2005; de Martino et al., 2020).

In principle, the mass discrepancy problem can be solved with a modification
of the theory of gravity rather than with dark matter (Sanders, 1990; Sanders
& McGaugh, 2002; Clifton et al., 2012; Nojiri et al., 2017; de Martino et al.,
2020). Modified Newtonian dynamics (MOND; Milgrom, 1983c, a, b) is one of the
most investigated modifications of Newtonian gravity (Famaey & McGaugh, 2012;
McGaugh, 2020; de Martino et al., 2020). Modified Newtonian dynamics suggests
that Newtonian gravity breaks down in the low acceleration regime, where the
gravitational field is smaller than the critical value $a_{0}\approx 10^{-10}$
m s-2. Due to its dependence on the acceleration, the formulation of a
covariant version of MOND, which is required to test the theory against the
properties of the large-scale structure on cosmic scales, is not unique and
remains difficult (Bekenstein & Milgrom, 1984; Bekenstein, 2004; Milgrom,
2009, 2015; Hernandez et al., 2019; Skordis & Zlosnik, 2020).

On the other hand, on the scale of galaxies, MOND has proved to be predictive
and successful (e.g., Merritt, 2020): Confirmed MOND predictions include the
baryonic Tully-Fisher relation of disk galaxies with virtually no systematic
scatter over $\sim 5$ orders of magnitude in baryonic mass (McGaugh et al.,
2000; Lelli et al., 2016; McGaugh, 2020), the dependence of the shape of the
rotation curves on the surface brightness of the galaxy (McGaugh, 2020), and
the large mass-to-light ratios of dwarf galaxies when interpreted in Newtonian
gravity (Aaronson, 1983; Kormendy, 1987; Mateo et al., 1991; Mateo, 1998); as
expected in MOND, the strong equivalence principle also appears to be invalid,
as suggested by the analysis of accurate rotation curves of disk galaxies
(Chae et al., 2020).

A number of MOND predictions concerning the Milky Way (MW; Famaey & McGaugh,
2012) and the nearby dwarfs (Hodson et al., 2020) still need to be tested with
upcoming astrometric data. Indeed, when interpreted in Newtonian gravity, the
MOND gravitational field requires the presence of “phantom dark matter” in
addition to the baryonic matter (Milgrom, 1983a; Bekenstein & Milgrom, 1984).
Specifically, MOND predicts the existence of a disk of phantom dark matter in
the MW (Bienaymé et al., 2009); therefore, the acceleration of the stars
perpendicular to the plane of the disk can constrain the total surface density
of the baryonic and phantom disks (Nipoti et al., 2007). At 1.1 kpc above the
Galactic plane at the distance of the Sun from the Galactic center, the total
surface density in MOND is expected to be 60% larger than the surface density
of the baryonic disk, whereas this enhancement is predicted to be 51% in
Newtonian gravity with a spherical dark matter halo (Bienaymé et al., 2009;
Famaey & McGaugh, 2012). Similarly, the scale length of the “baryonic +
phantom” disk can be constrained by the vertical acceleration profiles at
different radii. In the formulation of MOND suggested by Bekenstein & Milgrom
(1984), this scale length should be 1.25 times the scale length of the visible
stellar disk (Bienaymé et al., 2009; Famaey & McGaugh, 2012). Moreover, the
angle between the minor axis of the velocity ellipsoid of the stars in the
solar neighborhood and the vertical direction depends on the shape of the
gravitational potential (Cuddeford & Amendt, 1991): For axisymmetric
potentials, the angles expected in MOND and Newtonian gravity with a spherical
dark matter halo differ by 2 degrees at the distance of the Sun from the MW
center and at 2 kpc above the MW disk (Bienaymé et al., 2009).

In this work we propose a novel prediction of MOND concerning the positions
and proper motions of hypervelocity stars (HVSs) within the MW. The existence
of HVSs was predicted by Hills (1988) as a result of three-body interactions
between the supermassive black hole (SMBH) at the Galactic center and binary
stars. After the interactions, one of the binary stars is ejected with high
radial velocity while the other is captured by the SMBH. Some other mechanisms
of generating HVSs have been proposed, including the interactions between a
star and the binary black hole at the center (Yu & Tremaine, 2003), the
interactions between binary stars and binary black holes (Wang et al., 2018),
the interactions between the MW and a dwarf galaxy (Abadi et al., 2009), and
star formation in the outflows driven by an active galactic nucleus (Wang &
Loeb, 2018; Silk et al., 2012). Since the serendipitous discovery of the first
HVS (Brown et al., 2005), $\sim 90$ high-velocity stars have been identified
as candidate HVSs (see, e.g., Hirsch et al., 2005; Edelmann et al., 2005;
Brown et al., 2006, 2009, 2012, 2014, 2018; Tillich et al., 2011; Li et al.,
2012, 2015; Pereira et al., 2013; Zheng et al., 2014; Huang et al., 2017;
Marchetti et al., 2017; Neugent et al., 2018; Massey et al., 2018; Hattori et
al., 2018a; Erkal et al., 2019; Du et al., 2019; Luna et al., 2019; Koposov et
al., 2020; Li et al., 2021).

The HVSs are ejected from the Galactic center with a speed of several hundred
km s-1 or higher. These stars are sometimes separated into samples of bound
and unbound stars. This distinction is irrelevant for our purpose here, and we
consider as an HVS any star that is ejected from the Galactic center on a
purely radial orbit with null tangential velocity in the galactocentric
reference frame. Hypervelocity stars obtain nonzero tangential velocities due
to the nonspherical components of the Galaxy gravitational potential. As they
travel to large distances from the center, the distribution of their positions
and velocities carries signatures of the gravitational field of the Galaxy.
Several schemes to constrain the shape of the halo (Gnedin et al., 2005; Rossi
et al., 2017; Contigiani et al., 2019) and to measure the virial mass of the
MW (Fragione & Loeb, 2017) using the HVSs have been proposed. Perets et al.
(2009) suggested a novel method to discriminate between various models of
Galactic potential within CDM and MOND paradigms by using the asymmetry of the
velocity distributions of incoming and outgoing HVSs. Hattori et al. (2018b)
also proposed a method for estimating the position and velocity of the Sun
within the MW by using the fact that the HVSs have low tangential velocities.
Kenyon et al. (2018) showed that tangential velocities may also be caused by
systems external to the MW, such as the Large Magellanic Cloud (LMC).

In this work we simulate the kinematics of the HVSs in MOND as well as in
Newtonian gravity. We show that the azimuthal components, $v_{\phi}$, of the
tangential velocities of the HVSs may distinguish MOND from Newtonian gravity.
Section 2 illustrates the quasi-linear formulation of MOND (QUMOND) that we
adopt in this work. Section 3 describes our model of the distribution of the
MW baryonic matter that generates the QUMOND gravitational potential. In Sect.
4 we illustrate the model of the dark matter halo we adopt for comparison with
the QUMOND predictions. In Sect. 5 we illustrate and discuss our simulations.
In Sect. 6 we show the galactocentric tangential velocities of the HVSs in
QUMOND and detail our QUMOND predictions. We conclude in Sect. 7.

##  2 Quasi-linear modified Newtonian dynamics

Throughout this work, we adopted QUMOND, the quasi-linear formulation of MOND
(Milgrom, 2010), where the gravitational field is

$\@vec{g}=\nu\left(\frac{|\@vec{g}_{\rm N}|}{a_{0}}\right)\ \@vec{g}_{\rm
N}\,$ (1)

and $\@vec{g}_{\rm N}$ is the Newtonian gravitational field due to the
baryonic matter alone. The interpolating function $\nu(x)$ satisfies the
limits $\nu(x)\to 1$ when $x\gg 1$, and $\nu(x)\to x^{-1/2}$ when $x\ll 1$. We
adopt

$\nu(x)=\left[\frac{1}{2}\left(1+\sqrt{1+4x^{-\gamma}}\right)\right]^{\frac{1}{\gamma}}\,,$
(2)

with $\gamma$ = 1 or $\gamma=$ 2 (Famaey & McGaugh, 2012). The function with
$\gamma$ = 1 is known as the simple interpolation function. The acceleration
scale, below which Newtonian gravity modifies, is set by $a_{0}$. The value of
$a_{0}$ is found by fitting the rotation curve data (Begeman et al., 1991;
Bottema et al., 2002), the observed correlation between the mass discrepancy
and the acceleration (McGaugh, 2004), and the baryonic Tully-Fisher relation
(McGaugh, 2011). The best-fit value of $a_{0}$ varies from 3000 to 4000 km2
s-2 kpc-1 and also slightly depends on the chosen interpolation function
$\nu(x)$ (Eq. 2). We chose an intermediate value: $a_{0}=3600$ km2 s-2
kpc${}^{-1}=1.2\times 10^{-10}$ m s-2.

For the Galaxy, we first assumed a simple model with three baryonic
components: a central SMBH, a stellar disk and a stellar bulge. With this
model, the Newtonian acceleration entering Eq. (1) is

$\@vec{g}_{\rm N}=-\@vec{\nabla}(\Phi_{\rm BH}+\Phi_{\rm Bulge}+\Phi_{\rm
Disk})\,,$ (3)

where $\Phi_{\rm BH}$, $\Phi_{\rm Bulge}$, and $\Phi_{\rm Disk}$ are the
Newtonian gravitational potentials that we provide in the next section. A more
sophisticated model that includes the additional baryonic component of a hot
gaseous (HG) halo will be discussed in Sect. 6.1.2.

##  3 Newtonian gravitational potentials

We investigated two variants of the model for the baryonic gravitational
potential that enters Eq. (3): an axisymmetric model and a non-axisymmetric
model. In the latter model, the deviation from the axial symmetry originates
only from the presence of a triaxial bulge. We used the reference frame of the
Galaxy with the origin at the Galaxy center. We used cylindrical coordinates
$(R,\phi,z)$ for the axisymmetric model and Cartesian coordinates $(x,y,z)$
for the triaxial model; $R$ lies on the $x$-$y$ plane, which we take as the
equatorial plane of the Galactic disk. For the spherically symmetric
components, we used spherical polar coordinates $(r,\theta,\phi)$.

###  3.1 The axisymmetric model

We adopted simple analytical potentials for the three MW components. For the
stellar disk, we used the Miyamoto-Nagai model (Miyamoto & Nagai, 1975)

$\Phi_{\rm Disk}(R,z)=-\frac{GM_{\rm D}}{\sqrt{R^{2}+\left(a_{\rm
D}+\sqrt{z^{2}+b_{\rm D}^{2}}\right)^{2}}}\,,$ (4)

with $M_{\rm D}=1.0\times 10^{11}\ M_{\sun}$, $a_{\rm D}=6.5\ {\rm kpc}$, and
$b_{\rm D}=0.26\ {\rm kpc}$ (Kafle et al., 2014; Price-Whelan et al., 2014;
Rossi et al., 2017; Contigiani et al., 2019).

For the bulge, we took the Hernquist sphere (Hernquist, 1990)

$\Phi_{\rm Bulge}(r)=-\frac{GM_{\rm B}}{r_{\rm B}+r}\,,$ (5)

with $M_{\rm B}=3.4\times 10^{10}\ M_{\sun}$, $r_{\rm B}=0.7\ {\rm kpc}$
(Kafle et al., 2014; Price-Whelan et al., 2014; Rossi et al., 2017; Contigiani
et al., 2019), and a total bulge mass of $3.0\times 10^{10}M_{\sun}$.

Finally, the gravitational potential of the central SMBH is

$\Phi_{\rm BH}(r)=-\frac{GM_{\rm BH}}{r}\,,$ (6)

with $M_{\rm BH}=4.0\times 10^{6}\ M_{\sun}$. This mass of the SMBH is
comparable to various estimates reported in the literature (e.g., Eisenhauer
et al., 2005; Ghez et al., 2008; Boehle et al., 2016; Gillessen et al., 2017).

In this model, the axial symmetry of the MW originates only from the stellar
disk because the gravitational potentials of the SMBH and the bulge are
spherically symmetric.

Figure 1:  Density profiles of the spherical (red) and triaxial (blue) bulges
as functions of the distance from the center. For the triaxial bulge, we show
the density profiles along the $x$, $y$, and $z$ axes. Figure 2:  Magnitude of
the radial acceleration in the plane of the disk, $|a_{R}|$, due to the bulge
alone in Newtonian gravity: spherical bulge (red) and triaxial bulge (blue).
For the triaxial bulge, the acceleration varies with the azimuthal angle,
$\phi$, which is the angle with respect to the $x$ axis in the plane of the
disk. We show the results for various values of $\phi$; $\phi=0^{\circ}$ and
$\phi=90^{\circ}$ correspond to the $x$ and $y$ axes, respectively.

###  3.2 The non-axisymmetric model: The triaxial bulge

To estimate the impact of a non-axisymmetric distribution of the baryonic
matter on the HVSs kinematics in QUMOND, we considered the effects of a
triaxial bulge and kept unaltered the gravitational potentials of the SMBH and
the stellar disk of the axisymmetric model. A triaxial bulge is the primary
source of azimuthal angular momentum within the Galaxy (Gardner et al., 2020).
Our adopted density profile of the triaxial bulge is (Binney et al., 1997;
McGaugh, 2008):

$\rho_{\rm Bulge}(x,y,z)=\frac{M_{0}}{\eta\zeta b_{\rm m}^{3}}\
\frac{e^{-(b/b_{\rm m})^{2}}}{\left(1+\frac{b}{b_{0}}\right)^{1.8}}\,,$ (7)

where

$b=\sqrt{x^{2}+\frac{y^{2}}{\eta^{2}}+\frac{z^{2}}{\zeta^{2}}}\,,$ (8)

$\eta=0.5$, $\zeta=0.6$, $b_{\rm m}=1.9$ kpc, and $b_{0}=0.1$ kpc (Binney et
al., 1997). The parameters are determined from the observed luminosity
distribution (Binney et al., 1997) and indicate that the largest axis of the
bulge, the $x$ axis, is in the plane of the disk, and the second largest axis,
the $z$ axis, is in the vertical direction. It follows that two principal axes
of the bulge are in the plane of the disk and, therefore, the bulge is not
tilted with respect to the plane of the disk. We chose the value of $M_{0}$ so
that the Newtonian accelerations due to the triaxial and the spherical bulges
are comparable at $R\gtrsim 0.5$ kpc in the plane of the disk; in other words,
the axisymmetric and non-axisymmetric models yield comparable rotation curves
beyond $\sim 0.5$ kpc and the dynamical differences between the two models are
limited to the central region. We set $M_{0}=5.7\times 10^{11}M_{\sun}$, which
yields a total mass of the bulge of $2.1\times 10^{10}M_{\sun}$. Figure 1
compares the density profiles of the spherical and triaxial bulges in our two
models.

We derived the Newtonian gravitational potential of the triaxial bulge,
$\Phi_{\rm Bulge}$, by solving the standard Poisson’s equation in three
dimensions within a box of volume (16 kpc)3 about the center. Figure 2 shows
that the magnitudes of the Newtonian radial acceleration for the triaxial
bulge are comparable to the Newtonian radial acceleration generated by the
spherical bulge of Eq. (5) at $R\gtrsim 0.5$ kpc, as we require with our
choice of $M_{0}$. Within a distance of 0.1 kpc from the center, the spherical
bulge generates a larger acceleration due to its steeper density profile shown
in Fig. 1. In the intermediate regions between 0.1 kpc and $\sim 0.5$ kpc, the
density profile of the spherical bulge is smaller than the density profile of
the triaxial bulge; however, the radial acceleration of the spherical bulge is
still larger than the acceleration of the triaxial case in any direction
because, due to the steeper central density, the spherical bulge encloses a
larger mass within $\sim 0.5$ kpc. At large distances, $r\gtrsim 5$ kpc, the
triaxiality is not effective and the mass distribution within the bulge can be
treated as spherically symmetric.

McGaugh (2008) investigates the mass model of the MW in the context of MOND by
matching the predicted rotation curve against observations. The scale length
of the stellar disk and the mass of the bulge determine the relative
contributions of these two components of the Galaxy to the total rotation
curve. Our model is roughly comparable to the model of McGaugh (2008) that has
the stellar disk with scale length $R_{\rm d}=$ 4 kpc and bulge mass within 1%
of the bulge mass in our model. Indeed, the circular velocity $v_{\rm circ}$
in the McGaugh model is within 5% of our $v_{\rm circ}$ within $r=3$ kpc, and
$\sim 10\%$ smaller at $r>3$ kpc. For different values of $R_{\rm d}$, the
agreement slightly worsens. For example, in McGaugh’s model with $R_{\rm
d}=2.3$ kpc, $v_{\rm circ}$ is $\sim 30\%$ smaller than our $v_{\rm circ}$
within $r=3$ kpc, but remains within $\sim 5\%$ at $r>3$ kpc.

Figure 3:  Schematic diagrams of our models of the Galaxy in Newtonian
gravity. The baryonic components, the bulge and the disk, are shown in black
and the dark matter halo in gray. Different shapes of the dark matter halo are
shown: (a) spherical (edge-on view of the disk), (b) prolate with the major
axis on the plane of disk (face-on view of the disk), (c) oblate with the
minor axis perpendicular to the disk (edge-on view of the disk), and (d)
prolate with the major axis perpendicular to the disk (edge-on view of the
disk). Only in case (b) are the axes of symmetry of the disk and the halo not
aligned.

##  4 The dark matter halo in Newtonian gravity

Figure 4:  Components of the velocity of a star in the spherical polar
coordinates with origin at the Galactic center (GC):
$\@vec{v}=\@vec{v}_{r}+\@vec{v}_{\theta}+\@vec{v}_{\phi}$.

We needed to compare the dynamics of HVSs in QUMOND with the corresponding
expectations in Newtonian gravity. We thus also considered a Newtonian model
of the MW with the same axisymmetric baryonic components described in the
previous section, but with an additional surrounding dark matter halo.

For the dark matter halo, we adopted the triaxial generalization suggested by
Vogelsberger et al. (2008) of the spherical Navarro-Frenk-White (NFW)
gravitational potential (Navarro et al., 1997):

$\Phi_{\rm Halo}=-\frac{GM_{\rm 200}}{f(C_{\rm 200})}\ \frac{1}{\tilde{r}}\
\ln\left(1+\frac{\tilde{r}}{r_{\rm s}}\right)\ ,$ (9)

where $f(u)=\ln(1+u)-u/(1+u)$. $M_{\rm 200}=8.35\times 10^{11}M_{\sun}$ is the
mass within $r_{\rm 200}$,111$r_{\rm 200}$ is the radius of a spherical volume
within which the mean mass density is 200 times the critical density of the
Universe. $C_{\rm 200}=r_{\rm 200}/r_{\rm s}=10.82$ is the concentration
parameter, and $r_{\rm s}=18$ kpc is the scale radius. Our adopted values of
the parameters are those used in Hesp & Helmi (2018). They are consistent with
the estimates from the kinematics of halo stars (Xue et al., 2008; Deason et
al., 2012). The generalized radius is

$\tilde{r}=\frac{r_{E}(r+r_{\rm a})}{r_{E}+r_{\rm a}}\,,$ (10)

where the ellipsoidal radius $r_{E}$ is

$r_{E}=\sqrt{\left(\frac{x}{a}\right)^{2}+\left(\frac{y}{b}\right)^{2}+\left(\frac{z}{c}\right)^{2}}\,$
(11)

and $r_{\rm a}=1.2r_{\rm s}$ (Hesp & Helmi, 2018) is the length scale that
determines the transition from the triaxial to the spherical shape: In the
inner region ($r\ll r_{\rm a}$) the halo is triaxial ($\tilde{r}\approx
r_{E}$), whereas in the outer region ($r\gg r_{\rm a}$) the halo is almost
spherical ($\tilde{r}\approx r$). This transition in the shape of the halo is
a generic prediction of $\Lambda$CDM simulations (see, e.g., Hayashi et al.
2007). The parameters $a$, $b$ and $c$ satisfy the relation
$a^{2}+b^{2}+c^{2}=3$.

We defined the two triaxiality parameters $q_{y}=b/a$ and $q_{z}=c/a$. Once
$q_{y}$ and $q_{z}$ are specified, $a$ is given by

$a=\sqrt{\frac{3}{1+q_{y}^{2}+q_{z}^{2}}}\ \ .$ (12)

Since the shape of the halo is currently poorly constrained (Bland-Hawthorn &
Gerhard, 2016), we varied either $q_{y}$ or $q_{z}$ within 40% from unity.

In this paper we did not use a triaxial dark matter halo. We only explored
spheroidal halos with different axes of symmetry (see Fig. 3). To study the
effects of the dark matter halo shape on the azimuthal component of the HVS
velocity, $v_{\phi}$ (see Fig. 4), we chose a spheroidal halo that is
axisymmetric about the $y$ axis lying in the plane of the disk (case (b) of
Fig. 3) whereas the baryonic components are axisymmetric about the vertical
$z$ axis. In this model, we set $q_{y}>1$ and $q_{z}=1$. To study the effects
of the dark matter halo shape on the latitudinal component $v_{\theta}$ of the
HVS velocities, we considered a spheroidal halo that is axisymmetric about the
$z$ axis (i.e., $q_{y}=1$ and $q_{z}\neq 1$). The halo is oblate if $q_{z}<1$
and prolate if $q_{z}>1$ (cases (c) and (d) of Fig. 3). Indeed, the kinematics
of halo stars from the Sloan Digital Sky Survey suggests that the dark matter
halo might be oblate within $\sim 20$ kpc, with $q_{z}=0.7\pm 0.1$, based on
the estimate of constant potential surfaces (Loebman et al., 2014).

At the solar neighborhood ($R=8$ kpc), our models yield a circular velocity of
235 km s-1 in Newtonian gravity and 240 (210) km s-1 in QUMOND with $\gamma=1$
($\gamma=2$). In addition, in our Newtonian model, the total mass enclosed
within 120 pc is in perfect agreement with the observed value reported in
Table 3 of Kenyon et al. (2008), derived from the estimate of the mass of the
central region of the stellar disk of Launhardt et al. (2002). The MW central
region is where the HVSs experience the largest deceleration. We verified that
varying the radial acceleration profile of the MW central region within the
observational uncertainties does not affect our conclusions on the tangential
velocities of the HVSs.

##  5 Simulations of the kinematics of the HVSs

In this section we describe our simulations and our synthetic samples of HVSs.
Section 5.1 illustrates the initial conditions of the equations of motion of
the HVSs that we adopted. In Sect. 5.2 we show that only stars beyond a
minimum galactocentric distance, corresponding to a minimum ejection velocity,
are relevant to discriminate between QUMOND and Newtonian gravity. This
threshold depends on the star mass: for $4\leavevmode\nobreak\ M_{\sun}$
stars, the minimum galactocentric distance is 15 kpc, which corresponds to the
minimum ejection velocity $\sim 750$ km s-1 for models in Newtonian gravity
and QUMOND with $\gamma=1$. For QUMOND with $\gamma=2$, the ejection velocity
threshold is $\sim 710$ km s-1, which also corresponds to a minimum distance
of 15 kpc.

###  5.1 Simulation setup

We assumed that the HVSs are generated through Hills’ mechanism (Hills, 1988):
In a three-body interaction between the SMBH associated with SgrA⋆ and a
binary star, one star of the binary is ejected and the other is captured by
the black hole. We emphasize that the assumption of Hills’ mechanism for the
ejection of the HVSs does not affect our conclusions, as detailed at the end
of this section.

In Hills’ ejection scenario, we simulated the velocity distribution of the
HVSs ejected from the Galactic center by means of a three-body numerical code
that reproduces the encounter of a set of equal-mass binary stars with the
SMBH. The binaries’ orbital parameters and minimum approach distance to the
SMBH are drawn from appropriate distributions (Bromley et al., 2006, and
references therein). More details on the numerical code and on the ejection
velocity distribution will be provided in Gallo et al. (2021) and in an
additional, separate paper. Here, we only report the information that is
instrumental to the present analysis.

We adopted a mass $M_{\rm BH}=4\times 10^{6}M_{\sun}$ for the black hole (as
in our model in Sect. 3.1) and a mass of 4 $M_{\sun}$ for each star in the
binary. The choice of mass for the binary stars is consistent with the fact
that most of the observed unbound HVSs have masses between 2.5 and 4
$M_{\sun}$ (Brown, 2015). In addition, choosing the upper limit of this mass
range is a conservative choice that is appropriate for our investigation, as
we clarify in Sect. 6. We obtained a distribution of ejection speeds that
displays a prominent peak at $v_{\rm ej}\sim 510$ km s-1, extends to
velocities of $\sim 4000$ km s-1, and has a positive skewness; 16% (12%) of
the ejected stars have speed $v_{\rm ej}\gtrsim 710$ (750) km s-1, a value
that will be relevant for the present analysis (see Sect. 5.2). Our results on
the ejection velocities are comparable to the results obtained from the
analytical prescriptions provided by Bromley et al. (2006).

The distribution of ejection speeds formally corresponds to the distribution
of velocities that the ejected stars would have at infinite distance from the
SMBH in absence of other sources of gravitational potential. However, in the
context of the Galactic center, this distribution can be taken as the velocity
distribution of the ejected stars at the starting point of their trajectories
across the MW, which we set as the radius of the sphere of influence of the
SMBH (i.e., $r=3$ pc; Genzel et al., 2010). Although the distribution of
ejection velocities weakly depends on the binary mass (Bromley et al., 2006),
our results regarding the kinematics of the HVSs are independent of their
mass, because the gravitational acceleration acting on the HVSs is mass-
independent.

Figure 5:  Distributions of HVSs in three sections of phase space in QUMOND
with $\gamma=1.0$ (blue dots) and Newtonian gravity with a prolate
($q_{y}=1,q_{z}=1.1$) dark matter halo (red dots). The left, middle, and right
panels show the $r$-$v_{r}$, $z$-$v_{z}$, and $r$-$v_{\theta}$ sections,
respectively. Here, $r$ is the galactocentric distance, $z$ is the vertical
coordinate, and $v_{r}$, $v_{z}$, $v_{\theta}$ are the radial, vertical, and
latitudinal components of the velocity, respectively.

The direction $\hat{n}(\theta,\phi)$ of the ejection velocity is assigned
randomly to each star, because Hills’ mechanism yields isotropic ejections. We
used the fourth order Runge-Kutta method with adaptive step-size to integrate
the equation of motion of each star. We started with a predefined time-step of
5 kyr and the time-step was adjusted so that, at each step, the position and
velocity of the star were determined with an accuracy of 1 pc and 0.01 km s-1,
respectively. For stars with ejection velocities higher than 600 km s-1, the
conservation of total energy holds with a relative accuracy of $10^{-12}$ for
each time-step as well as over the entire simulation. Stars with ejection
velocities lower than 600 km s-1 are not relevant for our study.

The typical lifetime of an isolated star of 4 $M_{\sun}$ on the main sequence
is $\tau_{\rm L}\simeq 160$ Myr, for Solar metallicities (Schaller et al.,
1992; Brown et al., 2006). We assumed this lifetime as the total lifetime of
our simulated stars. At the time of ejection, each star was assigned a random
age of $\tau_{\rm ej}=\epsilon\tau_{\rm L}$, where $\epsilon$ is a random
number drawn from a uniform density distribution between 0 and 1. We took the
average ejection rate to be $10^{-4}\ {\rm yr}^{-1}$, which is consistent with
the estimates by Yu & Tremaine (2003) and Zhang et al. (2013) (see also Hills,
1988). Therefore, in the simulation, we ejected stars in intervals of $\Delta
t=0.01$ Myr. For our chosen lifetime and ejection rate, the distribution of
ejected stars in the MW reaches a steady state after 160 Myr. We started the
simulation at $t=0$ and the $i$-th star was ejected at $t=(i-1)\Delta t$ with
its age $\tau^{i}_{\rm ej}$ at the time of ejection. We chose the time of
observation, which is the total run-time of the simulation, to be $t_{\rm
obs}=400$ Myr when the distribution of the HVSs is in the stationary regime.
For each star, we tested the condition of its survival: A star survives long
enough to be observed if $[t_{\rm obs}-(i-1)\Delta t]<(\tau_{\rm
L}-\tau^{i}_{\rm ej})$. If a star satisfies this condition, its travel time is
determined by $[t_{\rm obs}-(i-1)\Delta t]$.

We simulated the kinematics of the HVSs through the Galaxy (i) for QUMOND,
with the baryonic components only, and (ii) for Newtonian gravity, with the
baryonic components and the dark matter halo. For each model of the Galactic
potential, we determined positions and velocities of the HVSs at the time of
observation, $t=t_{\rm obs}$. Tables 1 and 2 summarize the simulations that we
performed in QUMOND and in Newtonian gravity. The tables list both the
simulations of our simpler model for the baryonic Galactic components
considered in Eq. (3) and described in Sect. 3 and the simulations that
include the additional nonspherical HG halo surrounding the MW (see Sect.
6.1.2).

Table 1: List of the simulations for QUMOND.

$\gamma$ | Bulge | Hot gaseous halo | Number of stars
---|---|---|---
1 | spherical | no | 40000
1 | triaxial | no | individual stars
2 | spherical | no | 40000
2 | triaxial | no | individual stars
1 | spherical | yes | 40000
2 | spherical | yes | 40000
222To find the upper limit of $|v_{\phi}|$ for QUMOND with a triaxial bulge,
we performed many simulations of individual stars. In this case, we used the
gravitational field of the triaxial bulge within 5 kpc of the Galactic center
and that of the spherical bulge at larger radii.

Table 2: List of the simulations for Newtonian gravity.

$q_{y}$ | $q_{z}$ | Panel of Fig. 3 | Hot gaseous halo
---|---|---|---
1 | 0.9 | c | no
1 | 1 | a | no
1 | 1.1 | d | no
1.1 | 1 | b | no
1.2 | 1 | b | no
1.2 | 1 | b | yes
1.4 | 1 | b | yes
333Simulations in Newtonian gravity were performed with a spherical bulge.
Each simulation contains 40000 stars.

We compared the final phase-space distributions of the HVSs obtained in QUMOND
with the distributions obtained in Newtonian gravity. Due to the large radial
velocities of the HVSs, their phase-space distributions in the $(r,v_{r})$ or
$(z,v_{z})$ space in QUMOND are virtually indistinguishable from those in
Newtonian gravity (left and middle panels of Fig. 5). However, the
distributions of the galactocentric tangential velocity components, the
latitudinal component $v_{\theta}=r({\mathrm{d}}\theta/{\mathrm{d}}t)$, and
the azimuthal component
$v_{\phi}=r\sin\theta({\mathrm{d}}\phi/{\mathrm{d}}t)$, are distinctive and
can be used, in principle, to discriminate between QUMOND and Newtonian
gravity (right panel of Fig. 5). Our result is independent of the mechanism
responsible for the ejection of the HVSs from the Galactic center: any
mechanism that can expel stars from the Galactic center with radial velocities
$v_{\rm ej}\gtrsim 710$ km s-1, as we discuss below, and null tangential
velocities would be suitable to perform the analysis presented in this work,
leading to the same conclusions.

###  5.2 Evolution of the tangential velocity components

The tangential velocity components of the HVSs are excellent probes of the
nonspherical components of the Galactic gravitational potential. However, this
distinctive ability only holds for HVSs with ejection speeds higher than a
threshold: for stars with lower ejection speed, the tangential velocity
components may indeed be disproportionately high.

To find the threshold ejection velocity, we simulated the dynamics of 200 HVS
with mass 4 $M_{\sun}$, ejection velocities $v_{\mathrm{ej}}$ between 650 and
850 km s-1, and random initial directions. For each star, the simulation time
was taken to be its lifetime, namely 160 Myr for 4 $M_{\sun}$ stars. This time
is the maximum possible travel time. For each star, we determined the maximum
values of the magnitudes of the tangential velocity components, $|v_{\theta}|$
and $|v_{\phi}|$.

We considered both Newtonian gravity with different shapes of the dark matter
halo and QUMOND with $\gamma=1$ and $\gamma=2$. For both the Newtonian and the
QUMOND scenarios, the baryonic components (Eqs. 4 \- 6) were taken to be
axisymmetric about the $z$ axis. Due to the pull of the disk, the stars obtain
nonzero $v_{\theta}$ values for all the models. However, in QUMOND with a
spherical bulge, the stars do not obtain any $|v_{\phi}|$ because the baryonic
matter, the only component present, is axisymmetric. In our models of
Newtonian gravity, the stars acquire nonzero values of $v_{\phi}$ only when
the two axes of symmetry of the baryonic and dark matter components are
misaligned: for example, for $(q_{y},q_{z})=(1.1,1)$ (panel (b) in Fig. 3),
the dark halo is symmetric about the $y$ axis whereas the baryonic matter is
symmetric about the $z$ axis. When no such misalignment is present, as in the
cases sketched in panels (a), (c), and (d) of Fig. 3, $v_{\phi}$ vanishes, and
the only non-null component of the tangential velocity is the latitudinal
velocity, $v_{\theta}$.

The left panel of Fig. 6 shows that $|v_{\theta}|$ can be as large as 600 km
s-1 for stars with $v_{\rm ej}<710$ km s-1 in QUMOND with $\gamma=2$, and for
stars with $v_{\rm ej}<750$ km s-1 in all the other models. In contrast, for
stars with ejection speed higher than these thresholds, the maximum values of
$|v_{\theta}|$ are consistently lower than $\sim 100$ km s-1. The right panel
of Fig. 6 shows qualitatively similar results for $|v_{\phi}|$ in the
Newtonian models with the misaligned axes of symmetry ($q_{y}\neq 1$): in
these models, $|v_{\phi}|$ is lower than $\sim 20$ km s-1 when $v_{\rm
ej}\gtrsim 750$ km s-1. In these Newtonian models with $q_{y}\neq 1$, the
vertical scatter in the maximum values of $|v_{\theta}|$ (left panel) is
large. The stars ejected with smaller angle with respect to the disk (i.e.,
$(90^{\circ}-\theta)\lesssim 30^{\circ}$) undergo larger deflection and attain
larger $|v_{\theta}|$. On the other hand, the stars ejected almost
perpendicular to the disk do not bend significantly, hence have smaller
$|v_{\theta}|$. Conversely, in the same Newtonian models, the stars always
attain a $|v_{\phi}|\gtrsim 350$ km s-1 (right panel) irrespective of the
direction of the HVS ejection velocity.

Figure 6:  Left panel: Maximum values of the latitudinal velocity components,
$|v_{\theta}|$, of 4 $M_{\sun}$ HVSs as a function of their ejection
velocities, $v_{\rm ej}$, in Newtonian gravity with different shapes of the
dark matter halo and in QUMOND with $\gamma=1$ and $\gamma=2$; symbols and
models are detailed in the inset. The vertical scatter in the plot originates
from the random directions of the star ejection. The vertical dotted lines
indicate $v_{\rm ej}=710$ and $v_{\rm ej}=750$ km s-1. In all the models
except QUMOND with $\gamma=2$, the maximum $|v_{\theta}|$ substantially drops
for stars with $v_{\rm ej}\gtrsim 750$ km s-1. For QUMOND with $\gamma=2$, the
maximum $|v_{\theta}|$ drops for stars with $v_{\rm ej}\gtrsim 710$ km s-1.
Right panel: Maximum values of the azimuthal velocity components,
$|v_{\phi}|$, of 4 $M_{\sun}$ HVSs as a function of their ejection velocities,
$v_{\rm ej}$, in two Newtonian gravity models where the axes of symmetry of
the baryonic matter and of the dark matter halo are misaligned. The other
models considered in the left panel do not appear because they are
axisymmetric and have zero $|v_{\phi}|$ values. The vertical scatter
originates from the random directions of the star ejection. The vertical
dotted line indicates $v_{\rm ej}=750$ km s-1.

High tangential velocities of the stars with lower ejection speed are caused
by the exchange of kinetic energy between radial and angular degrees of
freedom, especially when the stars undergo an inner turnaround. Figure 7
illustrates that the inner turnaround does occur for stars with low ejection
velocity. In Fig. 7 we show $v_{\phi}$ as a function of the radial component
$v_{r}$ for three stars with ejection velocity $v_{\rm ej}=\\{745,747,749\\}$
km s-1 in the Newtonian models with the misaligned axes, with
$(q_{y},q_{z})=(1.1,1)$. In each case, $v_{r}=v_{\rm ej}$ and $v_{\phi}=0$,
initially; at later times, $v_{\phi}$ starts increasing as $v_{r}$ decreases.
When $v_{r}=0$ for the first time, as indicated by the vertical dotted lines,
the star undergoes the outer turnaround, namely it reaches its maximum
distance from the MW center and starts moving inward. Therefore, $v_{r}$
becomes negative while $v_{\phi}$ keeps increasing. When $v_{r}$ becomes zero
for the second time, the star undergoes the inner turnaround, namely the star
reaches the closest approach to the center and starts again moving outward.
Around this phase of the inner turnaround, $|v_{\phi}|$ reaches its maximum,
as shown for $v_{\rm ej}=745$ km s-1. Stars with higher ejection speeds take
longer time to undergo the outer turnaround and may not live long enough to
experience the inner turnaround, as illustrated by the other two stars shown
in Fig. 7. As a result, the maximum of $|v_{\phi}|$ falls sharply around
$v_{\rm ej}\approx 750$ km s-1 (right panel of Fig. 6).

We conclude that for stars with ejection velocities higher than a threshold
velocity, the maximum values of the magnitudes of the latitudinal and
azimuthal velocity components, $|v_{\theta}|$ and $|v_{\phi}|$, are
proportionate to the departure from the spherical symmetry of the potential.
For stars with ejection speeds lower than the threshold, this proportionality
disappears because these stars may experience the inner turnaround. Our
simulations suggest that, for 4 $M_{\sun}$ stars, this threshold ejection
velocity is $\sim 710$ km s-1 for QUMOND with $\gamma=2$ and $\sim 750$ km s-1
for all the other models we investigate.

The ejection velocity is not observable but it is correlated with the outer
turnaround, namely the maximum distance from the MW center that the star can
travel. Figure 8 shows the outer turnaround as a function of the ejection
velocity of a star in different models of the Galactic potential. This figure
shows that neither in the Newtonian models nor in QUMOND with $\gamma=1$ can
stars with $v_{\rm ej}\lesssim 750$ km s-1 travel beyond a distance of about
15 kpc from the Galactic center. Similarly, in QUMOND with $\gamma=2$, the
stars with $v_{\rm ej}\lesssim 710$ km s-1 cannot travel beyond 15 kpc.
Therefore, identifying the $4\leavevmode\nobreak\ M_{\sun}$ stars that have
not experienced the inner turnaround, and thus have tangential velocities
proportionate to the deviation of the gravitational potential from axial
symmetry, requires considering only $4\leavevmode\nobreak\ M_{\sun}$ HVSs at
galactocentric distance $r>15$ kpc.

The threshold ejection velocity increases with decreasing HVS mass, because of
the longer lifetimes of lower-mass HVSs: longer-lived stars do not experience
the inner turnaround only if they travel to larger galactocentric distances;
in other words, only if they have larger ejection speeds. For example, for
$3\leavevmode\nobreak\ M_{\sun}$ or $2.5\leavevmode\nobreak\ M_{\sun}$ stars
with lifetimes on the main sequence $\sim 350$ Myr and $\sim 580$ Myr,
respectively (Schaller et al., 1992), the threshold ejection velocities are
$\sim 790$ km s-1 and $\sim 815$ km s-1 in Newtonian models; these threshold
ejection speeds correspond to outer turnaround radii of $\sim 30$ kpc and
$\sim 50$ kpc, respectively (Fig. 8). In QUMOND with $\gamma=1$ ($\gamma=2$),
for $3\leavevmode\nobreak\ M_{\sun}$ or $2.5\leavevmode\nobreak\ M_{\sun}$
stars, the threshold ejection velocity is $\sim 800$ ($755$) km s-1 and $\sim
830$ ($785$) km s-1, respectively; the corresponding turnaround radius is
$\sim 30$ ($30$) kpc and $\sim 50$ ($50$) kpc, respectively. Therefore, the
$3\leavevmode\nobreak\ M_{\sun}$ or $2.5\leavevmode\nobreak\ M_{\sun}$ stars
are of interest for our test only if they are at galactocentric distance $r$
larger than $30$ kpc and $50$ kpc, respectively.

Figure 7:  Azimuthal velocity component, $v_{\phi}$, as a function of the
radial velocity component, $v_{r}$, for three HVSs with mass 4 $M_{\sun}$ and
ejection velocity $v_{\rm ej}=\\{745,747,749\\}$ km s-1 in Newtonian gravity
with $(q_{y},q_{z})=(1.1,1)$. In each case, the star starts with $v_{r}=v_{\rm
ej}$ and $v_{\phi}=0$. The maxima of $v_{\phi}$ are quite different in the
three cases. The stars with $v_{\rm ej}=747$ and 749 km s-1 do not live long
enough to undergo the inner turnaround, that is, they do not encounter
$v_{r}=0$ for the second time. On the contrary, the inner turnaround occurs
for the star with $v_{\rm ej}=745$ km s-1. Figure 8:  HVS outer turnaround as
a function of the ejection speed, $v_{\rm ej}$, in Newtonian gravity with
different shapes of the dark matter halo and in QUMOND with $\gamma=1$ and
$\gamma=2$; models and symbols are detailed in the inset. The scatter of the
points originates from the random initial directions of the HVS ejection
velocities. The Newtonian gravitational potential is dominated by the dark
matter halo that becomes approximately spherical at large radii and makes the
scatter of the points decrease with increasing $v_{\rm ej}$. The vertical
dotted lines mark the thresholds $v_{\rm ej}=710$ km s-1 for QUMOND with
$\gamma=2$ and $v_{\rm ej}=750$ km s-1 for all the other models; the
horizontal line marks the minimum galactocentric distance of 15 kpc,
appropriate for 4 $M_{\sun}$ HVSs. Less massive stars have larger velocity
thresholds and larger minimum galactocentric distances.

##  6 Tangential velocity in QUMOND and Newtonian gravity

In this section we show the distributions of the azimuthal and latitudinal
components of the tangential velocity, $v_{\phi}$ and $v_{\theta}$
respectively, in both QUMOND and Newtonian gravity. We show that, for HVSs
within $\sim 60$ kpc, QUMOND yields upper limits for $|v_{\phi}|$ values that
are substantially lower than the values that non-axisymmetric gravitational
potentials can generate in Newtonian gravity. The QUMOND scenario might thus
be challenged if values of $|v_{\phi}|$ higher than those upper limits are
observed.

In QUMOND, the symmetry of the Galactic potential is determined by the
distribution of the baryonic components alone, unlike the case of Newtonian
gravity, where the shape of the dark matter halo plays a crucial and dominant
role. In QUMOND, if the baryonic distribution is axially symmetric, the pull
of the stellar disk generates non-null $v_{\theta}$ values, whereas $v_{\phi}$
is always zero. Non-null values of $v_{\phi}$ only appear in a non-
axisymmetric distribution of baryons, which mostly originates from a triaxial
bulge (Eq. 7) or from a nonspherical HG halo (Sect. 6.1.2).

In Newtonian gravity with a perfectly axisymmetric distribution of the
baryonic matter, as we adopt here (Eqs. 4-6), non-null values of $v_{\phi}$
can still appear if the halo of dark matter is not axially symmetric about the
same axis of symmetry of the baryonic distribution, which is the $z$ axis in
our models. More importantly, a tilted or triaxial dark matter halo acts on
the variation in $v_{\phi}$ for most of the HVS trajectory, whereas in QUMOND
the action of the triaxiality of the bulge is limited to the initial phase of
the HVS trajectory within $\sim 5$ kpc of the center.

This fundamental difference suggests that we might expect substantially higher
values of $v_{\phi}$ in Newtonian gravity than in QUMOND, unless in Newtonian
gravity either the dark matter halo is perfectly spherical or the halo is
axisymmetric and its axis of symmetry is aligned with the axis of symmetry of
the baryonic components. These possibilities, however, appear unlikely in the
current dark matter scenarios (Hayashi et al., 2007; Debattista et al., 2013).

A detailed study on how the shape of the dark matter halo in Newtonian gravity
affects the velocity components of the HVSs is presented elsewhere (Gallo et
al., 2021). Here, we briefly discuss the distributions of the azimuthal and
latitudinal velocity components $v_{\phi}$ and $v_{\theta}$ in QUMOND and
Newtonian gravity, and how $v_{\phi}$ can discriminate between the two
theories of gravity.

###  6.1 The upper limit of the azimuthal component, $v_{\phi}$, in QUMOND

In QUMOND, if the baryonic components are non-axisymmetric, the HVSs may have
non-null $v_{\phi}$ values. To estimate these values of $|v_{\phi}|$, in Sect.
6.1.1 we consider the effects of a triaxial bulge, which is the Galactic
baryonic component that is expected to play the major role in affecting
$|v_{\phi}|$ (Gardner et al., 2020). In Sect. 6.1.2, we illustrate the effects
of including a nonspherical halo of hot gas surrounding a MW with a spherical
bulge.

####  6.1.1 The role of a triaxial bulge

Finding the gravitational field of the triaxial bulge over the entire
numerical domain requires a demanding amount of computational time. We thus
adopted a simplified approach that uses the fact that the triaxiality of the
bulge is not effective at distances $r\gtrsim$ 5 kpc, according to the density
profile in Eq. (7) and the acceleration profile shown in Fig. 2. For our
purpose, at these large radii, the gravitational potential of the bulge is
well approximated by the potential generated by a spherically symmetric
source. Therefore, we only used the field of a triaxial bulge to simulate the
motion of the HVSs within 5 kpc of the Galactic center; at larger radii, we
adopted the gravitational field of a spherical bulge. We matched the two
regimes by setting the components of the velocity at $r=5$ kpc from the inner
region as the initial conditions for the outer region. The gravitational
acceleration at $r=5$ kpc due to the triaxial bulge is about 7% smaller than
the gravitational acceleration at $r=5$ kpc due to the spherical bulge,
suggesting that our approach introduces a limited error in the integration of
the equation of motion of the HVSs. In fact, this approximation does not
affect our main result: it produces a slight overestimate of the QUMOND upper
limit we discuss below, because the acceleration of the spherical bulge is
higher than the actual triaxal bulge acceleration.

As argued in Sect. 5.2, $4\leavevmode\nobreak\ M_{\sun}$ stars with ejection
velocity $v_{\rm ej}\lesssim 710$ km s-1 in QUMOND with $\gamma=2$, or with
$v_{\rm ej}\lesssim 750$ km s-1 in the other QUMOND and Newtonian models, tend
to have disproportionately high tangential velocity components and cannot be
used to probe the nonspherical components in the Galaxy potential. The maximum
value of $v_{\phi}$ is inversely proportional to the ejection velocity, as
expected by the argument illustrated in Sect. 5.2.

In order to find the maximum $|v_{\phi}|$ in QUMOND, we chose the lowest
ejection velocity relevant for our analysis. Figure 9 shows $|v_{\phi}|$ as a
function of the galactocentric distance $r$ for three different ejection
velocities: $v_{\rm ej}=710$ km s-1 in QUMOND with $\gamma=2$ (red curves),
$v_{\rm ej}=750$ km s-1 and 1000 km s-1 in QUMOND with $\gamma=1$ (red and
orange curves). All the stars shown in Fig. 9 are ejected at a polar angle
$\theta=45^{\circ}$; the solid lines show $v_{\phi}$ for an azimuthal angle
$\phi=45^{\circ}$, whereas the dashed and dot-dashed lines show $v_{\phi}$ for
the slowest stars for different $\phi$. We remind the reader that the largest
semimajor axis of the triaxial bulge is along the $x$ axis (see Eq. 7 and Fig.
4). According to Fig. 9, at $r=5$ kpc, the maximum $v_{\phi}$ we should expect
is $\sim 8$ km s-1. We note that this result is not sensitive to the choice of
the interpolating function $\nu$ (Eq. 2), because within 5 kpc the
gravitational field is mostly Newtonian.

Figure 9:  Magnitude of the azimuthal velocity component, $v_{\phi}$, in
QUMOND as a function of the galactocentric distance, $r$, for three 4
$M_{\sun}$ HVSs with different ejection velocities: $v_{\rm ej}=710$ km s-1 in
QUMOND with $\gamma=2$ (red), $v_{\rm ej}=750$ km s-1 in QUMOND with
$\gamma=1$ (blue), and $v_{\rm ej}=1000$ km s-1 in QUMOND with $\gamma=1$
(orange). The solid lines are for the stars ejected along $\theta=45^{\circ}$
and $\phi=45^{\circ}$, whereas the dashed (dash-dotted) line is for
$\theta=45^{\circ}$ and $\phi=15^{\circ}$ ($75^{\circ}$). In these examples,
the QUMOND field is due to the triaxial bulge and the axisymmetric disk. The
stars acquire $v_{\phi}$ only due to the triaxiality of the bulge.  Figure 10:
Magnitudes of the azimuthal velocity components $v_{\phi}$ of
$4\leavevmode\nobreak\ M_{\sun}$ HVSs in Newtonian gravity as a function of
their radial coordinates, $r$, at the time of observation for a prolate dark
matter halo with axis of symmetry in the plane of the disk. The light blue and
purple dots show the HVSs for a halo with $(q_{y},q_{z})=(1.1,1)$ and
$(q_{y},q_{z})=(1.2,1)$, respectively. The black line shows the upper limit of
$|v_{\phi}|$ in QUMOND due to the triaxial bulge. The upper limit also holds
for HVSs with masses lower than 4 $M_{\sun}$, provided they are beyond a
certain galactocentric distance, e.g., $r\gtrsim 30$ kpc (50 kpc) for
$3\leavevmode\nobreak\ M_{\sun}$ ($2.5\leavevmode\nobreak\ M_{\sun}$) HVSs.

Figure 11: Magnitudes of the azimuthal velocity components, $v_{\phi}$, of
$4\leavevmode\nobreak\ M_{\sun}$ HVSs as a function of their radial
coordinates, $r$, at the time of observation in the presence of an HG halo
with total mass of $1.5\times 10^{11}M_{\odot}$ within 100 kpc in both
Newtonian and QUMOND scenarios. The left and right panels show the results for
two different oblatenesses, $a_{\rm HG}=$ 0.4 and 0.8, of the HG halo (see
Eqs. 14-15).

Figure 12: Same as Fig. 11 except for the HG halo with a lower mass:
$1.5\times 10^{10}M_{\odot}$ within 100 kpc. Figure 13:  Distributions of the
latitudinal component of the velocity, $|v_{\theta}|$, in Newtonian gravity
and in QUMOND for 4 $M_{\sun}$ HVSs. The left panel shows the distributions of
$|v_{\theta}|$ for different models of the Galactic potential: Newtonian
gravity with a spherical dark matter halo (gray shaded histogram) and QUMOND
with $\gamma=1$ (orange histogram) and $\gamma=2$ (red histogram). As the
gravitational pull of the baryonic disk is enhanced in QUMOND, the fraction of
HVSs with high $|v_{\theta}|$ is larger than in Newtonian gravity with a
spherical halo. The right panel shows the distributions of $|v_{\theta}|$ for
three different shapes of the dark matter halo in Newtonian gravity. The three
halos have $q_{y}=1$ but different $q_{z}$: a spherical halo with $q_{z}=1$
(gray shaded histogram), a prolate halo with $q_{z}=1.1$ (blue histogram), and
an oblate halo with $q_{z}=0.9$ (red histogram). An oblate halo enhances the
gravitational pull of the baryonic disk, as shown by the larger fraction of
HVSs with high $|v_{\theta}|$ and the smaller fraction of HVSs with low
$|v_{\theta}|$ (compare the shaded and red histograms). The opposite occurs
with a prolate halo (compare the shaded and blue histograms).

Assuming that the angular momentum per unit mass
$|\ell_{z}|=r\sin\theta|v_{\phi}|$ is conserved at larger radii, and
$|v_{\phi}|$ thus falls off as $r^{-1}$, we obtain

$|v_{\phi}^{\rm max}(r)|\sim 4\ \left(\frac{10\ {\rm kpc}}{r}\right){\rm km\
s}^{-1}\,$ (13)

for 4 $M_{\sun}$ stars. The upper-limit radial profile of $|v_{\phi}|$
reported in Eq. (13) is shown by the black solid line in Fig. 10. The above
equation provides a conservative upper limit for two reasons. First, the
equation was derived using the 4 $M_{\sun}$ HVSs with the lowest ejection
velocities relevant for our analysis, namely 710 and 750 km s-1. Any HVS that
travels beyond 15 kpc (see Fig. 8) has higher ejection velocity, hence lower
azimuthal velocity. Secondly, the star trajectories bend toward the Galactic
disk located in the $x$-$y$ plane ($\theta=90^{\circ}$; see Fig. 4), so that
$\theta$ approaches $90^{\circ}$.444As the disk is located on
$\theta=90^{\circ}$, $\theta$ always approaches $90^{\circ}$. For an HVS
ejected above the disk, the initial value of $\theta$ is between $0^{\circ}$
and $90^{\circ}$, and $\theta$ increases. When the HVS is ejected below the
disk, the initial $\theta$ is between $90^{\circ}$ and $180^{\circ}$, and
$\theta$ decreases. As a result, $\sin\theta$ always increases and causes a
further decrease in $|v_{\phi}|$, because $|\ell_{z}|=r\sin\theta|v_{\phi}|$
remains constant. Therefore, using the HVS with the lowest ejection velocity
and taking $\sin\theta$ as a constant yield the conservative upper limit of
$|v_{\phi}|$ given in Eq. (13).

The radial profile $|v^{\rm max}_{\phi}(r)|$ of Eq. (13) also provides a
conservative upper limit for stars with masses lower than 4 $M_{\sun}$.
Indeed, as discussed in Sect. 5.2, stars with lower mass, and thus longer
lifetimes, must be ejected with higher speeds in order to avoid the inner
turnaround. Therefore, for higher ejection speeds, the normalization of the
azimuthal velocity radial profile, $|v_{\phi}(r)|$ in Eq. (13), will be lower.

In Newtonian gravity with an axisymmetric distribution of baryons (Eqs. 4-6),
the HVSs obtain non-null $|v_{\phi}|$ values only when the dark matter halo is
non-axisymmetric about the $z$ axis, as in our models with $q_{y}\neq 1$.
Figure 10 shows the values of $|v_{\phi}|$ and $r$ of the HVSs for two
different shapes of the halo: $(q_{y},q_{z})=(1.1,1)$, and
$(q_{y},q_{z})=(1.2,1)$. In both cases, the halo is prolate with the semimajor
axis along the $y$ axis in the plane of the disk. The halo becomes more
spherical at larger radii, $r\gtrsim r_{a}=21.6$ kpc. Thus, the azimuthal
angular momentum, $\ell_{z}=r\sin\theta\ v_{\phi}$, of each HVS becomes a
constant as the star travels beyond $r\sim r_{a}=21.6$ kpc and $|v_{\phi}|$
decreases as $r^{-1}$. Near $r\sim 21.6$ kpc, $|v_{\phi}|$ of an HVS can be as
high as 10 km s-1 for $q_{y}=1.2$, and 6 km s-1 for $q_{y}=1.1$ (see Fig. 10).
These values are substantially higher than the QUMOND upper limit
$v_{\phi}^{\rm max}(r)\sim 1.8$ km s-1 derived with Eq. (13) at that distance
and shown by the solid line in Fig. 10.

####  6.1.2 The effect of a hot gaseous halo

In addition to the baryonic components considered in the previous section, we
included in our MW model a reservoir of baryons in the form of an HG halo
extending up to the virial radius of the MW. The existence of this HG halo
with a temperature $\sim 10^{6}$ K is suggested by the O VII and O VIII
emission and absorption lines in the soft X-ray band (Paerels & Kahn, 2003;
Gupta et al., 2012; Fang et al., 2013; Gatto et al., 2013; Salem et al.,
2015). The presence of such diffuse hot gas may alleviate the missing baryon
problem 555The baryonic census in the Local Group appears to add up to only
$\sim 15\%$ (Fukugita et al., 1998) of the baryonic mass expected from the
estimated baryonic abundance ($\Omega_{b}h^{2}\approx 0.022$) (Planck
Collaboration et al., 2020) from the Big Bang nucleosynthesis and the CMB
anisotropies.  (Fukugita et al., 1998; Gupta et al., 2012; Shull et al., 2012;
Fang et al., 2013; Planck Collaboration et al., 2020). An oblate gaseous halo
with the smallest principal axis lying on the Galactic plane may be part of
the Vast Polar Structure of the MW (Pawlowski et al., 2011; Zhao et al., 2013;
Hammer et al., 2013).

We modeled the HG halo using the density profile (Thomas et al., 2017)

$\rho_{\rm HG}(m)=\rho_{\rm 0,HG}\left(1+\frac{m}{r_{\rm
0,HG}}\right)^{-3}\exp\left(-\frac{m^{2}}{r_{\rm t,HG}^{2}}\right)\,,$ (14)

where $r_{\rm 0,HG}=100$ kpc is the core radius and $r_{\rm t,HG}=200$ kpc is
the truncation radius. The elliptical radius $m$ for the oblate halo is

$m=\sqrt{\frac{x^{2}}{a_{\rm HG}^{2}}+y^{2}+z^{2}}\ .$ (15)

Currently, the shape of the halo has no observational constraints; we thus
varied $a_{\rm HG}$ from 0.4 to 0.8. Similarly, its total mass is uncertain by
a factor of ten (Miller & Bregman, 2013; Gatto et al., 2013; Salem et al.,
2015), and we thus explored two different values of $\rho_{\rm 0,HG}$:
$3.0\times 10^{5}$ and $3.0\times 10^{4}$ $M_{\odot}$ kpc-3. These values
yield a total mass of $1.5\times 10^{11}M_{\odot}$ and $1.5\times
10^{10}M_{\odot}$, respectively, within 100 kpc.

We explored the effects of the HG halo on the azimuthal velocities of the HVSs
in both the Newtonian and QUMOND scenarios (Figs. 11-12). In both cases, we
considered the axisymmetric models for the central black hole, bulge and disk
of Eqs. (4-6).

In Newtonian gravity, we considered the non-axisymmetric dark matter halos
with $(q_{y},q_{z})=(1.2,1)$ and $(q_{y},q_{z})=(1.4,1)$ and assumed that the
principal axes of the dark halo coincided with that of the HG halo. Whereas
the dark halo is nearly spherical at large radii ($r\gtrsim 21.6$ kpc) and its
contribution to the $|v_{\phi}|$ values thus falls as $r^{-1}$ at large
distances, the constant oblateness of the HG halo increases the $|v_{\phi}|$
values at all radii. Nevertheless the $|v_{\phi}|$ values are dominated by the
shape of the non-axisymmetric dark matter halo, because its mass is at least
10 times the mass of the HG halo.

On the contrary, in QUMOND, the gravitational field is substantially enhanced
at large distances by the presence of the HG halo and its mass and shape have
a relevant effect on the values of $|v_{\phi}|$. Therefore, for the HG halo
with the highest mass $M(<100\,\mathrm{kpc})=1.5\times 10^{11}\,M_{\odot}$
(Fig. 11), although the $|v_{\phi}|$’s in Newtonian scenario can still exceed
the values in QUMOND at distances between 15 kpc and 60 kpc irrespective of
the shape of the HG halo, at larger distances, the maximum possible values of
$|v_{\phi}|$ in QUMOND are higher than or comparable to the maximum possible
values in Newtonian gravity, depending on the shape of the HG halo. On the
contrary, for the HG halo with the lowest mass
$M(<100\,\mathrm{kpc})=1.5\times 10^{10}\,M_{\odot}$ (Fig. 12), the
$|v_{\phi}|$ values in Newtonian gravity can be significantly higher than the
values in QUMOND at distances between 15 kpc and 100 kpc.

In summary, the presence of an oblate HG halo increases the upper limit of the
$|v_{\phi}|$’s in QUMOND at all distances, compared to the presence of a
triaxial bulge alone. However, at smaller distances (15 kpc $\lesssim
r\lesssim$ 60 kpc), the $|v_{\phi}|$ values in Newtonian gravity due to a non-
axisymmetric dark matter halo may still substantially exceed the QUMOND
values.

###  6.2 Latitudinal component, $v_{\theta}$

The difference in the latitudinal component $v_{\theta}$ of the velocity in
QUMOND and in Newtonian gravity is subtler than the difference in $v_{\phi}$.
In QUMOND, the source of the gravitational field is concentrated in the plane
of the disk and it bends the HVSs trajectories toward this plane. In Newtonian
gravity, the main role is played by the spheroidal dark matter halo: A prolate
halo, with its major axis perpendicular to the stellar disk, is likely to
generate $v_{\theta}$ lower than in QUMOND; on the contrary, an oblate halo
with its major axis in the plane of the disk will generate $v_{\theta}$
comparable or even higher than in QUMOND. Accurate predictions of these
differences clearly depend on the exact values of the axial ratios, in
addition to the actual mass and size of the dark matter halo.

From our knowledge of the baryonic components, we can predict the distribution
of $|v_{\theta}|$ in Newtonian gravity, assuming that the dark matter halo is
spherical and therefore it does not affect $v_{\theta}$. Similarly, in our
adopted model, the SMBH and the bulge have spherically symmetric potential and
do not contribute to $v_{\theta}$. The shaded histogram in the left panel of
Fig. 13 shows the distribution of $|v_{\theta}|$ in this Newtonian model,
effectively caused by the disk alone. In QUMOND, at larger distances, where
the Newtonian gravitational field approaches $a_{0}$, the gravitational pull
of the disk gets enhanced compared to Newtonian gravity. Consequently, the
HVSs have larger $|v_{\theta}|$ in QUMOND than in Newtonian gravity.

Figure 14 shows the magnitude of $v_{\theta}$ as a function of the
galactocentric distance $r$ due to the baryonic components in QUMOND with
$\gamma=1$ and 2 (orange and red dots, respectively), and in Newtonian gravity
with a spherical halo (black dots). In both models, the maximum possible
values of $|v_{\theta}|$ decrease as $\sim r^{-1}$, because, as each star
travels beyond the length scale of the stellar disk, its angular momentum,
$\sim rv_{\theta}$, becomes a constant. Because of the enhanced pull of
QUMOND, the maximum value of $v_{\theta}$ is higher than in Newtonian gravity
at all radii. Unfortunately, unlike the case of $v_{\phi}$, these differences
cannot be used to distinguish QUMOND from Newtonian gravity, because these
same differences can be generated by an appropriate shape of the dark matter
halo, as we illustrate below.

In Newtonian gravity with a dark matter halo, non-null $v_{\theta}$ values are
caused by the Galactic disk and by the halo with $q_{z}\neq 1$ (models
sketched in panels (c) and (d) of Fig. 3). For an oblate halo ($q_{z}<1$), the
gravitational pull of the halo enhances the pull of the disk. Hence, the
number of HVSs with high $|v_{\theta}|$ is larger than in the case of a
spherical halo, as shown by the comparison of the red and shaded histograms in
the right panel of Fig. 13. The reverse happens for a prolate halo with
$q_{z}>1$. Therefore, the difference between the QUMOND and the Newtonian
distributions shown in the left panel of Fig. 13 can be easily mimicked in
Newtonian gravity by an appropriate oblate dark matter halo.

This degeneracy was emphasized by Read & Moore (2005) when they attempted to
distinguish MOND from Newtonian gravity using the stellar streams of the
Sagittarius dwarf. The degeneracy between an oblate halo and QUMOND can be
broken if $q_{z}$ is sufficiently small. Indeed, similar to the case of
$v_{\phi}$, QUMOND sets an upper limit to $|v_{\theta}|$; on the contrary, in
Newtonian gravity, $|v_{\theta}|$ may be higher than this upper limit if
$q_{z}$ is sufficiently small. In our models, the Newtonian $|v_{\theta}|$
values are higher than the QUMOND upper limit if $q_{z}\lesssim 0.6$. However,
the difference between the QUMOND and Newtonian $|v_{\theta}|$ values is not
as prominent as for $|v_{\phi}|$ values: In our model, the transition scale of
the halo shape from oblate to spherical is $r_{\rm a}=21.6$ kpc; therefore,
the size of the portion of the halo that is actually oblate is comparable to
the size of the baryonic disk. On the contrary, the size of the triaxial bulge
causing nonzero $|v_{\phi}|$ values in QUMOND is much smaller than the scale
length of the non-axisymmetric halo responsible for non-null $|v_{\phi}|$’s in
Newtonian gravity. Therefore, $v_{\theta}$ is less effective than $v_{\phi}$
at distinguishing QUMOND from Newtonian gravity. Figure 15 quantifies this
difficulty by showing $|v_{\theta}|$ versus $r$ for QUMOND with $\gamma=1$ and
$\gamma=2$, and for Newtonian gravity with an oblate halo with $q_{z}=0.4$;
this dark matter halo is a rather extreme case, when compared to the current
estimates of the shape of the MW dark matter halo (Loebman et al., 2014).

Figure 14:  Latitudinal velocity components, $|v_{\theta}|$, as a function of
the galactocentric distance, $r$, of 4 $M_{\sun}$ HVSs in QUMOND with
$\gamma=1$ (orange dots) and $\gamma=2$ (red dots), and in Newtonian gravity
with a spherical dark matter halo ($q_{z}=1$) (black dots). Figure 15:
Latitudinal velocity components, $|v_{\theta}|$, as a function of the
galactocentric distance, $r$, of 4 $M_{\sun}$ HVSs in QUMOND with $\gamma=1$
(orange dots) and $\gamma=2$ (red dots), and in Newtonian gravity with an
oblate dark matter halo with $q_{z}=0.4$ (dark cyan dots).

###  6.3 Azimuthal velocities: A comparison with real data

We attempted a first comparison of the QUMOND upper limit $|v_{\phi}^{\rm
max}(r)|$ (Eq. 13) with the measured $|v_{\phi}^{\rm obs}|$ of nine stars
drawn from the HVS survey sample of Brown et al. (2014). These nine stars have
masses in the range $\sim 2.5-4\leavevmode\nobreak\ M_{\sun}$. Their
galactocentric distances are larger than the minimum mass-dependent radii
required for our test (see Sect. 5.2). In addition, Kenyon et al. (2018) show
that HVSs ejected within $25^{\circ}$ from the line joining the MW center and
the LMC are affected by the LMC when they reach a galactocentric distance of
$35$ to $65$ kpc. The nine HVSs of our sample are all located more than $\sim
84^{\circ}$ away from the LMC, and the LMC pull should thus be irrelevant.

To be used for our test, these stars must originate from the Galactic center
rather than being disk runaway stars. In principle, we could distinguish the
two kinds of stars by tracing their trajectories back in time, if we knew the
correct theory of gravity and the correct MW gravitational potential. When
this information is unknown, and it is actually what we wish to constrain,
this approach clearly generates a circularity problem. We can solve this
problem by tracing back the star trajectories in different gravitational
potentials and different theories of gravity to select those stars, if any,
that appear to be HVSs in all models. We will investigate this self-consistent
treatment, and thus the actual feasibility of our test, elsewhere. Here, we
simply wished to see whether, if we assumed that the observed allegedly HVSs
were indeed HVSs in all the models, the current data would have sufficed to
distinguish between QUMOND and Newtonian gravity.

We relied on the analysis of Brown et al. (2018) who, in Newtonian gravity,
assume an axisymmetric disk and a spherical dark matter halo (Kenyon et al.,
2014) to estimate a larger probability for the nine stars mentioned above to
come from the Galactic center than to be disk runaway stars.

We computed the azimuthal components, $v_{\phi}^{\rm obs}$, of their
galactocentric velocities from the proper motions available in the Gaia Early
Data Release 3 (EDR3; Gaia Collaboration et al., 2016, 2020). We adopted the
heliocentric distances derived from Brown et al. (2018), and the radial
velocities of Brown et al. (2014). We found that the magnitudes of the
azimuthal velocities are in the range $|v_{\phi}^{\rm obs}|\simeq 37-452$ km
s-1. For eight of these stars, the uncertainties on the azimuthal velocities
are dominated by the errors on the proper motions, whereas the errors on the
stars’ distances and radial velocities, as well as on the position and
velocity of the Sun in the galactocentric reference frame, appear to be
negligible. For these stars, the relative uncertainties on the azimuthal
velocities are in the range $\sim 50-340$%. For the ninth star, B329,
neglecting the error on its distance is inappropriate. B329 is a $(3.21\pm
0.24)\leavevmode\nobreak\ M_{\sun}$ star at galactocentric distance $r=(61\pm
13)$ kpc; by ignoring the uncertainty on the distance, we derived
$|v_{\phi}^{\rm obs}|=(371\pm 74)$ km s-1, whereas including the uncertainty
on the distance yields $|v_{\phi}^{\rm obs}|=(371\pm 131)$ km s-1.

We conclude that the current uncertainties on $|v_{\phi}^{\rm obs}|$ make the
azimuthal components consistent with zero within 3$\sigma$ for all of the nine
stars of our sample. We are thus unable to verify whether these measures would
be in principle consistent with the QUMOND limit. Indeed, even in Newtonian
gravity, the large values of $|v_{\phi}^{\rm obs}|$ would probably imply an
unrealistically large flatness or triaxiality of the dark matter halo (Fig.
10). Therefore, even if these nine stars were HVSs in both QUMOND and
Newtonian gravity, the comparison of their azimuthal components with the
QUMOND upper limit, $|v_{\phi}^{\rm max}(r)|$, would still be inconclusive.

##  7 Discussion and conclusions

We showed that measuring the galactocentric tangential velocity of HVSs in the
MW can effectively allow us to discriminate between MOND, in its QUMOND
formulation, and Newtonian gravity. Specifically, we demonstrated that HVSs
with sufficiently high ejection speed possess galactocentric azimuthal
velocities whose magnitude, $|v_{\phi}|$, cannot exceed a velocity threshold
in QUMOND, while $|v_{\phi}|$ has no upper bounds in Newtonian gravity. This
result could naturally translate into an observational test to discriminate
between these two theories of gravity.

Our findings follow from the fact that the HVSs ejected from the Galaxy center
on radial orbits acquire non-null azimuthal tangential speeds, $v_{\phi}$, due
to the non-axisymmetric components of the Galactic gravitational potential.
Hypervelocity stars with low ejection velocity can turn back toward the
Galactic center and move outward again, acquiring a high $|v_{\phi}|$ that is
not proportionate to the deviation from the axial symmetry of the potential.
In contrast, HVSs with high ejection velocities reach larger galactocentric
distances before turning back and thus die before experiencing the inner
turnaround. These HVSs acquire a substantially lower $|v_{\phi}|$, which is
proportional to the deviation from the axial symmetry of the potential.

In our models of the Galactic gravitational potential, the ejection velocity
threshold for $4\leavevmode\nobreak\ M_{\sun}$ stars is $\sim 710$ km s-1 in
QUMOND with $\gamma=2$, and it is $\sim 750$ km s-1 both in QUMOND with
$\gamma=1$ and in the Newtonian gravity models that we investigated here.
Hypervelocity stars with ejection velocities higher than this value reach
galactocentric distances larger than $\sim 15$ kpc if they live long enough.
We thus expect that the $|v_{\phi}|$ component of $4\leavevmode\nobreak\
M_{\sun}$ HVSs beyond this distance can be used to probe the deviation from
the axial symmetry of the Galactic potential.

The ejection velocity threshold and the corresponding minimum galactocentric
distance increase with decreasing HVS mass. For example, in Newtonian gravity,
$3\leavevmode\nobreak\ M_{\sun}$ and $2.5\leavevmode\nobreak\ M_{\sun}$ HVSs
have ejection velocity thresholds of $\sim 790$ km s-1 and $\sim 815$ km s-1,
which correspond to the minimum galactocentric distances of 30 kpc and 50 kpc,
respectively. In QUMOND with $\gamma=1$ ($\gamma=2$), $3\leavevmode\nobreak\
M_{\sun}$ and $2.5\leavevmode\nobreak\ M_{\sun}$ HVSs have ejection velocity
thresholds of $\sim 800$ ($755$) km s-1 and $\sim 830$ ($785$) km s-1, which
correspond to the minimum galactocentric distances of 30 ($30$) kpc and 50
($50$) kpc, respectively (see Sect. 5.2). Therefore, the mass of the star sets
the minimum distances beyond which the star must be located if we wish to use
it for the test we suggested here.

In Newtonian gravity, the symmetric features of the baryonic distribution can
be overcome by those of the dark matter halo surrounding the Galaxy. Here, we
explored the case where both the baryonic components and the dark matter halo
are axisymmetric but their axes are misaligned; the shape of the halo depends
on the distance from the Galactic center, and the halo approaches the
spherical symmetry at distances larger than $\sim 21.6$ kpc. It follows that
$v_{\phi}$ is affected by the non-axisymmetric features of the potential up to
this radius.

Among the baryonic components, only a triaxial bulge as well as the possible
presence of a nonspherical HG halo can affect $v_{\phi}$. Additional
variations in $v_{\phi}$ might derive from massive external objects, such as
the LMC (Kenyon et al., 2018), whereas other non-axisymmetric components
within the disk, such as spiral arms or density inhomogeneities, are expected
to affect only the stars within $\sim 1$ kpc of the disk (Gardner et al.,
2020). The contribution from the HG halo in $|v_{\phi}|$ strongly depends on
its shape and total mass, which are poorly constrained by observations.

In our QUMOND model with a triaxial bulge, the triaxiality of the bulge is
effective up to $r\sim 5$ kpc; beyond this radius, the bulge can be
approximated to be spherical and, because of the conservation of angular
momentum, the maximum value of $|v_{\phi}|$ is proportional to $r^{-1}$ (Eq.
(13)). Therefore, when compared with the Newtonian model affecting
$|v_{\phi}|$ out to $\sim 21.6$ kpc, the values of $|v_{\phi}|$ of HVSs at
$r\gtrsim 15$ kpc in QUMOND are substantially lower than the values allowed
for a non-axisymmetric potential in Newtonian gravity. For example, at $r\sim
20$ kpc, we find $|v_{\phi}|\lesssim$ 2 km s-1 in QUMOND, whereas $|v_{\phi}|$
can be as high as 6 or 10 km s-1 for a dark matter halo with axial ratios
$q_{y}=1.1$ or 1.2, respectively (see Fig. 10).

If a nonspherical HG halo is included in the MW model, it dominates the
$|v_{\phi}|$ values of the HVSs in QUMOND, and the $|v_{\phi}|$ values in this
case may be higher than those generated by the triaxial bulge alone. On the
contrary, in Newtonian gravity, the $|v_{\phi}|$ values are still dominated by
the dark matter halo and can still largely exceed the QUMOND values for HVSs
at distances of up to 60 kpc from the center.

We conclude that precise measurements of $v_{\phi}$ for a few
$4\leavevmode\nobreak\ M_{\sun}$ HVSs at galactocentric distances larger than
$\sim 15$ kpc (or $3\leavevmode\nobreak\ M_{\sun}$ HVSs at $r\gtrsim 30$ kpc,
or $2.5\leavevmode\nobreak\ M_{\sun}$ HVSs at $r\gtrsim 50$ kpc) and smaller
than $\sim 60$ kpc may in principle test the validity of QUMOND. Finding a few
HVSs with azimuthal components, $|v_{\phi}|$, above the QUMOND upper limits,
$|v_{\phi}^{\rm max}|$, given by Eq. (13) (black line in Fig. 10) or shown in
Figs. 11-12 in the presence of the HG halo, would suggest that MOND may not be
the correct theory of gravity, at least in its QUMOND formulation. Such tests
clearly require that the HVSs are confirmed to originate from the Galactic
center and that their trajectories are not perturbed by external objects, such
as the LMC.

Assessing the HVS nature of observed stars is far from trivial, however,
because we need to trace their trajectories back in time with the theory of
gravity and in the gravitational potential that we wish to constrain. In
addition, in QUMOND the perturbation of the trajectories by objects beyond the
MW is complicated by the external field effect that is absent in Newtonian
gravity (e.g., Haghi et al., 2016; Famaey et al., 2018; Hodson et al., 2020).
The external field effect is expected to be more relevant for the most distant
stars and needs to be quantified. We plan to tackle these issues elsewhere.

If we assume that the stars currently identified as HVSs are indeed HVSs in
both QUMOND and Newtonian gravity, we conclude that the current uncertainties
on their azimuthal velocities, $v_{\phi}$, mostly due to the large relative
uncertainties on the proper motion measurements, are too large to provide a
conclusive comparison of the data with our QUMOND limit. Indeed, in the
subsample of nine HVSs drawn from the sample of Brown et al. (2014) that we
considered here, the azimuthal velocities have relative uncertainties in the
range $35$-$340\%$. The precision on the estimates of $v_{\phi}$ thus needs to
be improved by at least a factor of $\sim 10$ to make our test decisive.
Future measurements from space-borne astrometric missions with expected
microarcsecond precision on star positions, such as Theia (The Theia
Collaboration et al., 2017; Malbet et al., 2019, 2021), are expected to allow
us to discriminate between the two theories of gravity.

###### Acknowledgements.

We are sincerely grateful to the referee who provided a number of insightful
and constructive suggestions that helped to clarify parts of our work that
were initially vague. We thank Warren Brown, Margaret Geller and Scott Kenyon
for stimulating discussions on the topic of hypervelocity stars, and for
constructive comments on the manuscript. SSC was supported by the grant “The
Milky Way and Dwarf Weights with Space Scales” funded by University of Torino
and Compagnia di S. Paolo (UniTO-CSP), by the grant no. IDROL 70541 IDRF
2020.0756 funded by Fondazione CRT, by INFN, and by the Departments of
Excellence grant L.232/2016 of the Italian Ministry of Education, University
and Research (MIUR). This last grant fully supported the PhD fellowship of AG.
We acknowledge partial support from the INFN grant InDark. The work of SE
included here was part of his Master Thesis project at the University of
Torino. This work has made use of data from the European Space Agency (ESA)
mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data
Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This research has made use
of NASA’s Astrophysics Data System Bibliographic Services. This work has made
use of Topcat (Taylor, 2005) and the following Python modules: Matplotlib
(Hunter, 2007), NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020),
AstroPy (Astropy Collaboration et al., 2013, 2018) GalPy (Bovy, 2015).

## References

  * Aaronson (1983) Aaronson, M. 1983, ApJ, 266, L11
  * Abadi et al. (2009) Abadi, M. G., Navarro, J. F., & Steinmetz, M. 2009, ApJ, 691, L63
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * 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
  * Bekenstein (2004) Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509
  * Bertone et al. (2005) Bertone, G., Hooper, D., & Silk, J. 2005, Phys. Rep, 405, 279
  * Bhattacharjee et al. (2014) Bhattacharjee, P., Chaudhury, S., & Kundu, S. 2014, ApJ, 785, 63
  * Bienaymé et al. (2009) Bienaymé, O., Famaey, B., Wu, X., Zhao, H. S., & Aubert, D. 2009, A&A, 500, 801
  * Binney et al. (1997) Binney, J., Gerhard, O., & Spergel, D. 1997, MNRAS, 288, 365
  * Biviano (2000) Biviano, A. 2000, in Constructing the Universe with Clusters of Galaxies, ed. F. Durret & D. Gerbal, 1
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J. & Gerhard, O. 2016, ARA&A, 54, 529
  * Blumenthal et al. (1984) Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Nature, 311, 517
  * Blumenthal et al. (1982) Blumenthal, G. R., Pagels, H., & Primack, J. R. 1982, Nature, 299, 37
  * Boehle et al. (2016) Boehle, A., Ghez, A. M., Schödel, R., et al. 2016, ApJ, 830, 17
  * Böhringer & Werner (2010) Böhringer, H. & Werner, N. 2010, A&A Rev., 18, 127
  * Bond et al. (1982) Bond, J. R., Silk, J., & Kolb, E. W. 1982, ApJ, 255, 341
  * Bosma (1978) Bosma, A. 1978, PhD thesis, -
  * Bottema et al. (2002) Bottema, R., Pestaña, J. L. G., Rothberg, B., & Sand ers, R. H. 2002, A&A, 393, 453
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29
  * Bromley et al. (2006) Bromley, B. C., Kenyon, S. J., Geller, M. J., et al. 2006, ApJ, 653, 1194
  * Brown (2015) Brown, W. R. 2015, ARA&A, 53, 15
  * Brown et al. (2009) Brown, W. R., Geller, M. J., & Kenyon, S. J. 2009, ApJ, 690, 1639
  * Brown et al. (2012) Brown, W. R., Geller, M. J., & Kenyon, S. J. 2012, ApJ, 751, 55
  * Brown et al. (2014) Brown, W. R., Geller, M. J., & Kenyon, S. J. 2014, ApJ, 787, 89
  * Brown et al. (2005) Brown, W. R., Geller, M. J., Kenyon, S. J., & Kurtz, M. J. 2005, ApJ, 622, L33
  * Brown et al. (2006) Brown, W. R., Geller, M. J., Kenyon, S. J., & Kurtz, M. J. 2006, ApJ, 647, 303
  * Brown et al. (2018) Brown, W. R., Lattanzi, M. G., Kenyon, S. J., & Geller, M. J. 2018, ApJ, 866, 39
  * Cappellari (2016) Cappellari, M. 2016, ARA&A, 54, 597
  * Cappellari et al. (2006) Cappellari, M., Bacon, R., Bureau, M., et al. 2006, MNRAS, 366, 1126
  * Chae et al. (2020) Chae, K.-H., Lelli, F., Desmond, H., et al. 2020, ApJ, 904, 51
  * Clerc et al. (2020) Clerc, N., Kirkpatrick, C. C., Finoguenov, A., et al. 2020, MNRAS, 497, 3976
  * Clifton et al. (2012) Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep, 513, 1
  * Contigiani et al. (2019) Contigiani, O., Rossi, E. M., & Marchetti, T. 2019, Mon. Not. Roy. Astron. Soc., 487, 4025
  * Cuddeford & Amendt (1991) Cuddeford, P. & Amendt, P. 1991, MNRAS, 253, 427
  * Davis et al. (1985) Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371
  * de Martino et al. (2020) de Martino, I., Chakrabarty, S. S., Cesare, V., et al. 2020, Universe, 6, 107
  * de Zeeuw & Franx (1991) de Zeeuw, T. & Franx, M. 1991, Annual Review of Astronomy and Astrophysics, 29, 239
  * Deason et al. (2012) Deason, A. J., Belokurov, V., Evans, N. W., et al. 2012, MNRAS, 425, 2840
  * Debattista et al. (2013) Debattista, V. P., Roškar, R., Valluri, M., et al. 2013, MNRAS, 434, 2971
  * Diaferio et al. (2008) Diaferio, A., Schindler, S., & Dolag, K. 2008, Space Sci. Rev., 134, 7
  * Du et al. (2019) Du, C., Li, H., Yan, Y., et al. 2019, ApJS, 244, 4
  * Edelmann et al. (2005) Edelmann, H., Napiwotzki, R., Heber, U., Christlieb, N., & Reimers, D. 2005, ApJ, 634, L181
  * Eisenhauer et al. (2005) Eisenhauer, F., Genzel, R., Alexander, T., et al. 2005, ApJ, 628, 246
  * Erkal et al. (2019) Erkal, D., Boubert, D., Gualandris, A., Evans, N. W., & Antonini, F. 2019, MNRAS, 483, 2007
  * Fall & Efstathiou (1980) Fall, S. M. & Efstathiou, G. 1980, MNRAS, 193, 189
  * Famaey et al. (2018) Famaey, B., McGaugh, S., & Milgrom, M. 2018, MNRAS, 480, 473
  * Famaey & McGaugh (2012) Famaey, B. & McGaugh, S. S. 2012, Living Reviews in Relativity, 15, 10
  * Fang et al. (2013) Fang, T., Bullock, J., & Boylan-Kolchin, M. 2013, ApJ, 762, 20
  * Fragione & Loeb (2017) Fragione, G. & Loeb, A. 2017, New Astronomy, 55, 32–38
  * Franx et al. (1991) Franx, M., Illingworth, G., & de Zeeuw, T. 1991, ApJ, 383, 112
  * Frenk & White (2012) Frenk, C. S. & White, S. D. M. 2012, Annalen der Physik, 524, 507
  * Fukugita et al. (1998) Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Luri, X., Chemin, L., et al. 2020, arXiv e-prints, arXiv:2012.01771
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gallo et al. (2021) Gallo, A., Ostorero, L., Chakrabarty, S. S., Ebagezio, S., & Diaferio, A. 2021, A&A, submitted, arXiv:2111.09657
  * Gardner et al. (2020) Gardner, S., Hinkel, A., & Yanny, B. 2020, Astrophys. J., 890, 110
  * Gatto et al. (2013) Gatto, A., Fraternali, F., Read, J. I., et al. 2013, MNRAS, 433, 2749
  * Geller et al. (2013) Geller, M. J., Diaferio, A., Rines, K. J., & Serra, A. L. 2013, ApJ, 764, 58
  * Genzel et al. (2010) Genzel, R., Eisenhauer, F., & Gillessen, S. 2010, Reviews of Modern Physics, 82, 3121
  * Ghez et al. (2008) Ghez, A. M., Salim, S., Weinberg, N. N., et al. 2008, ApJ, 689, 1044
  * Gillessen et al. (2017) Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30
  * Giodini et al. (2009) Giodini, S., Pierini, D., Finoguenov, A., et al. 2009, ApJ, 703, 982
  * Girardi et al. (2002) Girardi, M., Manzato, P., Mezzetti, M., Giuricin, G., & Limboz, F. 2002, ApJ, 569, 720
  * Gnedin et al. (2005) Gnedin, O. Y., Gould, A., Miralda-Escudé, J., & Zentner, A. R. 2005, ApJ, 634, 344
  * Gupta et al. (2012) Gupta, A., Mathur, S., Krongold, Y., Nicastro, F., & Galeazzi, M. 2012, ApJ, 756, L8
  * Haghi et al. (2016) Haghi, H., Bazkiaei, A. E., Zonoozi, A. H., & Kroupa, P. 2016, MNRAS, 458, 4172
  * Hammer et al. (2013) Hammer, F., Yang, Y., Fouquet, S., et al. 2013, MNRAS, 431, 3543
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hattori et al. (2018a) Hattori, K., Valluri, M., Bell, E. F., & Roederer, I. U. 2018a, ApJ, 866, 121
  * Hattori et al. (2018b) Hattori, K., Valluri, M., & Castro, N. 2018b, ApJ, 869, 33
  * Hayashi et al. (2007) Hayashi, E., Navarro, J. F., & Springel, V. 2007, MNRAS, 377, 50
  * Hernandez et al. (2019) Hernandez, X., Sussman, R. A., & Nasser, L. 2019, MNRAS, 483, 147
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
  * Hesp & Helmi (2018) Hesp, C. & Helmi, A. 2018, arXiv e-prints, arXiv:1804.03670
  * Hills (1988) Hills, J. G. 1988, Nature, 331, 687
  * Hirsch et al. (2005) Hirsch, H. A., Heber, U., O’Toole, S. J., & Bresolin, F. 2005, A&A, 444, L61
  * Hodson et al. (2020) Hodson, A. O., Diaferio, A., & Ostorero, L. 2020, A&A, 640, A26
  * Hohl (1971) Hohl, F. 1971, ApJ, 168, 343
  * Huang et al. (2017) Huang, Y., Liu, X. W., Zhang, H. W., et al. 2017, ApJ, 847, L9
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Irwin & Hatzidimitriou (1995) Irwin, M. & Hatzidimitriou, D. 1995, MNRAS, 277, 1354
  * Kafle et al. (2014) Kafle, P. R., Sharma, S., Lewis, G. F., & Bland-Hawthorn, J. 2014, ApJ, 794, 59
  * Kenyon et al. (2014) Kenyon, S. J., Bromley, B. C., Brown, W. R., & Geller, M. J. 2014, ApJ, 793, 122
  * Kenyon et al. (2018) Kenyon, S. J., Bromley, B. C., Brown, W. R., & Geller, M. J. 2018, ApJ, 864, 130
  * Kenyon et al. (2008) Kenyon, S. J., Bromley, B. C., Geller, M. J., & Brown, W. R. 2008, ApJ, 680, 312
  * Kneib et al. (1996) Kneib, J. P., Ellis, R. S., Smail, I., Couch, W. J., & Sharples, R. M. 1996, ApJ, 471, 643
  * Koposov et al. (2020) Koposov, S. E., Boubert, D., Li, T. S., et al. 2020, MNRAS, 491, 2465
  * Kormendy (1987) Kormendy, J. 1987, in Dark matter in the universe, ed. J. Kormendy & G. R. Knapp, Vol. 117, 139
  * Launhardt et al. (2002) Launhardt, R., Zylka, R., & Mezger, P. G. 2002, A&A, 384, 112
  * Lelli et al. (2016) Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, ApJ, 816, L14
  * Li et al. (2012) Li, Y., Luo, A., Zhao, G., et al. 2012, ApJ, 744, L24
  * Li et al. (2021) Li, Y.-B., Luo, A. L., Lu, Y.-J., et al. 2021, ApJS, 252, 3
  * Li et al. (2015) Li, Y.-B., Luo, A. L., Zhao, G., et al. 2015, Research in Astronomy and Astrophysics, 15, 1364
  * Loebman et al. (2014) Loebman, S. R., Ivezić, Ž., Quinn, T. R., et al. 2014, The Astrophysical Journal, 794, 151
  * Loebman et al. (2014) Loebman, S. R., Ivezić, Ž., Quinn, T. R., et al. 2014, ApJ, 794, 151
  * Łokas (2009) Łokas, E. L. 2009, MNRAS, 394, L102
  * Luna et al. (2019) Luna, A., Minniti, D., & Alonso-García, J. 2019, ApJ, 887, L39
  * Malbet et al. (2019) Malbet, F., Abbas, U., Alves, J., et al. 2019, arXiv e-prints, arXiv:1910.08028
  * Malbet et al. (2021) Malbet, F., Boehm, C., Krone-Martins, A., et al. 2021, Experimental Astronomy
  * Marchetti et al. (2017) Marchetti, T., Rossi, E. M., Kordopatis, G., et al. 2017, MNRAS, 470, 1388
  * Martinsson et al. (2013) Martinsson, T. P. K., Verheijen, M. A. W., Westfall, K. B., et al. 2013, A&A, 557, A131
  * Massey et al. (2018) Massey, P., Levine, S. E., Neugent, K. F., et al. 2018, AJ, 156, 265
  * Mateo et al. (1991) Mateo, M., Olszewski, E., Welch, D. L., Fischer, P., & Kunkel, W. 1991, AJ, 102, 914
  * Mateo (1998) Mateo, M. L. 1998, ARA&A, 36, 435
  * McGaugh (2020) McGaugh, S. 2020, Galaxies, 8, 35
  * McGaugh (2004) McGaugh, S. S. 2004, ApJ, 609, 652
  * McGaugh (2008) McGaugh, S. S. 2008, ApJ, 683, 137
  * McGaugh (2011) McGaugh, S. S. 2011, Phys. Rev. Lett., 106, 121303
  * McGaugh et al. (2000) McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. 2000, Astrophys. J. Lett., 533, L99
  * Meneghetti et al. (2020) Meneghetti, M., Davoli, G., Bergamini, P., et al. 2020, Science, 369, 1347
  * Merritt (2020) Merritt, D. 2020, A Philosophical Approach to MOND: Assessing the Milgromian Research Program in Cosmology
  * Milgrom (1983a) Milgrom, M. 1983a, ApJ, 270, 371
  * Milgrom (1983b) Milgrom, M. 1983b, ApJ, 270, 384
  * Milgrom (1983c) Milgrom, M. 1983c, ApJ, 270, 365
  * Milgrom (2009) Milgrom, M. 2009, Phys. Rev. D, 80, 123536
  * Milgrom (2010) Milgrom, M. 2010, MNRAS, 403, 886
  * Milgrom (2015) Milgrom, M. 2015, Canadian Journal of Physics, 93, 107
  * Miller & Bregman (2013) Miller, M. J. & Bregman, J. N. 2013, ApJ, 770, 118
  * Miyamoto & Nagai (1975) Miyamoto, M. & Nagai, R. 1975, PASJ, 27, 533
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Neugent et al. (2018) Neugent, K. F., Massey, P., Morrell, N. I., Skiff, B., & Georgy, C. 2018, AJ, 155, 207
  * Nipoti et al. (2007) Nipoti, C., Londrillo, P., Zhao, H., & Ciotti, L. 2007, MNRAS, 379, 597
  * Nojiri et al. (2017) Nojiri, S., Odintsov, S. D., & Oikonomou, V. K. 2017, Phys. Rep, 692, 1
  * Ostriker & Peebles (1973) Ostriker, J. P. & Peebles, P. J. E. 1973, ApJ, 186, 467
  * Paerels & Kahn (2003) Paerels, F. B. S. & Kahn, S. M. 2003, ARA&A, 41, 291
  * Pawlowski et al. (2011) Pawlowski, M. S., Kroupa, P., & de Boer, K. S. 2011, A&A, 532, A118
  * Peebles (1982) Peebles, P. J. E. 1982, ApJ, 263, L1
  * Pereira et al. (2013) Pereira, C. B., Jilinski, E. G., Drake, N. A., Ortega, V. G., & Roig, F. 2013, A&A, 559, A12
  * Perets et al. (2009) Perets, H. B., Wu, X., Zhao, H. S., et al. 2009, ApJ, 697, 2096
  * Persic et al. (1996) Persic, M., Salucci, P., & Stel, F. 1996, MNRAS, 281, 27
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Postman et al. (2012) Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25
  * Pota et al. (2013) Pota, V., Forbes, D. A., Romanowsky, A. J., et al. 2013, MNRAS, 428, 389
  * Price-Whelan et al. (2014) Price-Whelan, A. M., Hogg, D. W., Johnston, K. V., & Hendel, D. 2014, ApJ, 794, 4
  * Proctor et al. (2015) Proctor, R. N., Mendes de Oliveira, C., Azanha, L., Dupke, R., & Overzier, R. 2015, MNRAS, 449, 2345
  * Pulsoni et al. (2018) Pulsoni, C., Gerhard, O., Arnaboldi, M., et al. 2018, A&A, 618, A94
  * Raha et al. (1991) Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, Nature, 352, 411
  * Read & Moore (2005) Read, J. I. & Moore, B. 2005, MNRAS, 361, 971
  * Rines et al. (2013) Rines, K., Geller, M. J., Diaferio, A., & Kurtz, M. J. 2013, ApJ, 767, 15
  * Rines et al. (2004) Rines, K., Geller, M. J., Diaferio, A., Kurtz, M. J., & Jarrett, T. H. 2004, AJ, 128, 1078
  * Rosati et al. (2002) Rosati, P., Borgani, S., & Norman, C. 2002, ARA&A, 40, 539
  * Rossi et al. (2017) Rossi, E. M., Marchetti, T., Cacciato, M., Kuiack, M., & Sari, R. 2017, MNRAS, 467, 1844
  * Rubin et al. (1985) Rubin, V. C., Burstein, D., Ford, W. K., J., & Thonnard, N. 1985, ApJ, 289, 81
  * Rubin et al. (1982) Rubin, V. C., Ford, W. K., J., Thonnard, N., & Burstein, D. 1982, ApJ, 261, 439
  * Salem et al. (2015) Salem, M., Besla, G., Bryan, G., et al. 2015, ApJ, 815, 77
  * Sanders (1990) Sanders, R. H. 1990, A&A Rev., 2, 1
  * Sanders & McGaugh (2002) Sanders, R. H. & McGaugh, S. S. 2002, ARA&A, 40, 263
  * Sarazin (1986) Sarazin, C. L. 1986, Reviews of Modern Physics, 58, 1
  * Schaller et al. (1992) Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269
  * Sellwood (2014) Sellwood, J. A. 2014, Reviews of Modern Physics, 86, 1
  * Sellwood et al. (2019) Sellwood, J. A., Shen, J., & Li, Z. 2019, MNRAS, 486, 4710
  * Sereno et al. (2018) Sereno, M., Giocoli, C., Izzo, L., et al. 2018, Nature Astronomy, 2, 744
  * Shull et al. (2012) Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23
  * Silk (1967) Silk, J. 1967, Nature, 215, 1155
  * Silk et al. (2012) Silk, J., Antonuccio-Delogu, V., Dubois, Y., et al. 2012, A&A, 545, L11
  * Simon (2019) Simon, J. D. 2019, ARA&A, 57, 375
  * Skordis & Zlosnik (2020) Skordis, C. & Zlosnik, T. 2020, arXiv e-prints, arXiv:2007.00082
  * Sofue & Rubin (2001) Sofue, Y. & Rubin, V. 2001, ARA&A, 39, 137
  * Sohn et al. (2019) Sohn, J., Geller, M. J., & Zahid, H. J. 2019, ApJ, 880, 142
  * Sohn et al. (2017) Sohn, J., Geller, M. J., Zahid, H. J., et al. 2017, ApJS, 229, 20
  * Spergel et al. (2007) Spergel, D. N., Bean, R., Doré, O., et al. 2007, ApJS, 170, 377
  * Strigari (2013) Strigari, L. E. 2013, Phys. Rep, 531, 1
  * Tanabashi et al. (2018) Tanabashi, M., Hagiwara, K., Hikasa, K., et al. 2018, Phys. Rev. D, 98, 030001
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * The Theia Collaboration et al. (2017) The Theia Collaboration, Boehm, C., Krone-Martins, A., et al. 2017, arXiv e-prints, arXiv:1707.01348
  * Thomas et al. (2017) Thomas, G. F., Famaey, B., Ibata, R., Lüghausen, F., & Kroupa, P. 2017, A&A, 603, A65
  * Tian et al. (2020) Tian, Y., Yu, P.-C., Li, P., McGaugh, S. S., & Ko, C.-M. 2020, arXiv e-prints, arXiv:2010.00992
  * Tillich et al. (2011) Tillich, A., Heber, U., Geier, S., et al. 2011, A&A, 527, A137
  * Tyson et al. (1990) Tyson, J. A., Valdes, F., & Wenk, R. A. 1990, ApJ, 349, L1
  * Umetsu (2020) Umetsu, K. 2020, A&A Rev., 28, 7
  * van den Bosch & Swaters (2001) van den Bosch, F. C. & Swaters, R. A. 2001, MNRAS, 325, 1017
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2008) Vogelsberger, M., White, S. D. M., Helmi, A., & Springel, V. 2008, MNRAS, 385, 236
  * Voit (2005) Voit, G. M. 2005, Reviews of Modern Physics, 77, 207
  * Walker et al. (2009) Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2009, ApJ, 704, 1274
  * Walker et al. (2019) Walker, S., Simionescu, A., Nagai, D., et al. 2019, Space Sci. Rev., 215, 7
  * Wang & Loeb (2018) Wang, X. & Loeb, A. 2018, New A, 61, 95
  * Wang et al. (2018) Wang, Y.-H., Leigh, N., Yuan, Y.-F., & Perna, R. 2018, MNRAS, 475, 4595
  * Wechsler & Tinker (2018) Wechsler, R. H. & Tinker, J. L. 2018, ARA&A, 56, 435
  * White & Frenk (1991) White, S. D. M. & Frenk, C. S. 1991, ApJ, 379, 52
  * Xue et al. (2008) Xue, X. X., Rix, H. W., Zhao, G., et al. 2008, ApJ, 684, 1143
  * Yu & Tremaine (2003) Yu, Q. & Tremaine, S. 2003, ApJ, 599, 1129
  * Zavala & Frenk (2019) Zavala, J. & Frenk, C. S. 2019, Galaxies, 7, 81
  * Zhang et al. (2013) Zhang, F., Lu, Y., & Yu, Q. 2013, ApJ, 768, 153
  * Zhao et al. (2013) Zhao, H., Famaey, B., Lüghausen, F., & Kroupa, P. 2013, A&A, 557, L3
  * Zheng et al. (2014) Zheng, Z., Carlin, J. L., Beers, T. C., et al. 2014, ApJ, 785, L23