Issue 
A&A
Volume 657, January 2022



Article Number  A115  
Number of page(s)  16  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/202141136  
Published online  20 January 2022 
Probing modified Newtonian dynamics with hypervelocity stars
^{1}
Dipartimento di Fisica, Università di Torino, Via P. Giuria 1, 10125 Torino, Italy
email: sankhasubhra.chakrabarty@unito.it
^{2}
Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy
^{3}
I. Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, 50937 Köln, Germany
Received:
20
April
2021
Accepted:
1
October
2021
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_{θ}, and nonnull azimuthal components, v_{ϕ}, originate from nonaxisymmetric 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 ∼750 km s^{−1} for 4 M_{⊙} stars and increases with decreasing HVS mass. We determine the upper limit of v_{ϕ} as a function of the galactocentric distance for these highspeed HVSs if MOND, in its quasilinear 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_{ϕ}. In Newtonian gravity, the HVSs within 60 kpc of the Galactic center may easily have v_{ϕ} 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 highspeed 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 ∼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: gravitation / dark matter / Galaxy: general / Galaxy: structure / Galaxy: kinematics and dynamics
© ESO 2022
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 masstolight 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. 2021), the Xray 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 masstolight ratios of order ∼100 − 400 M_{⊙}/L_{⊙} (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 VI 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 nonbaryonic 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 VI 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 1983a,b,c) 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} ≈ 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 largescale structure on cosmic scales, is not unique and remains difficult (Bekenstein & Milgrom 1984; Bekenstein 2004; Milgrom 2009, 2015; Hernandez et al. 2019; Skordis & Zlosnik 2021).
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 ∼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 masstolight 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 1983b; 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 threebody 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), ∼90 highvelocity 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, 2021; 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).
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_{ϕ}, of the tangential velocities of the HVSs may distinguish MOND from Newtonian gravity. Section 2 illustrates the quasilinear 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. Quasilinear modified Newtonian dynamics
Throughout this work, we adopted QUMOND, the quasilinear formulation of MOND (Milgrom 2010), where the gravitational field is
and g_{N} is the Newtonian gravitational field due to the baryonic matter alone. The interpolating function ν(x) satisfies the limits ν(x)→1 when x ≫ 1, and ν(x)→x^{−1/2} when x ≪ 1. We adopt
with γ = 1 or γ = 2 (Famaey & McGaugh 2012). The function with γ = 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 bestfit value of a_{0} varies from 3000 to 4000 km^{2} s^{−2} kpc^{−1} and also slightly depends on the chosen interpolation function ν(x) (Eq. (2)). We chose an intermediate value: a_{0} = 3600 km^{2} s^{−2} kpc^{−1} = 1.2 × 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
where Φ_{BH}, Φ_{Bulge}, and Φ_{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 nonaxisymmetric 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, ϕ, 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, θ, ϕ).
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)
with M_{D} = 1.0 × 10^{11} M_{⊙}, a_{D} = 6.5 kpc, and b_{D} = 0.26 kpc (Kafle et al. 2014; PriceWhelan et al. 2014; Rossi et al. 2017; Contigiani et al. 2019).
For the bulge, we took the Hernquist sphere (Hernquist 1990)
with M_{B} = 3.4 × 10^{10} M_{⊙}, r_{B} = 0.7 kpc (Kafle et al. 2014; PriceWhelan et al. 2014; Rossi et al. 2017; Contigiani et al. 2019), and a total bulge mass of 3.0 × 10^{10} M_{⊙}.
Finally, the gravitational potential of the central SMBH is
with M_{BH} = 4.0 × 10^{6} M_{⊙}. 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.
3.2. The nonaxisymmetric model: The triaxial bulge
To estimate the impact of a nonaxisymmetric 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):
where
η = 0.5, ζ = 0.6, b_{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 ≳ 0.5 kpc in the plane of the disk; in other words, the axisymmetric and nonaxisymmetric models yield comparable rotation curves beyond ∼0.5 kpc and the dynamical differences between the two models are limited to the central region. We set M_{0} = 5.7 × 10^{11} M_{⊙}, which yields a total mass of the bulge of 2.1 × 10^{10} M_{⊙}. Figure 1 compares the density profiles of the spherical and triaxial bulges in our two models.
Fig. 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. 
We derived the Newtonian gravitational potential of the triaxial bulge, Φ_{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 ≳ 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 ∼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 ∼0.5 kpc. At large distances, r ≳ 5 kpc, the triaxiality is not effective and the mass distribution within the bulge can be treated as spherically symmetric.
Fig. 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, ϕ, which is the angle with respect to the x axis in the plane of the disk. We show the results for various values of ϕ; ϕ = 0° and ϕ = 90° correspond to the x and y axes, respectively. 
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_{d} = 4 kpc and bulge mass within 1% of the bulge mass in our model. Indeed, the circular velocity v_{circ} in the McGaugh model is within 5% of our v_{circ} within r = 3 kpc, and ∼10% smaller at r > 3 kpc. For different values of R_{d}, the agreement slightly worsens. For example, in McGaugh’s model with R_{d} = 2.3 kpc, v_{circ} is ∼30% smaller than our v_{circ} within r = 3 kpc, but remains within ∼5% at r > 3 kpc.
4. The dark matter halo in Newtonian gravity
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 NavarroFrenkWhite (NFW) gravitational potential (Navarro et al. 1997):
where f(u) = ln(1 + u)−u/(1 + u). M_{200} = 8.35 × 10^{11} M_{⊙} is the mass within r_{200}^{1}, C_{200} = r_{200}/r_{s} = 10.82 is the concentration parameter, and r_{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
where the ellipsoidal radius r_{E} is
and r_{a} = 1.2 r_{s} (Hesp & Helmi 2018) is the length scale that determines the transition from the triaxial to the spherical shape: In the inner region (r ≪ r_{a}) the halo is triaxial (), whereas in the outer region (r ≫ r_{a}) the halo is almost spherical (). This transition in the shape of the halo is a generic prediction of Λ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
Since the shape of the halo is currently poorly constrained (BlandHawthorn & 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_{ϕ} (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_{θ} of the HVS velocities, we considered a spheroidal halo that is axisymmetric about the z axis (i.e., q_{y} = 1 and q_{z} ≠ 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 ∼20 kpc, with q_{z} = 0.7 ± 0.1, based on the estimate of constant potential surfaces (Loebman et al. 2014).
Fig. 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 (edgeon view of the disk), (b) prolate with the major axis on the plane of disk (faceon view of the disk), (c) oblate with the minor axis perpendicular to the disk (edgeon view of the disk), and (d) prolate with the major axis perpendicular to the disk (edgeon view of the disk). Only in case (b) are the axes of symmetry of the disk and the halo not aligned. 
Fig. 4. Components of the velocity of a star in the spherical polar coordinates with origin at the Galactic center (GC): v = v_{r} + v_{θ} + v_{ϕ}. 
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 γ = 1 (γ = 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 M_{⊙} stars, the minimum galactocentric distance is 15 kpc, which corresponds to the minimum ejection velocity ∼750 km s^{−1} for models in Newtonian gravity and QUMOND with γ = 1. For QUMOND with γ = 2, the ejection velocity threshold is ∼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 threebody 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 threebody numerical code that reproduces the encounter of a set of equalmass 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_{BH} = 4 × 10^{6} M_{⊙} for the black hole (as in our model in Sect. 3.1) and a mass of 4 M_{⊙} 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_{⊙} (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_{ej} ∼ 510 km s^{−1}, extends to velocities of ∼4000 km s^{−1}, and has a positive skewness; 16% (12%) of the ejected stars have speed v_{ej} ≳ 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 massindependent.
The direction of the ejection velocity is assigned randomly to each star, because Hills’ mechanism yields isotropic ejections. We used the fourth order RungeKutta method with adaptive stepsize to integrate the equation of motion of each star. We started with a predefined timestep of 5 kyr and the timestep 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 timestep 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_{⊙} on the main sequence is τ_{L} ≃ 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 τ_{ej} = ϵτ_{L}, where ϵ is a random number drawn from a uniform density distribution between 0 and 1. We took the average ejection rate to be 10^{−4} 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 Δ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 ith star was ejected at t = (i − 1)Δt with its age at the time of ejection. We chose the time of observation, which is the total runtime of the simulation, to be t_{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 . If a star satisfies this condition, its travel time is determined by [t_{obs} − (i − 1)Δ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_{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).
Simulations for QUMOND.
Simulations for Newtonian gravity.
We compared the final phasespace distributions of the HVSs obtained in QUMOND with the distributions obtained in Newtonian gravity. Due to the large radial velocities of the HVSs, their phasespace 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_{θ} = r(dθ/dt), and the azimuthal component v_{ϕ} = r sin θ(dϕ/dt), 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_{ej} ≳ 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.
Fig. 5. Distributions of HVSs in three sections of phase space in QUMOND with γ = 1.0 (blue dots) and Newtonian gravity with a prolate (q_{y} = 1, q_{z} = 1.1) dark matter halo (red dots). Left, middle, and right panels: r − v_{r}, z − v_{z}, and r − v_{θ} sections, respectively. Here, r is the galactocentric distance, z is the vertical coordinate, and v_{r}, v_{z}, v_{θ} are the radial, vertical, and latitudinal components of the velocity, respectively. 
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_{⊙}, ejection velocities v_{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_{⊙} 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_{θ} and v_{ϕ}.
We considered both Newtonian gravity with different shapes of the dark matter halo and QUMOND with γ = 1 and γ = 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_{θ} values for all the models. However, in QUMOND with a spherical bulge, the stars do not obtain any v_{ϕ} because the baryonic matter, the only component present, is axisymmetric. In our models of Newtonian gravity, the stars acquire nonzero values of v_{ϕ} 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_{ϕ} vanishes, and the only nonnull component of the tangential velocity is the latitudinal velocity, v_{θ}.
The left panel of Fig. 6 shows that v_{θ} can be as large as 600 km s^{−1} for stars with v_{ej} < 710 km s^{−1} in QUMOND with γ = 2, and for stars with v_{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_{θ} are consistently lower than ∼100 km s^{−1}. The right panel of Fig. 6 shows qualitatively similar results for v_{ϕ} in the Newtonian models with the misaligned axes of symmetry (q_{y} ≠ 1): in these models, v_{ϕ} is lower than ∼20 km s^{−1} when v_{ej} ≳ 750 km s^{−1}. In these Newtonian models with q_{y} ≠ 1, the vertical scatter in the maximum values of v_{θ} (left panel) is large. The stars ejected with smaller angle with respect to the disk (i.e., (90° −θ)≲30°) undergo larger deflection and attain larger v_{θ}. On the other hand, the stars ejected almost perpendicular to the disk do not bend significantly, hence have smaller v_{θ}. Conversely, in the same Newtonian models, the stars always attain a v_{ϕ}≳350 km s^{−1} (right panel) irrespective of the direction of the HVS ejection velocity.
Fig. 6. Left panel: maximum values of the latitudinal velocity components, v_{θ}, of 4 M_{⊙} HVSs as a function of their ejection velocities, v_{ej}, in Newtonian gravity with different shapes of the dark matter halo and in QUMOND with γ = 1 and γ = 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_{ej} = 710 and v_{ej} = 750 km s^{−1}. In all the models except QUMOND with γ = 2, the maximum v_{θ} substantially drops for stars with v_{ej} ≳ 750 km s^{−1}. For QUMOND with γ = 2, the maximum v_{θ} drops for stars with v_{ej} ≳ 710 km s^{−1}. Right panel: maximum values of the azimuthal velocity components, v_{ϕ}, of 4 M_{⊙} HVSs as a function of their ejection velocities, v_{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_{ϕ} values. The vertical scatter originates from the random directions of the star ejection. The vertical dotted line indicates v_{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_{ϕ} as a function of the radial component v_{r} for three stars with ejection velocity v_{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_{ej} and v_{ϕ} = 0, initially; at later times, v_{ϕ} 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_{ϕ} 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_{ϕ} reaches its maximum, as shown for v_{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_{ϕ} falls sharply around v_{ej} ≈ 750 km s^{−1} (right panel of Fig. 6).
Fig. 7. Azimuthal velocity component, v_{ϕ}, as a function of the radial velocity component, v_{r}, for three HVSs with mass 4 M_{⊙} and ejection velocity v_{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_{ej} and v_{ϕ} = 0. The maxima of v_{ϕ} are quite different in the three cases. The stars with v_{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_{ej} = 745 km s^{−1}. 
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_{θ} and v_{ϕ}, 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_{⊙} stars, this threshold ejection velocity is ∼710 km s^{−1} for QUMOND with γ = 2 and ∼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 γ = 1 can stars with v_{ej} ≲ 750 km s^{−1} travel beyond a distance of about 15 kpc from the Galactic center. Similarly, in QUMOND with γ = 2, the stars with v_{ej} ≲ 710 km s^{−1} cannot travel beyond 15 kpc. Therefore, identifying the 4 M_{⊙} 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 M_{⊙} HVSs at galactocentric distance r > 15 kpc.
Fig. 8. HVS outer turnaround as a function of the ejection speed, v_{ej}, in Newtonian gravity with different shapes of the dark matter halo and in QUMOND with γ = 1 and γ = 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_{ej}. The vertical dotted lines mark the thresholds v_{ej} = 710 km s^{−1} for QUMOND with γ = 2 and v_{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_{⊙} HVSs. Less massive stars have larger velocity thresholds and larger minimum galactocentric distances. 
The threshold ejection velocity increases with decreasing HVS mass, because of the longer lifetimes of lowermass HVSs: longerlived 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 M_{⊙} or 2.5 M_{⊙} stars with lifetimes on the main sequence ∼350 Myr and ∼580 Myr, respectively (Schaller et al. 1992), the threshold ejection velocities are ∼790 km s^{−1} and ∼815 km s^{−1} in Newtonian models; these threshold ejection speeds correspond to outer turnaround radii of ∼30 kpc and ∼50 kpc, respectively (Fig. 8). In QUMOND with γ = 1 (γ = 2), for 3 M_{⊙} or 2.5 M_{⊙} stars, the threshold ejection velocity is ∼800 (755) km s^{−1} and ∼830 (785) km s^{−1}, respectively; the corresponding turnaround radius is ∼30 (30) kpc and ∼50 (50) kpc, respectively. Therefore, the 3 M_{⊙} or 2.5 M_{⊙} stars are of interest for our test only if they are at galactocentric distance r larger than 30 kpc and 50 kpc, respectively.
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_{ϕ} and v_{θ} respectively, in both QUMOND and Newtonian gravity. We show that, for HVSs within ∼60 kpc, QUMOND yields upper limits for v_{ϕ} values that are substantially lower than the values that nonaxisymmetric gravitational potentials can generate in Newtonian gravity. The QUMOND scenario might thus be challenged if values of v_{ϕ} 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 nonnull v_{θ} values, whereas v_{ϕ} is always zero. Nonnull values of v_{ϕ} only appear in a nonaxisymmetric 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) and (6)), nonnull values of v_{ϕ} 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_{ϕ} 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 ∼5 kpc of the center.
This fundamental difference suggests that we might expect substantially higher values of v_{ϕ} 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_{ϕ} and v_{θ} in QUMOND and Newtonian gravity, and how v_{ϕ} can discriminate between the two theories of gravity.
6.1. The upper limit of the azimuthal component, v_{ϕ}, in QUMOND
In QUMOND, if the baryonic components are nonaxisymmetric, the HVSs may have nonnull v_{ϕ} values. To estimate these values of v_{ϕ}, 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_{ϕ} (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 ≳ 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 M_{⊙} stars with ejection velocity v_{ej} ≲ 710 km s^{−1} in QUMOND with γ = 2, or with v_{ej} ≲ 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_{ϕ} is inversely proportional to the ejection velocity, as expected by the argument illustrated in Sect. 5.2.
In order to find the maximum v_{ϕ} in QUMOND, we chose the lowest ejection velocity relevant for our analysis. Figure 9 shows v_{ϕ} as a function of the galactocentric distance r for three different ejection velocities: v_{ej} = 710 km s^{−1} in QUMOND with γ = 2 (red curves), v_{ej} = 750 km s^{−1} and 1000 km s^{−1} in QUMOND with γ = 1 (red and orange curves). All the stars shown in Fig. 9 are ejected at a polar angle θ = 45°; the solid lines show v_{ϕ} for an azimuthal angle ϕ = 45°, whereas the dashed and dotdashed lines show v_{ϕ} for the slowest stars for different ϕ. 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_{ϕ} we should expect is ∼8 km s^{−1}. We note that this result is not sensitive to the choice of the interpolating function ν (Eq. (2)), because within 5 kpc the gravitational field is mostly Newtonian.
Fig. 9. Magnitude of the azimuthal velocity component, v_{ϕ}, in QUMOND as a function of the galactocentric distance, r, for three 4 M_{⊙} HVSs with different ejection velocities: v_{ej} = 710 km s^{−1} in QUMOND with γ = 2 (red), v_{ej} = 750 km s^{−1} in QUMOND with γ = 1 (blue), and v_{ej} = 1000 km s^{−1} in QUMOND with γ = 1 (orange). The solid lines are for the stars ejected along θ = 45° and ϕ = 45°, whereas the dashed (dashdotted) line is for θ = 45° and ϕ = 15° (75°). In these examples, the QUMOND field is due to the triaxial bulge and the axisymmetric disk. The stars acquire v_{ϕ} only due to the triaxiality of the bulge. 
Assuming that the angular momentum per unit mass ℓ_{z}=r sin θv_{ϕ} is conserved at larger radii, and v_{ϕ} thus falls off as r^{−1}, we obtain
for 4 M_{⊙} stars. The upperlimit radial profile of v_{ϕ} 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_{⊙} 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 (θ = 90°; see Fig. 4), so that θ approaches 90°^{2}. As a result, sin θ always increases and causes a further decrease in v_{ϕ}, because ℓ_{z}=r sin θv_{ϕ} remains constant. Therefore, using the HVS with the lowest ejection velocity and taking sin θ as a constant yield the conservative upper limit of v_{ϕ} given in Eq. (13).
Fig. 10. Magnitudes of the azimuthal velocity components v_{ϕ} of 4 M_{⊙} 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_{ϕ} in QUMOND due to the triaxial bulge. The upper limit also holds for HVSs with masses lower than 4 M_{⊙}, provided they are beyond a certain galactocentric distance, e.g., r ≳ 30 kpc (50 kpc) for 3 M_{⊙} (2.5 M_{⊙}) HVSs. 
The radial profile of Eq. (13) also provides a conservative upper limit for stars with masses lower than 4 M_{⊙}. 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_{ϕ}(r) in Eq. (13), will be lower.
In Newtonian gravity with an axisymmetric distribution of baryons (Eqs. (4)–(6)), the HVSs obtain nonnull v_{ϕ} values only when the dark matter halo is nonaxisymmetric about the z axis, as in our models with q_{y} ≠ 1. Figure 10 shows the values of v_{ϕ} 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 ≳ r_{a} = 21.6 kpc. Thus, the azimuthal angular momentum, ℓ_{z} = r sin θv_{ϕ}, of each HVS becomes a constant as the star travels beyond r ∼ r_{a} = 21.6 kpc and v_{ϕ} decreases as r^{−1}. Near r ∼ 21.6 kpc, v_{ϕ} 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 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 ∼10^{6} K is suggested by the O VII and O VIII emission and absorption lines in the soft Xray 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^{3} (Fukugita et al. 1998; Gupta et al. 2012; Shull et al. 2012; Fang et al. 2013; Planck Collaboration VI 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)
where r_{0, HG} = 100 kpc is the core radius and r_{t, HG} = 200 kpc is the truncation radius. The elliptical radius m for the oblate halo is
Currently, the shape of the halo has no observational constraints; we thus varied a_{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 ρ_{0, HG}: 3.0 × 10^{5} and 3.0 × 10^{4} M_{⊙} kpc^{−3}. These values yield a total mass of 1.5 × 10^{11} M_{⊙} and 1.5 × 10^{10} M_{⊙}, 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 and 12). In both cases, we considered the axisymmetric models for the central black hole, bulge and disk of (Eqs. (4)–(6)).
Fig. 11. Magnitudes of the azimuthal velocity components, v_{ϕ}, of 4 M_{⊙} 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 × 10^{11} M_{⊙} within 100 kpc in both Newtonian and QUMOND scenarios. Left and right panels: results for two different oblatenesses, a_{HG} = 0.4 and 0.8, of the HG halo (see Eqs. (14) and (15)). 
Fig. 12. Same as Fig. 11 except for the HG halo with a lower mass: 1.5 × 10^{10} M_{⊙} within 100 kpc. 
In Newtonian gravity, we considered the nonaxisymmetric 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 ≳ 21.6 kpc) and its contribution to the v_{ϕ} values thus falls as r^{−1} at large distances, the constant oblateness of the HG halo increases the v_{ϕ} values at all radii. Nevertheless the v_{ϕ} values are dominated by the shape of the nonaxisymmetric 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_{ϕ}. Therefore, for the HG halo with the highest mass M(< 100 kpc) = 1.5 × 10^{11} M_{⊙} (Fig. 11), although the v_{ϕ}’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_{ϕ} 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 kpc) = 1.5 × 10^{10} M_{⊙} (Fig. 12), the v_{ϕ} 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_{ϕ}’s in QUMOND at all distances, compared to the presence of a triaxial bulge alone. However, at smaller distances (15 kpc ≲ r ≲ 60 kpc), the v_{ϕ} values in Newtonian gravity due to a nonaxisymmetric dark matter halo may still substantially exceed the QUMOND values.
6.2. Latitudinal component, v_{θ}
The difference in the latitudinal component v_{θ} of the velocity in QUMOND and in Newtonian gravity is subtler than the difference in v_{ϕ}. 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_{θ} lower than in QUMOND; on the contrary, an oblate halo with its major axis in the plane of the disk will generate v_{θ} 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_{θ} in Newtonian gravity, assuming that the dark matter halo is spherical and therefore it does not affect v_{θ}. Similarly, in our adopted model, the SMBH and the bulge have spherically symmetric potential and do not contribute to v_{θ}. The shaded histogram in the left panel of Fig. 13 shows the distribution of v_{θ} 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_{θ} in QUMOND than in Newtonian gravity.
Fig. 13. Distributions of the latitudinal component of the velocity, v_{θ}, in Newtonian gravity and in QUMOND for 4 M_{⊙} HVSs. Left panel: distributions of v_{θ} for different models of the Galactic potential: Newtonian gravity with a spherical dark matter halo (gray shaded histogram) and QUMOND with γ = 1 (orange histogram) and γ = 2 (red histogram). As the gravitational pull of the baryonic disk is enhanced in QUMOND, the fraction of HVSs with high v_{θ} is larger than in Newtonian gravity with a spherical halo. Right panel: distributions of v_{θ} 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_{θ} and the smaller fraction of HVSs with low v_{θ} (compare the shaded and red histograms). The opposite occurs with a prolate halo (compare the shaded and blue histograms). 
Figure 14 shows the magnitude of v_{θ} as a function of the galactocentric distance r due to the baryonic components in QUMOND with γ = 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_{θ} decrease as ∼r^{−1}, because, as each star travels beyond the length scale of the stellar disk, its angular momentum, ∼rv_{θ}, becomes a constant. Because of the enhanced pull of QUMOND, the maximum value of v_{θ} is higher than in Newtonian gravity at all radii. Unfortunately, unlike the case of v_{ϕ}, 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.
Fig. 14. Latitudinal velocity components, v_{θ}, as a function of the galactocentric distance, r, of 4 M_{⊙} HVSs in QUMOND with γ = 1 (orange dots) and γ = 2 (red dots), and in Newtonian gravity with a spherical dark matter halo (q_{z} = 1) (black dots). 
In Newtonian gravity with a dark matter halo, nonnull v_{θ} values are caused by the Galactic disk and by the halo with q_{z} ≠ 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_{θ} 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_{ϕ}, QUMOND sets an upper limit to v_{θ}; on the contrary, in Newtonian gravity, v_{θ} may be higher than this upper limit if q_{z} is sufficiently small. In our models, the Newtonian v_{θ} values are higher than the QUMOND upper limit if q_{z} ≲ 0.6. However, the difference between the QUMOND and Newtonian v_{θ} values is not as prominent as for v_{ϕ} values: In our model, the transition scale of the halo shape from oblate to spherical is r_{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_{ϕ} values in QUMOND is much smaller than the scale length of the nonaxisymmetric halo responsible for nonnull v_{ϕ}’s in Newtonian gravity. Therefore, v_{θ} is less effective than v_{ϕ} at distinguishing QUMOND from Newtonian gravity. Figure 15 quantifies this difficulty by showing v_{θ} versus r for QUMOND with γ = 1 and γ = 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).
Fig. 15. Latitudinal velocity components, v_{θ}, as a function of the galactocentric distance, r, of 4 M_{⊙} HVSs in QUMOND with γ = 1 (orange dots) and γ = 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 (Eq. (13)) with the measured of nine stars drawn from the HVS survey sample of Brown et al. (2014). These nine stars have masses in the range ∼2.5 − 4 M_{⊙}. Their galactocentric distances are larger than the minimum massdependent radii required for our test (see Sect. 5.2). In addition, Kenyon et al. (2018) show that HVSs ejected within 25° 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 ∼84° 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 selfconsistent 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, , of their galactocentric velocities from the proper motions available in the Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2016, 2021). 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 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 ∼50 − 340%. For the ninth star, B329, neglecting the error on its distance is inappropriate. B329 is a (3.21 ± 0.24) M_{⊙} star at galactocentric distance r = (61 ± 13) kpc; by ignoring the uncertainty on the distance, we derived km s^{−1}, whereas including the uncertainty on the distance yields km s^{−1}.
We conclude that the current uncertainties on make the azimuthal components consistent with zero within 3σ 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 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, , 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_{ϕ}, cannot exceed a velocity threshold in QUMOND, while v_{ϕ} 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 nonnull azimuthal tangential speeds, v_{ϕ}, due to the nonaxisymmetric 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_{ϕ} 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_{ϕ}, 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 M_{⊙} stars is ∼710 km s^{−1} in QUMOND with γ = 2, and it is ∼750 km s^{−1} both in QUMOND with γ = 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 ∼15 kpc if they live long enough. We thus expect that the v_{ϕ} component of 4 M_{⊙} 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 M_{⊙} and 2.5 M_{⊙} HVSs have ejection velocity thresholds of ∼790 km s^{−1} and ∼815 km s^{−1}, which correspond to the minimum galactocentric distances of 30 kpc and 50 kpc, respectively. In QUMOND with γ = 1 (γ = 2), 3 M_{⊙} and 2.5 M_{⊙} HVSs have ejection velocity thresholds of ∼800 (755) km s^{−1} and ∼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 ∼21.6 kpc. It follows that v_{ϕ} is affected by the nonaxisymmetric 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_{ϕ}. Additional variations in v_{ϕ} might derive from massive external objects, such as the LMC (Kenyon et al. 2018), whereas other nonaxisymmetric components within the disk, such as spiral arms or density inhomogeneities, are expected to affect only the stars within ∼1 kpc of the disk (Gardner et al. 2020). The contribution from the HG halo in v_{ϕ} 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 ∼ 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_{ϕ} is proportional to r^{−1} (Eq. (13)). Therefore, when compared with the Newtonian model affecting v_{ϕ} out to ∼21.6 kpc, the values of v_{ϕ} of HVSs at r ≳ 15 kpc in QUMOND are substantially lower than the values allowed for a nonaxisymmetric potential in Newtonian gravity. For example, at r ∼ 20 kpc, we find v_{ϕ}≲2 km s^{−1} in QUMOND, whereas v_{ϕ} 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_{ϕ} values of the HVSs in QUMOND, and the v_{ϕ} values in this case may be higher than those generated by the triaxial bulge alone. On the contrary, in Newtonian gravity, the v_{ϕ} 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_{ϕ} for a few 4 M_{⊙} HVSs at galactocentric distances larger than ∼15 kpc (or 3 M_{⊙} HVSs at r ≳ 30 kpc, or 2.5 M_{⊙} HVSs at r ≳ 50 kpc) and smaller than ∼60 kpc may in principle test the validity of QUMOND. Finding a few HVSs with azimuthal components, v_{ϕ}, above the QUMOND upper limits, , given by Eq. (13) (black line in Fig. 10) or shown in Figs. 11 and 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_{ϕ}, 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_{ϕ} thus needs to be improved by at least a factor of ∼10 to make our test decisive. Future measurements from spaceborne astrometric missions with expected microarcsecond precision on star positions, such as Theia (The Theia Collaboration 2017; Malbet et al. 2019, 2021), are expected to allow us to discriminate between the two theories of gravity.
The baryonic census in the Local Group appears to add up to only ∼15% (Fukugita et al. 1998) of the baryonic mass expected from the estimated baryonic abundance (Ω_{b}h^{2} ≈ 0.022) (Planck Collaboration VI 2020) from the Big Bang nucleosynthesis and the CMB anisotropies.
Acknowledgments
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 (UniTOCSP), 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 2013, 2018) GALPY (Bovy 2015).
References
 Aaronson, M. 1983, ApJ, 266, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Abadi, M. G., Navarro, J. F., & Steinmetz, M. 2009, ApJ, 691, L63 [Google Scholar]
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 Begeman, K. G., Broeils, A. H., & Sanders, R. H. 1991, MNRAS, 249, 523 [Google Scholar]
 Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509 [Google Scholar]
 Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Bertone, G., Hooper, D., & Silk, J. 2005, Phys. Rep., 405, 279 [Google Scholar]
 Bhattacharjee, P., Chaudhury, S., & Kundu, S. 2014, ApJ, 785, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Bienaymé, O., Famaey, B., Wu, X., Zhao, H. S., & Aubert, D. 2009, A&A, 500, 801 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Binney, J., Gerhard, O., & Spergel, D. 1997, MNRAS, 288, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Biviano, A. 2000, in Constructing the Universe with Clusters of Galaxies, ed. F. Durret, & D. Gerbal, 1 [Google Scholar]
 BlandHawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [Google Scholar]
 Blumenthal, G. R., Pagels, H., & Primack, J. R. 1982, Nature, 299, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Nature, 311, 517 [Google Scholar]
 Boehle, A., Ghez, A. M., Schödel, R., et al. 2016, ApJ, 830, 17 [Google Scholar]
 Böhringer, H., & Werner, N. 2010, A&ARv, 18, 127 [Google Scholar]
 Bond, J. R., Silk, J., & Kolb, E. W. 1982, ApJ, 255, 341 [NASA ADS] [CrossRef] [Google Scholar]
 Bosma, A. 1978, PhD Thesis, University of Groningen, the Netherlands [Google Scholar]
 Bottema, R., Pestaña, J. L. G., Rothberg, B., & Sanders, R. H. 2002, A&A, 393, 453 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Bromley, B. C., Kenyon, S. J., Geller, M. J., et al. 2006, ApJ, 653, 1194 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, W. R. 2015, ARA&A, 53, 15 [Google Scholar]
 Brown, W. R., Geller, M. J., Kenyon, S. J., & Kurtz, M. J. 2005, ApJ, 622, L33 [Google Scholar]
 Brown, W. R., Geller, M. J., Kenyon, S. J., & Kurtz, M. J. 2006, ApJ, 647, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, W. R., Geller, M. J., & Kenyon, S. J. 2009, ApJ, 690, 1639 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, W. R., Geller, M. J., & Kenyon, S. J. 2012, ApJ, 751, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, W. R., Geller, M. J., & Kenyon, S. J. 2014, ApJ, 787, 89 [Google Scholar]
 Brown, W. R., Lattanzi, M. G., Kenyon, S. J., & Geller, M. J. 2018, ApJ, 866, 39 [Google Scholar]
 Cappellari, M. 2016, ARA&A, 54, 597 [Google Scholar]
 Cappellari, M., Bacon, R., Bureau, M., et al. 2006, MNRAS, 366, 1126 [Google Scholar]
 Chae, K.H., Lelli, F., Desmond, H., et al. 2020, ApJ, 904, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Clerc, N., Kirkpatrick, C. C., Finoguenov, A., et al. 2020, MNRAS, 497, 3976 [NASA ADS] [CrossRef] [Google Scholar]
 Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep., 513, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Contigiani, O., Rossi, E. M., & Marchetti, T. 2019, MNRAS, 487, 4025 [NASA ADS] [CrossRef] [Google Scholar]
 Cuddeford, P., & Amendt, P. 1991, MNRAS, 253, 427 [NASA ADS] [Google Scholar]
 Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [Google Scholar]
 Deason, A. J., Belokurov, V., Evans, N. W., et al. 2012, MNRAS, 425, 2840 [Google Scholar]
 Debattista, V. P., Roškar, R., Valluri, M., et al. 2013, MNRAS, 434, 2971 [NASA ADS] [CrossRef] [Google Scholar]
 de Martino, I., Chakrabarty, S. S., Cesare, V., et al. 2020, Universe, 6, 107 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, T., & Franx, M. 1991, ARA&A, 29, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Diaferio, A., Schindler, S., & Dolag, K. 2008, Space Sci. Rev., 134, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Du, C., Li, H., Yan, Y., et al. 2019, ApJS, 244, 4 [Google Scholar]
 Edelmann, H., Napiwotzki, R., Heber, U., Christlieb, N., & Reimers, D. 2005, ApJ, 634, L181 [Google Scholar]
 Eisenhauer, F., Genzel, R., Alexander, T., et al. 2005, ApJ, 628, 246 [Google Scholar]
 Erkal, D., Boubert, D., Gualandris, A., Evans, N. W., & Antonini, F. 2019, MNRAS, 483, 2007 [NASA ADS] [CrossRef] [Google Scholar]
 Fall, S. M., & Efstathiou, G. 1980, MNRAS, 193, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Relativ., 15, 10 [Google Scholar]
 Famaey, B., McGaugh, S., & Milgrom, M. 2018, MNRAS, 480, 473 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, T., Bullock, J., & BoylanKolchin, M. 2013, ApJ, 762, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Fragione, G., & Loeb, A. 2017, New Astron., 55, 32 [CrossRef] [Google Scholar]
 Franx, M., Illingworth, G., & de Zeeuw, T. 1991, ApJ, 383, 112 [CrossRef] [Google Scholar]
 Frenk, C. S., & White, S. D. M. 2012, Ann. Phys., 524, 507 [NASA ADS] [CrossRef] [Google Scholar]
 Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Luri, X., et al.) 2021, A&A, 649, A7 [EDP Sciences] [Google Scholar]
 Gallo, A., Ostorero, L., Chakrabarty, S. S., et al. 2021, A&A, submitted [arXiv:2111.09657] [Google Scholar]
 Gardner, S., Hinkel, A., & Yanny, B. 2020, ApJ, 890, 110 [CrossRef] [Google Scholar]
 Gatto, A., Fraternali, F., Read, J. I., et al. 2013, MNRAS, 433, 2749 [NASA ADS] [CrossRef] [Google Scholar]
 Geller, M. J., Diaferio, A., Rines, K. J., & Serra, A. L. 2013, ApJ, 764, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Genzel, R., Eisenhauer, F., & Gillessen, S. 2010, Rev. Mod. Phys., 82, 3121 [Google Scholar]
 Ghez, A. M., Salim, S., Weinberg, N. N., et al. 2008, ApJ, 689, 1044 [Google Scholar]
 Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30 [Google Scholar]
 Giodini, S., Pierini, D., Finoguenov, A., et al. 2009, ApJ, 703, 982 [NASA ADS] [CrossRef] [Google Scholar]
 Girardi, M., Manzato, P., Mezzetti, M., Giuricin, G., & Limboz, F. 2002, ApJ, 569, 720 [NASA ADS] [CrossRef] [Google Scholar]
 Gnedin, O. Y., Gould, A., MiraldaEscudé, J., & Zentner, A. R. 2005, ApJ, 634, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Gupta, A., Mathur, S., Krongold, Y., Nicastro, F., & Galeazzi, M. 2012, ApJ, 756, L8 [CrossRef] [Google Scholar]
 Haghi, H., Bazkiaei, A. E., Zonoozi, A. H., & Kroupa, P. 2016, MNRAS, 458, 4172 [NASA ADS] [CrossRef] [Google Scholar]
 Hammer, F., Yang, Y., Fouquet, S., et al. 2013, MNRAS, 431, 3543 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Hattori, K., Valluri, M., Bell, E. F., & Roederer, I. U. 2018a, ApJ, 866, 121 [CrossRef] [Google Scholar]
 Hattori, K., Valluri, M., & Castro, N. 2018b, ApJ, 869, 33 [CrossRef] [Google Scholar]
 Hayashi, E., Navarro, J. F., & Springel, V. 2007, MNRAS, 377, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Hernandez, X., Sussman, R. A., & Nasser, L. 2019, MNRAS, 483, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
 Hesp, C., & Helmi, A. 2018, ArXiv eprints [arXiv:1804.03670] [Google Scholar]
 Hills, J. G. 1988, Nature, 331, 687 [Google Scholar]
 Hirsch, H. A., Heber, U., O’Toole, S. J., & Bresolin, F. 2005, A&A, 444, L61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hodson, A. O., Diaferio, A., & Ostorero, L. 2020, A&A, 640, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hohl, F. 1971, ApJ, 168, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, Y., Liu, X. W., Zhang, H. W., et al. 2017, ApJ, 847, L9 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Irwin, M., & Hatzidimitriou, D. 1995, MNRAS, 277, 1354 [NASA ADS] [CrossRef] [Google Scholar]
 Kafle, P. R., Sharma, S., Lewis, G. F., & BlandHawthorn, J. 2014, ApJ, 794, 59 [Google Scholar]
 Kenyon, S. J., Bromley, B. C., Geller, M. J., & Brown, W. R. 2008, ApJ, 680, 312 [NASA ADS] [CrossRef] [Google Scholar]
 Kenyon, S. J., Bromley, B. C., Brown, W. R., & Geller, M. J. 2014, ApJ, 793, 122 [Google Scholar]
 Kenyon, S. J., Bromley, B. C., Brown, W. R., & Geller, M. J. 2018, ApJ, 864, 130 [CrossRef] [Google Scholar]
 Kneib, J. P., Ellis, R. S., Smail, I., Couch, W. J., & Sharples, R. M. 1996, ApJ, 471, 643 [Google Scholar]
 Koposov, S. E., Boubert, D., Li, T. S., et al. 2020, MNRAS, 491, 2465 [Google Scholar]
 Kormendy, J. 1987, in Dark Matter in the Universe, eds. J. Kormendy, & G. R. Knapp, 117, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Launhardt, R., Zylka, R., & Mezger, P. G. 2002, A&A, 384, 112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, ApJ, 816, L14 [Google Scholar]
 Li, Y., Luo, A., Zhao, G., et al. 2012, ApJ, 744, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Y.B., Luo, A. L., Zhao, G., et al. 2015, Res. Astron. Astrophys., 15, 1364 [Google Scholar]
 Li, Y.B., Luo, A. L., Lu, Y.J., et al. 2021, ApJS, 252, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Loebman, S. R., Ivezić, Ž., Quinn, T. R., et al. 2014, ApJ, 794, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Łokas, E. L. 2009, MNRAS, 394, L102 [NASA ADS] [CrossRef] [Google Scholar]
 Luna, A., Minniti, D., & AlonsoGarcía, J. 2019, ApJ, 887, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Malbet, F., Abbas, U., Alves, J., et al. 2019, ArXiv eprints [arXiv:1910.08028] [Google Scholar]
 Malbet, F., Boehm, C., KroneMartins, A., et al. 2021, Exp. Astron., 51, 845 [NASA ADS] [CrossRef] [Google Scholar]
 Marchetti, T., Rossi, E. M., Kordopatis, G., et al. 2017, MNRAS, 470, 1388 [NASA ADS] [CrossRef] [Google Scholar]
 Martinsson, T. P. K., Verheijen, M. A. W., Westfall, K. B., et al. 2013, A&A, 557, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Massey, P., Levine, S. E., Neugent, K. F., et al. 2018, AJ, 156, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Mateo, M. L. 1998, ARA&A, 36, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Mateo, M., Olszewski, E., Welch, D. L., Fischer, P., & Kunkel, W. 1991, AJ, 102, 914 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. 2020, Galaxies, 8, 35 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2004, ApJ, 609, 652 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2008, ApJ, 683, 137 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2011, Phys. Rev. Lett., 106, 121303 [CrossRef] [Google Scholar]
 McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. 2000, ApJ, 533, L99 [NASA ADS] [CrossRef] [Google Scholar]
 Meneghetti, M., Davoli, G., Bergamini, P., et al. 2020, Science, 369, 1347 [Google Scholar]
 Merritt, D. 2020, A Philosophical Approach to MOND: Assessing the Milgromian Research Program in Cosmology (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Milgrom, M. 1983a, ApJ, 270, 365 [Google Scholar]
 Milgrom, M. 1983b, ApJ, 270, 371 [Google Scholar]
 Milgrom, M. 1983c, ApJ, 270, 384 [Google Scholar]
 Milgrom, M. 2009, Phys. Rev. D, 80, 123536 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2010, MNRAS, 403, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2015, Can. J. Phys., 93, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. J., & Bregman, J. N. 2013, ApJ, 770, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
 Neugent, K. F., Massey, P., Morrell, N. I., Skiff, B., & Georgy, C. 2018, AJ, 155, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Nipoti, C., Londrillo, P., Zhao, H., & Ciotti, L. 2007, MNRAS, 379, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Nojiri, S., Odintsov, S. D., & Oikonomou, V. K. 2017, Phys. Rep., 692, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467 [Google Scholar]
 Paerels, F. B. S., & Kahn, S. M. 2003, ARA&A, 41, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Pawlowski, M. S., Kroupa, P., & de Boer, K. S. 2011, A&A, 532, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peebles, P. J. E. 1982, ApJ, 263, L1 [Google Scholar]
 Pereira, C. B., Jilinski, E. G., Drake, N. A., Ortega, V. G., & Roig, F. 2013, A&A, 559, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perets, H. B., Wu, X., Zhao, H. S., et al. 2009, ApJ, 697, 2096 [NASA ADS] [CrossRef] [Google Scholar]
 Persic, M., Salucci, P., & Stel, F. 1996, MNRAS, 281, 27 [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25 [Google Scholar]
 Pota, V., Forbes, D. A., Romanowsky, A. J., et al. 2013, MNRAS, 428, 389 [NASA ADS] [CrossRef] [Google Scholar]
 PriceWhelan, A. M., Hogg, D. W., Johnston, K. V., & Hendel, D. 2014, ApJ, 794, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Proctor, R. N., Mendes de Oliveira, C., Azanha, L., Dupke, R., & Overzier, R. 2015, MNRAS, 449, 2345 [NASA ADS] [CrossRef] [Google Scholar]
 Pulsoni, C., Gerhard, O., Arnaboldi, M., et al. 2018, A&A, 618, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, Nature, 352, 411 [Google Scholar]
 Read, J. I., & Moore, B. 2005, MNRAS, 361, 971 [NASA ADS] [CrossRef] [Google Scholar]
 Rines, K., Geller, M. J., Diaferio, A., Kurtz, M. J., & Jarrett, T. H. 2004, AJ, 128, 1078 [NASA ADS] [CrossRef] [Google Scholar]
 Rines, K., Geller, M. J., Diaferio, A., & Kurtz, M. J. 2013, ApJ, 767, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Rosati, P., Borgani, S., & Norman, C. 2002, ARA&A, 40, 539 [Google Scholar]
 Rossi, E. M., Marchetti, T., Cacciato, M., Kuiack, M., & Sari, R. 2017, MNRAS, 467, 1844 [NASA ADS] [Google Scholar]
 Rubin, V. C., Ford, W. K., Jr., Thonnard, N., & Burstein, D. 1982, ApJ, 261, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Rubin, V. C., Burstein, D., Ford, W. K., Jr., & Thonnard, N. 1985, ApJ, 289, 81 [CrossRef] [Google Scholar]
 Salem, M., Besla, G., Bryan, G., et al. 2015, ApJ, 815, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H. 1990, A&ARv, 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H., & McGaugh, S. S. 2002, ARA&A, 40, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1 [Google Scholar]
 Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269 [Google Scholar]
 Sellwood, J. A. 2014, Rev. Mod. Phys., 86, 1 [Google Scholar]
 Sellwood, J. A., Shen, J., & Li, Z. 2019, MNRAS, 486, 4710 [NASA ADS] [CrossRef] [Google Scholar]
 Sereno, M., Giocoli, C., Izzo, L., et al. 2018, Nat. Astron., 2, 744 [NASA ADS] [CrossRef] [Google Scholar]
 Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Silk, J. 1967, Nature, 215, 1155 [NASA ADS] [CrossRef] [Google Scholar]
 Silk, J., AntonuccioDelogu, V., Dubois, Y., et al. 2012, A&A, 545, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, J. D. 2019, ARA&A, 57, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Skordis, C., & Zlosnik, T. 2021, Phys. Rev. Lett., 127, 161302 [NASA ADS] [CrossRef] [Google Scholar]
 Sofue, Y., & Rubin, V. 2001, ARA&A, 39, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Sohn, J., Geller, M. J., Zahid, H. J., et al. 2017, ApJS, 229, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Sohn, J., Geller, M. J., & Zahid, H. J. 2019, ApJ, 880, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D. N., Bean, R., Doré, O., et al. 2007, ApJS, 170, 377 [CrossRef] [Google Scholar]
 Strigari, L. E. 2013, Phys. Rep., 531, 1 [Google Scholar]
 Tanabashi, M., Hagiwara, K., Hikasa, K., et al. 2018, Phys. Rev. D, 98, 030001 [CrossRef] [Google Scholar]
 Taylor, M. B. 2005, in Astronomical Data Analysis Software and Systems XIV, eds. P. Shopbell, M. Britton, & R. Ebert, ASP Conf. Ser., 347, 29 [Google Scholar]
 The Theia Collaboration (Boehm, C., et al.) 2017, ArXiv eprints [arXiv:1707.01348] [Google Scholar]
 Thomas, G. F., Famaey, B., Ibata, R., Lüghausen, F., & Kroupa, P. 2017, A&A, 603, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tian, Y., Yu, P.C., Li, P., McGaugh, S. S., & Ko, C.M. 2021, ApJ, 910, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Tillich, A., Heber, U., Geier, S., et al. 2011, A&A, 527, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tyson, J. A., Valdes, F., & Wenk, R. A. 1990, ApJ, 349, L1 [CrossRef] [Google Scholar]
 Umetsu, K. 2020, A&ARv, 28, 7 [Google Scholar]
 van den Bosch, F. C., & Swaters, R. A. 2001, MNRAS, 325, 1017 [NASA ADS] [CrossRef] [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
 Vogelsberger, M., White, S. D. M., Helmi, A., & Springel, V. 2008, MNRAS, 385, 236 [NASA ADS] [CrossRef] [Google Scholar]
 Voit, G. M. 2005, Rev. Mod. Phys., 77, 207 [Google Scholar]
 Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2009, ApJ, 704, 1274 [Google Scholar]
 Walker, S., Simionescu, A., Nagai, D., et al. 2019, Space Sci. Rev., 215, 7 [Google Scholar]
 Wang, X., & Loeb, A. 2018, New Astron., 61, 95 [CrossRef] [Google Scholar]
 Wang, Y.H., Leigh, N., Yuan, Y.F., & Perna, R. 2018, MNRAS, 475, 4595 [NASA ADS] [CrossRef] [Google Scholar]
 Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435 [NASA ADS] [CrossRef] [Google Scholar]
 White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52 [Google Scholar]
 Xue, X. X., Rix, H. W., Zhao, G., et al. 2008, ApJ, 684, 1143 [Google Scholar]
 Yu, Q., & Tremaine, S. 2003, ApJ, 599, 1129 [Google Scholar]
 Zavala, J., & Frenk, C. S. 2019, Galaxies, 7, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, F., Lu, Y., & Yu, Q. 2013, ApJ, 768, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, H., Famaey, B., Lüghausen, F., & Kroupa, P. 2013, A&A, 557, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zheng, Z., Carlin, J. L., Beers, T. C., et al. 2014, ApJ, 785, L23 [Google Scholar]
All Tables
All Figures
Fig. 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. 

In the text 
Fig. 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, ϕ, which is the angle with respect to the x axis in the plane of the disk. We show the results for various values of ϕ; ϕ = 0° and ϕ = 90° correspond to the x and y axes, respectively. 

In the text 
Fig. 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 (edgeon view of the disk), (b) prolate with the major axis on the plane of disk (faceon view of the disk), (c) oblate with the minor axis perpendicular to the disk (edgeon view of the disk), and (d) prolate with the major axis perpendicular to the disk (edgeon view of the disk). Only in case (b) are the axes of symmetry of the disk and the halo not aligned. 

In the text 
Fig. 4. Components of the velocity of a star in the spherical polar coordinates with origin at the Galactic center (GC): v = v_{r} + v_{θ} + v_{ϕ}. 

In the text 
Fig. 5. Distributions of HVSs in three sections of phase space in QUMOND with γ = 1.0 (blue dots) and Newtonian gravity with a prolate (q_{y} = 1, q_{z} = 1.1) dark matter halo (red dots). Left, middle, and right panels: r − v_{r}, z − v_{z}, and r − v_{θ} sections, respectively. Here, r is the galactocentric distance, z is the vertical coordinate, and v_{r}, v_{z}, v_{θ} are the radial, vertical, and latitudinal components of the velocity, respectively. 

In the text 
Fig. 6. Left panel: maximum values of the latitudinal velocity components, v_{θ}, of 4 M_{⊙} HVSs as a function of their ejection velocities, v_{ej}, in Newtonian gravity with different shapes of the dark matter halo and in QUMOND with γ = 1 and γ = 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_{ej} = 710 and v_{ej} = 750 km s^{−1}. In all the models except QUMOND with γ = 2, the maximum v_{θ} substantially drops for stars with v_{ej} ≳ 750 km s^{−1}. For QUMOND with γ = 2, the maximum v_{θ} drops for stars with v_{ej} ≳ 710 km s^{−1}. Right panel: maximum values of the azimuthal velocity components, v_{ϕ}, of 4 M_{⊙} HVSs as a function of their ejection velocities, v_{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_{ϕ} values. The vertical scatter originates from the random directions of the star ejection. The vertical dotted line indicates v_{ej} = 750 km s^{−1}. 

In the text 
Fig. 7. Azimuthal velocity component, v_{ϕ}, as a function of the radial velocity component, v_{r}, for three HVSs with mass 4 M_{⊙} and ejection velocity v_{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_{ej} and v_{ϕ} = 0. The maxima of v_{ϕ} are quite different in the three cases. The stars with v_{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_{ej} = 745 km s^{−1}. 

In the text 
Fig. 8. HVS outer turnaround as a function of the ejection speed, v_{ej}, in Newtonian gravity with different shapes of the dark matter halo and in QUMOND with γ = 1 and γ = 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_{ej}. The vertical dotted lines mark the thresholds v_{ej} = 710 km s^{−1} for QUMOND with γ = 2 and v_{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_{⊙} HVSs. Less massive stars have larger velocity thresholds and larger minimum galactocentric distances. 

In the text 
Fig. 9. Magnitude of the azimuthal velocity component, v_{ϕ}, in QUMOND as a function of the galactocentric distance, r, for three 4 M_{⊙} HVSs with different ejection velocities: v_{ej} = 710 km s^{−1} in QUMOND with γ = 2 (red), v_{ej} = 750 km s^{−1} in QUMOND with γ = 1 (blue), and v_{ej} = 1000 km s^{−1} in QUMOND with γ = 1 (orange). The solid lines are for the stars ejected along θ = 45° and ϕ = 45°, whereas the dashed (dashdotted) line is for θ = 45° and ϕ = 15° (75°). In these examples, the QUMOND field is due to the triaxial bulge and the axisymmetric disk. The stars acquire v_{ϕ} only due to the triaxiality of the bulge. 

In the text 
Fig. 10. Magnitudes of the azimuthal velocity components v_{ϕ} of 4 M_{⊙} 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_{ϕ} in QUMOND due to the triaxial bulge. The upper limit also holds for HVSs with masses lower than 4 M_{⊙}, provided they are beyond a certain galactocentric distance, e.g., r ≳ 30 kpc (50 kpc) for 3 M_{⊙} (2.5 M_{⊙}) HVSs. 

In the text 
Fig. 11. Magnitudes of the azimuthal velocity components, v_{ϕ}, of 4 M_{⊙} 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 × 10^{11} M_{⊙} within 100 kpc in both Newtonian and QUMOND scenarios. Left and right panels: results for two different oblatenesses, a_{HG} = 0.4 and 0.8, of the HG halo (see Eqs. (14) and (15)). 

In the text 
Fig. 12. Same as Fig. 11 except for the HG halo with a lower mass: 1.5 × 10^{10} M_{⊙} within 100 kpc. 

In the text 
Fig. 13. Distributions of the latitudinal component of the velocity, v_{θ}, in Newtonian gravity and in QUMOND for 4 M_{⊙} HVSs. Left panel: distributions of v_{θ} for different models of the Galactic potential: Newtonian gravity with a spherical dark matter halo (gray shaded histogram) and QUMOND with γ = 1 (orange histogram) and γ = 2 (red histogram). As the gravitational pull of the baryonic disk is enhanced in QUMOND, the fraction of HVSs with high v_{θ} is larger than in Newtonian gravity with a spherical halo. Right panel: distributions of v_{θ} 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_{θ} and the smaller fraction of HVSs with low v_{θ} (compare the shaded and red histograms). The opposite occurs with a prolate halo (compare the shaded and blue histograms). 

In the text 
Fig. 14. Latitudinal velocity components, v_{θ}, as a function of the galactocentric distance, r, of 4 M_{⊙} HVSs in QUMOND with γ = 1 (orange dots) and γ = 2 (red dots), and in Newtonian gravity with a spherical dark matter halo (q_{z} = 1) (black dots). 

In the text 
Fig. 15. Latitudinal velocity components, v_{θ}, as a function of the galactocentric distance, r, of 4 M_{⊙} HVSs in QUMOND with γ = 1 (orange dots) and γ = 2 (red dots), and in Newtonian gravity with an oblate dark matter halo with q_{z} = 0.4 (dark cyan dots). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.