© A. Gallo et al. 2022

1. Introduction

The Lambda cold dark matter cosmological model predicts the existence of dark matter (DM) halos surrounding the Milky Way (MW) and the external galaxies. Although at first order the spherically symmetric Navarro–Frenk–White (NFW) profile (Navarro et al. 1997) can provide a good approximation of the shape of the DM halos, the first numerical N-body simulations (Frenk et al. 1988; Dubinski & Carlberg 1991; Warren et al. 1992; Cole & Lacey 1996) found the shape of the DM halos to be triaxial, and subsequent works (e.g., Jing & Suto 2002; Bailin & Steinmetz 2005; Hayashi et al. 2007; Vera-Ciro et al. 2011) confirmed this result. In these simulations, the DM halos are triaxial with a tendency to prolateness in the center, while they become more and more triaxial or oblate in the outer part, as a consequence of an accretion history that is more anisotropic at early times and becomes more isotropic at later times.

The results obtained from these DM-only simulations are affected by the absence of the baryonic component, whose inclusion is of fundamental importance for properly describing the formation of small-scale systems such as the MW and other galaxies. The inclusion of the gas dynamics in the simulations was a first attempt to account for the baryonic effects on the shape of the DM halo in the 1990s (e.g., Katz & Gunn 1991; Katz & White 1993; Dubinski 1994). In the following years, cosmological hydrodynamic simulations were used to investigate the impact of radiative cooling, star formation, and supernova feedback on the final shape of the DM halos (e.g., Gnedin et al. 2004; Kazantzidis et al. 2004; Gustafsson et al. 2006; Tissera et al. 2010; Abadi et al. 2010; Zemp et al. 2012; Bryan et al. 2013; Butsky et al. 2016). The inclusion of these effects leads to a larger sphericity of the DM halos, which are now predicted to be, in the central regions of galaxies, rounder or more oblate than previously thought.

In the last decades, predictions regarding the shape of the DM halos have been largely tested on the MW, with the aim of constraining the shape of its DM halo; those tests used different tracers of the Galactic gravitational potential, such as the distribution and kinematics of halo stars (e.g., Smith et al. 2009a; Loebman et al. 2014), the kinematics of disk stars in the solar neighborhood (Olling & Merrifield 2000), the tidal streams from satellite galaxies and from globular clusters (e.g., Helmi 2004; Johnston et al. 2005; Fellhauer et al. 2006; Růžička et al. 2007; Law et al. 2009; Koposov et al. 2010; Law & Majewski 2010; Vera-Ciro & Helmi 2013; Küpper et al. 2015; Bovy et al. 2016; Malhan & Ibata 2019), the distribution of both globular clusters (e.g., Posti & Helmi 2019) and the MW satellite galaxies (Zentner et al. 2005), and the flaring of the HI distribution (e.g., Olling & Merrifield 2000; Banerjee & Jog 2011).

Because of the use of different tracers, the results of these studies hold on different spatial scales. However, even on comparable scales, these results are not all consistent with one another, partly because of the use of different techniques and partly because of the different working assumptions.

The Galactic DM halo is suggested to be overall close to spherical on the basis of the tilt of the velocity ellipsoid of a sample of halo subdwarf stars located at galactocentric cylindrical radii of 7–10 kpc and depths ≲4.5 kpc below the Galactic plane, in the 250 deg² sky area covered by Stripe 82 of the Sloan Digital Sky Survey (Smith et al. 2009a). Under the assumptions that the MW DM halo is a spheroid and the full Galactic gravitational potential is axisymmetric, a variety of results are found. The GD-1 stellar stream excludes a significantly oblate DM halo at the GD-1 location, r ∼ 14 kpc, where the vertical-to-planar axis ratio (hereafter referred to as “flattening”) of the gravitational equipotential surfaces is constrained to q_Φ > 0.89, according to Koposov et al. (2010); the same stellar stream is found to provide a stronger constraint, yielding a prolate DM halo with mass density flattening q_ρ = 1.27, by Bovy et al. (2016). On the other hand, the stellar stream Pal 5 constrains the DM halo to be mildly oblate at r ∼ 19 kpc, with either potential flattening q_Φ = 0.95 (Küpper et al. 2015) or density flattening q_ρ = 0.9 (Bovy et al. 2016), suggesting a radial-dependent flattening for the DM halo. Finally, the combination of Pal 5 and GD-1, together with constraints on the force field near the Galactic disk, returns a nearly spherical DM halo, with q_ρ = 1.05, (Bovy et al. 2016) within the inner 20 kpc.

These results are at odds with those found with different probes on similar scales: the combination of the kinematics of disk stars in the vicinity of the Sun with the flaring of the HI disk is found to consistently constrain the DM halo to be oblate, with density flattening q_ρ ∼ 0.8 (Olling & Merrifield 2000); the kinematics of halo stars also constrains the dark halo to be significantly oblate, with flattening q_ρ = 0.4 (or, equivalently, q_Φ = 0.7; Loebman et al. 2014). Conversely, the distribution of globular clusters suggests a prolate DM halo with density flattening q_ρ = 1.3 (Posti & Helmi 2019).

At larger galactocentric distances, 20 ≲ r ≲ 60 kpc, the tidal streams of the Sagittarius dwarf spheroidal (Sgr dSph) galaxy lead to the conclusion that the DM halo potential has to be mildly oblate (q_Φ = 0.90−0.95; Johnston et al. 2005) or nearly spherical (q_Φ = 0.92−0.97; Fellhauer et al. 2006) to explain the precession of Sgr dSph’s orbit, while extremely oblate halos with density flattening q_ρ < 0.7 are ruled out (Ibata et al. 2001); on the other hand, the potential of the dark halo has to be prolate, with q_Φ = 1.25−1.5, to explain the kinematics of Sgr dSph’s older, leading stream (Helmi 2004). Finally, on scales r ≲ 200 kpc, the modeling of the Magellanic stellar streams generated by the interaction of the MW with the Magellanic system favors a DM halo that has a globally oblate potential with q_Φ < 1 (Růžička et al. 2007). The uncertainties on the quoted flattening range from ∼5% to ∼20%.

Some of the apparent inconsistencies among the results on the shape of the dark halo can be solved by assuming that the DM halo is triaxial and as such that the Galactic potential is globally non-axisymmetric. For instance, a triaxial DM potential with an intermediate-to-major axis ratio (b/a)_Φ = 0.99 and a minor-to-major axis ratio (c/a)_Φ = 0.72 on scales 20 ≲ r ≲ 60 kpc enabled Law & Majewski (2010) to explain, at the same time, both the angular precession and the radial velocities of the stars in the Sgr dSph leading stream (see also Law et al. 2009).

Other inconsistencies regarding the shape of the dark halo, especially those on different spatial scales, may in principle be resolved by discarding the simplifying assumption of a radial-independent shape for the halo: although common to all the abovementioned models, this assumption is not supported by the results of N-body simulations. In this context, by assuming an axisymmetric dark halo to model the flaring of the HI disk, Banerjee & Jog (2011) find the halo to be prolate, with a density flattening q_ρ increasing from 1 to 2 in the range 9 ≲ r ≲ 24 kpc. Conversely, with a non-axisymmetric model applied to the Sgr dSph streams, and accounting for the effects of the Large Magellanic Cloud (LMC), Vera-Ciro & Helmi (2013) constrain the DM halo potential to be mildly oblate at r ≲ 10 kpc, where q_Φ = 0.9, and smoothly translating to a triaxial shape at larger radii, where the intermediate-to-major axis ratio is (b/a)_Φ = 0.9 and the minor-to-major axis ratio is (c/a)_Φ = 0.8, in agreement with cosmological simulations. All these results clearly show that, despite significant efforts, the shape of the MW DM halo remains uncertain and that further work is necessary.

In this paper we use the hypervelocity stars (HVSs) as tracers of the shape of the MW DM halo. The existence of HVSs was postulated by Hills (1988), who defined as an HVS a star ejected from the Galactic Center after a close encounter between a stellar binary and the supermassive black hole (SMBH) associated with SgrA^⋆: the ejected star is characterized by a present-day speed exceeding the Galactic escape velocity, while the companion star becomes gravitationally bound to the SMBH. Alternative mechanisms for the generation of the HVSs were subsequently proposed, both before and after the first observational evidence for HVSs (Brown et al. 2005): a three-body interaction between a single star and a binary black hole (BH; Yu & Tremaine 2003), the scattering of a star off a stellar-mass BH (O’Leary & Loeb 2008), the star formation along the path of an outflow driven by an active galactic nucleus (Silk et al. 2012), the interaction between a globular cluster and a BH binary (Fragione & Capuzzo-Dolcetta 2016), and a four-body interaction between a binary star and a binary BH (Wang et al. 2018).

On the observational side, it was W. R. Brown who serendipitously discovered the first HVS candidate: a B-type star escaping the MW with a galactocentric velocity of ∼700 km s⁻¹ (Brown et al. 2005). Many HVS candidates were later found in both targeted and un-targeted surveys (e.g., Hirsch et al. 2005; Edelmann et al. 2005; Brown et al. 2006a,b, 2007a,b, 2009, 2012, 2014; Tillich et al. 2011; Li et al. 2012, 2015, 2021; Pereira et al. 2013; Zheng et al. 2014; Huang et al. 2017; Neugent et al. 2018; Du et al. 2019; Luna et al. 2019; Koposov et al. 2020). Some of these observed stars have speeds lower than the Galactic escape velocity: they are referred to as “bound HVSs”. As noted by Brown (2015), among the mechanisms proposed for HVS production, the Hills mechanism has the unique ability to generate a large number of unbound main-sequence stars and explain the presence of the so-called S stars in close orbit around the central SMBH (e.g., Ghez et al. 2003, 2005; Gillessen et al. 2009, 2017; Genzel et al. 2010).

We live in the Gaia era (e.g., Gaia Collaboration 2016a,b,c; Gaia Collaboration 2018a,b, 2021), and for the few last years several studies have been conducted to forecast the number of HVSs that Gaia will be able to detect by the end of the mission (Marchetti et al. 2018), to search for new HVSs candidates (e.g., Marchetti et al. 2017, 2019, 2021; Bromley et al. 2018), and to revalue the classification of a star as an HVS candidate based on the improved accuracy on proper motions provided by Gaia (e.g., Boubert et al. 2018; Brown et al. 2018; Irrgang et al. 2018, 2021; Erkal et al. 2019; Kreuzer et al. 2020). The birthplace of many stars is indeed still uncertain, and some of the HVS candidates may actually be: (i) runaway stars ejected from the Galactic disk (e.g., Blaauw 1961; Poveda et al. 1967; Leonard 1991; ii) halo star outliers whose velocities can reach 450 km s⁻¹ (Smith et al. 2009b) or fast halo stars, which can be unbound from the MW, produced by a tidal interaction between a dwarf galaxy and the MW near the Galactic Center (e.g., Abadi et al. 2009; Huang et al. 2021); or (iii) stars ejected as either HVSs or runaway stars from the nearest satellite galaxies of the MW, such as the LMC (Boubert & Evans 2016; Boubert et al. 2017) and the Sgr dSph galaxy (Huang et al. 2021). From the available literature, we estimate a sample of ∼70 HVS candidates for which a galactocentric origin has not been unambiguously ruled out. This sample includes both unbound and bound HVS candidates, which make up ∼40% and ∼60% of the full sample, respectively. These HVS candidates are the result of heterogeneous classification methods, and the number of true HVSs will remain uncertain until the galactocentric origin of these stars is unambiguously confirmed.

Since the discovery of the very first HVS candidate, HVSs have been recognized as a powerful tool for probing the shape of the Galactic DM halo (Gnedin et al. 2005; Yu & Madau 2007), its mass (Rossi et al. 2017; Fragione & Loeb 2017), or both (Contigiani et al. 2019). They have also been used to discriminate between different models of the Galactic gravitational potential in Newtonian gravity (Perets et al. 2009) and between different theories of gravity (Perets et al. 2009; Chakrabarty et al. 2022). In this paper, we assume that Newtonian gravity holds on galactic scales, we fix the mass of the DM halo, and we use the kinematical properties of the HVSs to constrain the shape of the dark halo. We account for neither the gravitational effects of the LMC on the HVS trajectories (Kenyon et al. 2018) nor the time dependence of the gravitational potential of the MW due to its interaction with the LMC (Boubert et al. 2020). Because the HVSs are ejected radially with high speed and may cross the entire Galaxy before dying out, in an isolated MW the small deviations from straight lines of their trajectories (typically of only a few per cent) are well determined by the asphericity of the MW gravitational potential, dominated by the DM halo at large galactocentric distances.

Previous studies aimed at constraining the shape of the gravitational potential of the dark halo made use of (i) sufficiently accurate (σ_μ ≲ 10 μas yr⁻¹) proper motion measurements of either one observed HVS with known distance or a set of two or more HVSs with unknown distance, with the constraints becoming tighter for larger samples of observed HVSs (Gnedin et al. 2005; ii) a triaxiality estimator that is a function of the components of the specific angular momentum of HVSs located at galactocentric distances r ≳ 50 kpc (Yu & Madau 2007); or (iii) a likelihood function constructed by back-propagating the phase-space position of each HVS to the Galactic Center in order to reproduce its observed phase-space coordinates and mass (Contigiani et al. 2019).

The above methods have been applied to triaxial (Gnedin et al. 2005; Yu & Madau 2007) and axisymmetric (Contigiani et al. 2019) DM halos. The techniques used by Gnedin et al. (2005) and Contigiani et al. (2019) do not depend on the model of the potential assumed for the DM halo, but require the integration of the HVS trajectories; conversely, the angular-momentum technique designed by Yu & Madau (2007) does depend on the dark halo potential but does not require the trajectory integration. In all these models, the shape of the DM halo potential is constant with radius. Furthermore, all the techniques require intermediate steps that involve each HVS of the sample individually, even though the final result depends on the contribution of all the HVSs of the sample.

In this work we propose a statistical method for constraining the axis ratios of a radial-dependent, triaxial gravitational potential of the DM halo from the distributions of the HVS phase-space coordinates that are mostly affected by the asphericity of this potential, namely the components of the galactocentric tangential velocities. Unlike the techniques illustrated above, our method does not require intermediate steps that involve each star of the sample, such as the trajectory integration or the evaluation of the angular momenta of the HVSs. Our method can be applied to different models of the MW gravitational potential, as is the case for the methods by Gnedin et al. (2005) and Contigiani et al. (2019).

The paper is organized as follows. In Sect. 2 we describe our numerical simulations of the initial velocity distribution of a sample of HVSs ejected according to the Hills mechanism and the simulations of the HVS trajectories in a Galactic gravitational potential generated by DM halos with different shapes; we also illustrate the construction of our HVS phase-space mock catalogs. In Sect. 3 we show how the asphericity of the DM halo mostly affects the HVS tangential velocity: we identify this velocity as the key variable for statistically discriminating between different shapes of the DM halo, and we select the appropriate HVS sample to pursue this goal. In Sect. 4 we present our statistical method for recovering the shape of the DM halo from a distribution of HVS tangential velocities. In Sects. 5 and 6 we show the results of the application of our method to an ideal sample of mock observed HVSs with null uncertainties and no observational limitations that traveled in an axisymmetric and non-axisymmetric Galactic gravitational potential, respectively. In Sect. 7 we investigate the effect of the size of the ideal sample of mock observed HVSs on the success rate of our method. We finally discuss our results and conclude in Sects. 8 and 9, respectively.

2. Numerical simulations and mock catalogs

In Sect. 2.1 we illustrate the generation of the distribution of the initial velocities of HVSs ejected with the Hills mechanism. In Sect. 2.2 we describe our model for the gravitational potential of the MW, generated by the baryonic components and a DM halo with different shapes, and in Sect. 2.3 we illustrate the simulation of the trajectories of the ejected HVSs across the Galaxy. Finally, in Sect. 2.4 we describe the mock catalogs that we built from the HVS phase-space distribution.

2.1. Ejected stars: velocity distribution

Following Hills (1988) (see also Bromley et al. 2006), we simulated the ejection of stars from the Galactic Center with a three-body numerical code that reproduces the gravitational interaction of a binary star system with the SMBH associated with SgrA^⋆. For the ejected stars, our code provides the distribution of the ejection velocities, v_ej. A detailed description of the code will be provided in a separate paper.

We set the mass of the SMBH to 4 × 10⁶ M_⊙, consistent with different estimates (e.g., Boehle et al. 2016; Gillessen et al. 2017). For simplicity, we restricted our simulations to equal-mass stellar binaries of 4 + 4 M_⊙ on hyperbolic orbits. The mass of 4 M_⊙ is representative of the upper end of the mass distribution of the HVS candidates observed in currently available surveys (see references in Sect. 1). We further discuss our mass choice in Sect. 8.

For a fixed mass of the binary members, the velocity distribution of the ejected stars depends on a series of parameters: (i) the stellar binary semimajor axis, a; (ii) the minimum approach distance, R_min, between the center of mass of the binary and the SMBH; (iii) the inclination angle, i, between the orbital plane of the binary star and the orbital plane of the binary’s center of mass and the SMBH; and (iv) the initial phase, ϕ, of the binary star. We randomly sample i and ϕ in the interval [0, 2π] with uniform probability density function. As for a and R_min, which determine the probability of ejection of the primary star, we randomly sampled a from the interval 0.05−4 AU with probability density function p(a)∝1/a, and we randomly drew R_min from the interval 1−700 AU with probability density function p(R_min)∝R_min (Bromley et al. 2006, and references therein).

Figure 1 shows the ejection velocity distribution for a simulation of N_int ∼ 240 000 three-body interactions that produces N_ej = 60 000 ejection events. Most of the ejected stars possess ejection velocities in the range ∼250−4000 km s⁻¹; the velocity distribution has a major peak at v_ej ∼ 510 km s⁻¹ and has positive skewness. A minority of stars are ejected with speeds lower than ∼250 km s⁻¹: these stars are the result of rare, double-ejection events, where neither of the binary members becomes gravitationally bound to the SMBH, and both stars are ejected. Because we assumed that the binary system starts its approach trajectory toward the SMBH with a velocity at infinity v_inf = 250 km s⁻¹ (Hills 1988), energy conservation requires that one of these stars is ejected with velocity lower than 250 km s⁻¹. Our distribution of ejection velocities is comparable to the ejection speed distribution obtained from the analytical prescription by Bromley et al. (2006).

Fig. 1. Distribution of ejection velocities for the ejected star(s) of a 4 + 4 M⊙ binary star system, after a close encounter with a 4 × 106 M⊙ SMBH. The distribution consists of Nej = 60 000 stars ejected in Nint ∼ 240 000 three-body interactions.

We note that the ejection velocity v_ej is defined as the speed that an ejected star would have at infinite distance from the SMBH, in absence of other gravitational sources (Hills 1988; Bromley et al. 2006), as in our three-body simulation. However, the ejection velocity distribution shown in Fig. 1 can be used as a distribution of initial velocities when simulating the trajectories of the HVSs from the Galactic Center to the outer halo, as we explain in Sect. 2.3.

2.2. Milky Way gravitational potential

We modeled the MW gravitational potential Φ as the superposition of the potentials generated by three distributions of baryonic matter and one distribution of DM:

$$\begin{array}{c}\hfill \mathrm{\Phi }={\mathrm{\Phi }}_{\mathrm{BH}}+{\mathrm{\Phi }}_{\mathrm{b}}+{\mathrm{\Phi }}_{\mathrm{d}}+{\mathrm{\Phi }}_{\mathrm{h}},\end{array}$$(1)

where Φ_BH is the potential generated by the SMBH located in the Galactic Center, Φ_b is generated by the Galactic bulge, Φ_d is the disk potential, and Φ_h is the potential of the DM halo embedding the Galaxy. We consider the Galaxy as isolated, neglecting both the presence of the LMC (Kenyon et al. 2018) and any time dependence of the Galactic potential deriving from the interaction of the MW with the LMC (Boubert et al. 2020).

In a reference frame with the origin at the Galactic Center, we used spherical coordinates (r, ϑ, φ) for the spherically symmetric components of the potential (with −90° ≤ϑ ≤ 90°), cylindrical coordinates (R, φ, z) for the axisymmetric components, and Cartesian coordinates (x, y, z) for the triaxial component. We took the x − y plane as the equatorial plane of the disk, with the x axis corresponding to the direction from the Sun to the Galactic Center, and we took the z axis as the vertical axis.

We included the contribution of the SMBH to the gravitational potential as

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{BH}}(r)=-\frac{G{M}_{\mathrm{BH}}}{r},\end{array}$$(2)

where M_BH is the mass of the SMBH.

For the bulge component, we adopted the Hernquist (1990) potential:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{b}}(r)=-\frac{G{M}_{\mathrm{b}}}{r+r_{\mathrm{b}}},\end{array}$$(3)

where M_b = 3.4 × 10¹⁰ M_⊙ and r_b = 0.7 kpc (Kafle et al. 2014; Price-Whelan et al. 2014; Rossi et al. 2017; Contigiani et al. 2019) are the scale mass and the scale radius of the model, respectively.

For the disk component, we adopted the axisymmetric potential by Miyamoto & Nagai (1975):

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{d}}(R,z)=-\frac{G{M}_{\mathrm{d}}}{\sqrt{R^2+{(a_{\mathrm{d}}+\sqrt{z^2+b_{\mathrm{d}}^2})}^2}},\end{array}$$(4)

where M_d = 1.0 × 10¹¹ M_⊙, a_d = 6.5 kpc and b_d = 0.26 kpc (Kafle et al. 2014; Price-Whelan et al. 2014; Rossi et al. 2017; Contigiani et al. 2019) are the scale mass and the scale lengths of the model, respectively.

Finally, we modeled the contribution of the DM halo with the triaxial generalization of the spherically symmetric NFW (Navarro et al. 1997) potential proposed by Vogelsberger et al. (2008):

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{h}}(\stackrel{\sim }{r})=-\frac{G{M}_{200}}{f(C_{200})}\frac{\ln (1+\frac{\stackrel{\sim }{r}}{r_{\mathrm{s}}})}{\stackrel{\sim }{r}},\end{array}$$(5)

where f(w) = ln(1 + w)−w/(1 + w), M₂₀₀ = 8.35 × 10¹¹ M_⊙ is the DM halo mass enclosed within r₂₀₀¹, C₂₀₀ = r₂₀₀/r_s = 10.82 is the halo concentration parameter, and r_s = 18 kpc is a generalized scale radius. For the above parameters, we adopted the values used by Hesp & Helmi (2018), which are consistent with the values derived from halo stars (Xue et al. 2008), blue horizontal branch stars (Deason et al. 2012), the massive satellite population of the MW (Cautun et al. 2014), and Cepheids (Ablimit et al. 2020). The coordinate $\stackrel{\sim }{r}$ is a generalized radius that replaces the radius r of the NFW spherical potential,

$$\begin{array}{c}\hfill \stackrel{\sim }{r}=\frac{(r_{\mathrm{a}}+r)r_{\mathrm{E}}}{r_{\mathrm{a}}+r_{\mathrm{E}}}.\end{array}$$(6)

Here, r_E is an “ellipsoidal radius”,

$$\begin{array}{c}\hfill r_{\mathrm{E}}=\sqrt{\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}},\end{array}$$(7)

where the three ellipsoid parameters, a, b, and c have to fulfill the condition a² + b² + c² = 3, and their combination defines the degree of triaxiality of the potential well. Specifically, the axis ratio of the equipotential surfaces on the x–y plane is q_y = b/a, whereas the axis ratio of the equipotential surfaces on the x–z plane is q_z = c/a. The parameter r_a is the scale length where the smooth transition from a triaxial potential to a nearly spherical potential occurs; we took it to be 1.2r_s, as in Hesp & Helmi (2018): the halo is triaxial ($\stackrel{\sim }{r}\approx r_E$) in the inner region (r ≪ r_a), whereas it is approximately spherical ($\stackrel{\sim }{r}\approx r$) in the outer region (r ≫ r_a).

In general, the potential ${\mathrm{\Phi }}_{\mathrm{h}}(\stackrel{\sim }{r})$ given by Eq. (5) is triaxial with q_y ≠ 1, q_z ≠ 1, and q_y ≠ q_z. However, this potential becomes spheroidal when either q_y or q_z is equal to 1 or when q_y = q_z ≠ 1: when q_y = 1, the potential is axisymmetric about the z axis; when q_z = 1, the potential is axisymmetric about the y axis; when q_y = q_z ≠ 1 the potential is axisymmetric about the x axis. When both q_y and q_z are equal to 1, the DM halo potential is spherically symmetric.

Because Φ_BH and Φ_b are spherically symmetric, and Φ_d is axisymmetric about the z axis, our gravitational potential Φ (Eq. (1)) is globally axisymmetric about the z axis when the DM halo is either spherical or spheroidal and axisymmetric about the z axis. On the other hand, the Galactic gravitational potential is non-axisymmetric when the DM halo is either triaxial or spheroidal with a symmetry axis misaligned with respect to the z axis.

In all of our simulations, we adopted the same parameters for the SMBH, the bulge, and the disk potentials (Eqs. (2)–(4)). We also adopted the same parameters for the potential of the DM halo (Eq. (5)), with the exception of (i) the triaxiality parameter q_z, which was set to a different value in each of the simulations of the axisymmetric Galactic potential, while q_y was kept fixed to 1, and (ii) both triaxiality parameters, q_y and q_z, in simulations of a non-axisymmetric Galactic potential. By varying the triaxiality parameter of the DM halo q_y and q_z in appropriate ranges of values, we explored the effect of the halo shape on the HVS observables.

Figure 2 shows the magnitude of the radial gravitational acceleration in the plane of the disk associated with our MW potential Φ as a function of the cylindrical coordinate R, for a spherical DM halo (q_y = 1, q_z = 1). The chosen gravitational potential generates masses enclosed within 120 pc and within 100 kpc that agree with the observed values derived by Launhardt et al. (2002) and reported in Dehnen & Binney (1998), respectively. It also reproduces a circular velocity of 235 km s⁻¹ at the solar neighborhood, in agreement with the observational values of 231 ± 6 km s⁻¹ by Bobylev (2017) and of 230 ± 12 km s⁻¹ by Bobylev & Bajkova (2016).

Fig. 2. Magnitude of the radial gravitational acceleration in the plane of the disk due to our MW potential Φ (Eq. (1)), as a function of the radial cylindrical coordinate R. The solid black line represents the total gravitational acceleration. The dashed and dotted lines represent the contributions of the SMBH (dashed orange line), the bulge (dotted green line), the disk (dot-dashed blue line), and a spherical DM halo (long-dashed magenta line).

2.3. Orbit integration

In the Galactic gravitational potential illustrated in Sect. 2.2, we simulated the time evolution of the position, velocity, and acceleration of a sample of stars ejected from the Galactic Center according to the Hills mechanism (see Sect. 2.1). For each star, we numerically integrated Newton’s equation of motion in a galactocentric reference frame and in Cartesian coordinates, using the Velocity Verlet algorithm (e.g., Frenkel & Smit 2001). We traced the star trajectory until the star death. In our simulations, the total energy of a star is conserved with a relative accuracy of ∼10⁻⁸.

In Sect. 2.1, we showed that our sample of ejected star consists of N_ej = 60 000 HVSs; we adopted this sample as our full sample of initial velocity magnitudes. This choice is legitimate: even though v_ej is defined as the speed that an ejected star would have at infinite distance from the SMBH in a three-body interaction (see Sect. 2.1), our SMBH is embedded in the Galactic mass distribution; thus, the ejected stars can be assumed to move at v_ej at the outer edge of the sphere of influence of the SMBH, just before the Galactic gravitational potential starts to overcome the SMBH potential. The Hills ejection mechanism yields an isotropic distribution of ejected stars. We thus assigned each of the N_ej stars an initial position (r₀, ϑ, φ), with r₀ = 3 pc the radius of the sphere of influence of the SMBH (Genzel et al. 2010), and (ϑ, φ) randomly drawn from a uniform distribution over the surface of a sphere; we assigned each star a velocity direction $\widehat{n}(\vartheta ,\phi )$. We end up with a sample of N_ej initial conditions.

In our simulations below, we computed the trajectories of the HVSs in a large number of different subsamples of size smaller than N_ej. Each subsample was generated by adopting a random subsample of the N_ej initial conditions. Each subsample was selected as follows.

At an average rate R = 10⁻⁴ yr⁻¹, roughly consistent with the estimate by Bromley et al. (2012) and Zhang et al. (2013) (see also Hills 1988; Yu & Tremaine 2003), all the N_ej HVSs would be ejected in 600 Myr. We thus assigned the i-th HVS of each subsample an ejection time t_ej, i uniformly sampled in the range 0−600 Myr. The stars have a finite lifetime and each star can thus have a different residual lifetime at t_ej, i. A 4 M_⊙ star with solar metallicity has a main-sequence lifetime τ_ms ≃ 160 Myr (Schaller et al. 1992; Brown et al. 2006b): we took this lifetime as the total lifetime of the star, τ_L. Thus, at the time of ejection, the i-th star is also assigned an age τ_ej, i randomly sampled from a uniform distribution between 0 and τ_L.

Figure 3 shows that the number N(t) of observable stars, whose lifetime is τ_L, reaches a steady state N(t)≃N_s ≃ 8000 at t ≃ τ_L. The steady state is the result of the combination of (i) a continuous ejection of stars, with average rate R, and (ii) the finite lifetime of the stars, τ_L. We chose to observe our star sample at the observation time, t_obs = 400 Myr. All the stars whose ejection time was larger than t_obs were discarded from the sample. Among the remaining HVSs, we selected the N_obs stars that are alive at the observation time t_obs, namely the stars that satisfy the condition t_obs − t_ej, i + τ_ej, i < τ_L. For these stars, we computed the trajectory through the Galaxy for a travel time t_travel, i = t_obs − t_ej, i.

Fig. 3. Number of living 4 M⊙ HVSs as a function of the time of observation. The ejection begins at t = 0 with a rate 10−4 yr−1. The number of living HVSs stops increasing beyond t ≃ τL = 160 Myr. The vertical dashed line marks the time of the steady state when we chose to observe the system, tobs = 400 Myr.

We computed the orbits of the sample of N_obs ≃ N_s ≃ 8000 HVSs in a series of 64 simulations, each of them corresponding to a different combination (q_y, q_z) of the triaxiality parameters of the DM halo (see Sect. 2.2). The different combinations were obtained by varying both q_y and q_z in the range 0.7−1.4 in steps of 0.1. A summary of the simulated shapes of the DM halo is reported in Table 1.

2.4. Mock catalogs A and B

With the simulations described in Sect. 2.3, we created two series of mock catalogs, hereafter referred to as mock catalogs A and B. In both series, each mock catalog includes the phase-space distribution of N_obs HVSs at the observation time t_obs in a DM halo with a given shape. However, the two series differ in the initial conditions of the HVSs.

Each mock catalog A was generated with the same set of N_obs initial conditions. In other words, all the stars with the same identification index i in all catalogs A are given the same combination Sⁱ={v_ej, $\widehat{n}(\vartheta ,\phi )$, t_ej, τ_ej}ⁱ of initial conditions, namely the same ejection velocity, ejection time t_ej, and age τ_ej at ejection. The only difference among these mock catalogs is the set of triaxiality parameters of the DM halo, (q_y, q_z). We used mock catalogs A to highlight the effect of the degree of triaxiality of the DM halo on the distribution of the phase-space coordinates of the HVSs, and to identify the phase-space coordinates that can serve as shape indicators (Sect. 3.1).

On the other hand, mock catalogs B were generated with different sets of N_obs initial conditions. In other words, the i-th stars in different catalogs B are given a different combination Sⁱ of initial conditions. Therefore, mock catalogs B differ from one another not only for the shape of the DM halo, but also for the sample of ejected stars.

We used mock catalogs B to explore the effect of the variation of the initial conditions on the detection of deviations from spherical symmetry of the shape of the DM halo potential (Sect. 3.3). We also used catalogs B to implement our method to recover the shape of the DM halo of the MW from real HVS samples (Sect. 4).

Each mock catalog A or B includes the positions and the velocities of the N_obs bound and unbound HVSs in the galactocentric, the Galactic heliocentric, and the equatorial reference frames, at the observation time t = t_obs. In our simulations, the star position in the galactocentric reference frame is (r, ϑ, φ) in spherical coordinates and (x, y, z) in Cartesian coordinates, with φ = 0° on the x axis, and ϑ = 0° on the x–y plane; the star galactocentric velocity is (v_r, v_ϑ, v_φ), with v_r the radial velocity, v_ϑ the latitudinal velocity, and v_φ the azimuthal velocity. We defined as tangential velocity the vector whose components are the latitudinal and the azimuthal velocities: v_t = (v_ϑ, v_φ); its magnitude is $v_{\mathrm{t}}\equiv |{\mathit{v}}_{\mathrm{t}}|={(v_{\vartheta }^2+v_{\phi }^2)}^{\frac{1}{2}}$. We show the distribution of the galactocentric distance, as well as the distributions of the galactocentric radial, latitudinal, and azimuthal velocities of the HVSs for one of our mock catalogs in Appendix A.

As will be shown at the end of Sect. 3.1, we evaluated the possible use of galactocentric, heliocentric, and equatorial phase-space coordinates to infer the shape of the dark halo. To this aim, we converted the star phase-space coordinates from the galactocentric to the Galactic heliocentric reference frame, by choosing the Sun to be located at (x, y, z) = (−R_⊙, 0, 0) and to have velocity (U_⊙, V_⊙ + Θ₀, W_⊙) in the galactocentric reference frame. We used R_⊙ = 8.277 kpc (Abuter et al. 2021), and our model’s rotational velocity at R_⊙, Θ₀ = 235 km s⁻¹, in agreement with Bobylev (2017) and Bobylev & Bajkova (2016) (see Sect. 2.2). For the velocity of the Sun with respect to the local standard of rest, we used U_⊙ = 11.1 km s⁻¹, V_⊙ = 12.24 km s⁻¹, and W_⊙ = 7.25 km s⁻¹ (Schönrich et al. 2010). The equations for the transformation of coordinates and velocities from the galactocentric to the Galactic heliocentric reference frame, and from the Galactic heliocentric to the equatorial reference frame, are reported in Appendix B.

3. Indicators of the shape of the DM halo

We used mock catalogs A described in Sect. 2.4 to investigate the effect of the triaxiality parameters of the gravitational potential of the DM halo on the HVS phase-space distribution. Our final goal was to identify the HVS phase-space coordinates (Sect. 3.1) and define the HVS samples (Sects. 3.1 and 3.2) that are best suited to discriminate between different shapes of the DM halo.

3.1. Effect of the shape of the DM halo on the HVS phase-space distribution

We explored the effect of the triaxiality parameters of the gravitational potential of the DM halo on the phase-space distribution of the simulated HVSs, with the goal of identifying the phase-space coordinates that are most sensitive to the shape of the DM halo, and can thus be identified as triaxiality indicators. For this purpose, we resorted to mock catalogs A: because these catalogs are all generated with a specific sample of stars, characterized by fixed distributions of initial conditions (Sect. 2.4), any significant difference among two of these mock catalogs can only be ascribed to the different triaxiality parameters of the gravitational potential of the DM halo.

For our investigation, from each mock catalog A of N_obs stars (see Sect. 2.4), we selected a subsample of stars that, at the time of observation t_obs, are located at galactocentric distances r > 10 kpc, where the gravitational effect of the DM halo starts to be relevant (see Fig. 2). We also required the stars to possess positive galactocentric radial velocity, v_r > 0, to match the observational HVS selection criterion. Applying these selection criteria returns, for each mock catalog, a subsample of N ≃ N_obs/10 ≃ 800 stars. In Sect. 3.2 we demonstrate that this reasonable sample selection turns out to be the best selection on the basis of kinematic arguments.

We performed our analysis by means of a statistical approach. Specifically, we made use of the two sample Kolmogorov-Smirnov test (Press et al. 2007), hereafter referred to as “KS test”, to check whether the null hypothesis H₀ is true, namely whether the distribution of a given phase-space coordinate in a spherical DM halo and the distribution of the same coordinate in a nonspherical DM halo are consistent with being drawn from the same parent population, and are thus indistinguishable. If the p-value of the KS test is p ≤ α, then H₀ is rejected at the adopted significance level α: in this case, we considered the two distributions as significantly different from each other, and we investigated the use of that phase-space coordinate as an indicator of the shape of the DM halo. We adopted a significance level α = 5%.

For our KS tests, we considered the distributions of the components of the star position and velocity in the galactocentric reference frame, in both spherical and Cartesian coordinates, for a series of pairs of DM halos; each halo pair is composed of a spherical DM halo and a nonspherical DM halo with different triaxiality parameters. For the sake of simplicity, we illustrate here the details of our investigation for spheroidal DM halos with q_y = 1, namely for spheroidal halos that are axisymmetric about the z axis and yield a global axisymmetric Galactic potential; we only comment on the case of a non-axisymmetric Galactic potential. This restriction will not imply a loss of generality of our main result. Both the axisymmetric and the non-axisymmetric cases are carefully explored in Sects. 5 and 6.

We created 50 series of mock catalogs A. Each series consisted of a spherical DM halo and six spheroidal DM halos axisymmetric about the z axis and with q_z ranging from q_z = 0.7 (extremely oblate halo) to q_z = 1.4 (extremely prolate halo). The set of initial conditions is the same for each halo in the same series, and differs from one series to the other. We computed the p-values of the KS tests between the distributions of the phase-space coordinates in the spherical halo and those in the spheroidal halos of the same series. Table 2 shows the range of these p-values for the 50 series of mock catalogs.

The distributions of the magnitude of the latitudinal velocity, |v_ϑ|², are the only distributions to be significantly different in a spherical DM halo and in a spheroidal DM halo axisymmetric about the z axis: for DM halos whose gravitational potential well displays a deviation from the spherical shape |q_z − 1|≥0.1, we always got p < α, implying that the null-hypothesis of the KS test could be rejected at the significance level α. In particular, for the largest deviations from the spherical DM halo considered in this work, the p-value is < 10⁻¹⁰. On the other hand, for a very mild deviation of the potential well of the DM halo from the spherical shape (i.e., for |q_z − 1|≤0.05) the distributions of |v_ϑ| are never significantly different from the spherical case. Finally, for DM halos whose potential well displays a deviation from the spherical shape |q_z − 1| in the range (0.05−0.1), the distributions of |v_ϑ| can either be or not be significantly different from the spherical case, depending on the set of stars’ initial conditions used for the generation of mock catalogs A. For example, when |q_z − 1| = 0.075, we found p in the range 0.01−0.14 for the case of spherical versus oblate (with q_z = 0.925) DM halo, and p in the range 0.01−0.11 for the case of spherical versus prolate (with q_z = 1.075) DM halo.

For all the distributions of phase-space coordinates different from |v_ϑ|, we found p > 0.7: we could not reject the null-hypothesis of the KS test at the significance level α = 5%, regardless of the shape of the DM halo. Therefore, we conclude that the magnitude of the HVS latitudinal velocity, |v_ϑ|, is the only phase-space coordinate whose distribution can be used to discriminate between a spherical and a spheroidal DM halo axisymmetric about the z axis.

This result is not surprising. HVSs are ejected radially outward from the Galactic Center, but they attain nonzero tangential velocities, v_t, as they travel through a non-spherically symmetric potential. In our model (see Eq. (1)), the potentials of both the SMBH and the bulge are spherically symmetric (see Eqs. (2)–(3)). However, the disk potential (see Eq. (4)) is axially symmetric about the z axis. Hence, it contributes to the component of the tangential velocity along the polar angle, that is, the latitudinal velocity v_ϑ. When the DM halo is either spherical or spheroidal with axial symmetry about the z axis, v_ϑ is still the only non-null component of the tangential velocity. In particular, when the DM halo is spherical, only the gravitational pull of the disk induces non-null v_ϑ; on the other hand, when the DM halo is spheroidal with axial symmetry about the z axis, the gravitational pull of the disk combines with that of the DM halo: in the case of prolate spheroidal DM halo (q_z > 1) the pull of the disk is opposite to that exerted by the DM halo, while in the case of oblate spheroidal DM halo (q_z < 1) both the disk and the halo attract the stars toward the Galactic plane, leading to higher tangential velocities.

This situation is illustrated in Fig. 4, where we show a comparison of the distributions of the magnitude of the latitudinal velocity, |v_ϑ|, for simulated HVSs in a Galaxy with a spheroidal DM halo that is axisymmetric about the z axis and characterized by an extremely prolate (q_z = 1.4), spherical (q_z = 1), and extremely oblate (q_z = 0.7) shape. The fraction of HVSs with higher |v_ϑ| is larger in the case of the oblate DM halo (green histogram) than in the case of the spherical DM halo (gray, shaded histogram), because of the concordant gravitational pull of disk and DM halo. Conversely, the fraction of stars with lower |v_ϑ| is larger in the case of the prolate DM halo (blue histogram) than in the case of the spherical DM halo, because of the opposite pull of disk and DM halo.

Fig. 4. Distribution of the magnitude of the HVS galactocentric latitudinal velocity, |vϑ|, in a Galaxy with an extremely prolate DM halo (qz = 1.4; blue histogram), a spherical DM halo (qz = 1; gray, shaded histogram), and an extremely oblate DM halo (qz = 0.7; green histogram). The distributions were generated with the same initial conditions (mock catalogs A) to highlight the effect of the different geometries of the DM halo.

The results reported in Table 2 show that, while the latitudinal velocity v_ϑ is an indicator of the shape of the DM halo, the polar angle ϑ is not, even though any difference in the distributions of ϑ for stars that have traveled in different DM halos is only determined by the halo shape, as v_ϑ is. This higher sensitivity of v_ϑ to the changes of the shape of the DM halo has two explanations: (i) v_ϑ is the time derivative of ϑ; any variation in v_ϑ implies a variation Δϑ of the coordinate ϑ over a finite time interval Δt; however, significant Δϑ can be achieved only over Δt that, on average, are larger than the HVS travel time; (ii) the v_ϑ’s are null at the ejection, and the final distribution of v_ϑ is determined by the gravitational potential alone; instead, the distribution of ϑ is uniform in cos ϑ at the start: hence, its final distribution is the result of the combination of the initial randomness and of the action of the gravitational potential.

Table 2 also shows that neither the Cartesian spatial coordinates nor the Cartesian components of the HVS velocity are significantly affected by the deviation from spherical of the shape of the DM halo; therefore, they are not useful indicators of the triaxiality of the DM halo. This result is expected for the Cartesian spatial coordinates x and y, as well as for the velocity components v_x and v_y, because the gravitational potential is symmetric about the z axis. On the other hand, for the same reason, the result is not fully expected for z and v_z. However, our simulations show that the projection, on the z axis, of any non-null v_ϑ induced by a nonspherical gravitational potential in the star motion is overwhelmed by the projection, on the same axis, of the radial velocity component of the star. If v_z is not a good indicator of the triaxiality of the DM halo, the spatial coordinate z cannot be a good indicator either: as ϑ is less sensitive than v_ϑ to the shape of the DM halo, z is less sensitive to the halo shape than v_z.

The results that we obtained for a spheroidal DM halo symmetric about the z axis (i.e., for a DM halo with q_y = 1 and q_z ≠ 1) can be extended to the more complex cases of a fully triaxial DM halo and of a spheroidal DM halo with a symmetry axis misaligned with respect to the z axis (i.e., for a DM halo with q_y ≠ 1). In those cases, the stars acquire both a non-null latitudinal velocity, v_ϑ, and a non-null azimuthal velocity, v_φ. Therefore, both components of v_t = (v_ϑ, v_φ) can be used as indicators of the triaxiality parameters of the DM halo.

We note that significant differences between the distributions of HVS tangential velocity components in a spherical and in a nonspherical DM halo only emerge when those velocity components are computed in the galactocentric reference frame. In the Galactic heliocentric reference frame, no phase-space coordinate is characterized by distributions that significantly differ in the cases of spherical and nonspherical DM halos, according to the KS test. Specifically, this is the case for each of the components of the star velocity. Indeed, all the components of the heliocentric velocity (v_d, v_l, v_b) are a composition of v_r and v_t = (v_ϑ, v_φ) (see Eqs. (B.4)–(B.6)), and the information on the triaxiality parameters stored in the galactocentric v_t’s is diluted in the velocity transformation from the galactocentric to the heliocentric system. The same results hold for the phase-space coordinates in the equatorial reference frame.

3.2. Star kinematics and sample selection

As illustrated in Sect. 3.1, the components v_ϑ and v_φ of the tangential velocity of the HVSs can effectively probe the nonspherical components of the gravitational potential. This result was demonstrated for a subsample of mock HVSs characterized by r > 10 kpc and v_r > 0. Here, we show that these sample selection criteria turn out to be the most suitable selection criteria also on the basis of the star kinematics, for HVSs of 4 M_⊙. Indeed, these criteria enable us to select those stars whose tangential velocity is not affected by effects other than the shape of the gravitational potential well. We adopted an analogous HVS sample selection based on stellar kinematics in Chakrabarty et al. (2022).

In our model of the Galactic potential (see Sect. 2.2), stars with ejection speed v_ej ≳ 800 km s⁻¹ always possess tangential velocities that are independent of their radial velocities, and are induced only by the nonspherical shape of the gravitational potential well. Those stars are robust indicators of the shape of the DM halo in any stage of their trajectory. Conversely, for stars with ejection speed v_ej ≲ 800 km s⁻¹, the tangential velocity, v_t, can be strongly coupled with the radial velocity, v_r, for a significant fraction of their trajectory, and v_t can be very high regardless of the shape of the DM halo. Indeed, if any of these stars is sufficiently young at ejection, it may experience the outer turnaround before dying out. At the outer turnaround, the star starts falling back toward the Galactic Center, its radial velocity becomes negative, and its tangential velocity starts growing; this growth continues until the star experiences the inner turnaround, and then starts a new outward trajectory, with positive radial velocity and a tangential velocity that decreases with time. Therefore, those stars may serve as indicators of the halo shape only in specific parts of their outward trajectory, namely only during specific periods of their lifetime.

The situation is illustrated in Fig. 5 for the case of a spherical DM halo, and for a set of HVSs ejected radially outward in the direction $\widehat{n}(\vartheta ,\phi )=\widehat{n}(45°,45°)$ and with representative values of the ejection speed, |v_ej|={700, 720, 740, 800, 900} km s⁻¹. All the stars are ejected with null age, τ_ej = 0, to illustrate the evolution of the stars’ observables during the largest possible travel time (i.e., 160 Myr, for the 4 M_⊙ stars considered in this work). Stars not ejected with null age would experience only part of the evolution shown in the figure. We note that the choice of the ejection direction of the HVSs does not affect any of the results that are presented in the following.

Fig. 5. Kinematics of HVSs with different ejection speeds. Panel a: magnitude of the galactocentric latitudinal velocity, |vϑ|, as a function of the galactocentric radial velocity, vr, for HVSs ejected with different speeds, vej, and traveling in a spherical DM halo. Initially, the star velocity is purely radial, and |vϑ| = 0; as time goes on, |vϑ| becomes non-null due to the nonspherical potential of the Galactic disk (see Eq. (4)). For stars with lower ejection speeds (vej ≲ 800 km s−1), radial and angular dynamics are coupled: this coupling results in a fast growth of |vϑ| after the first turnaround of the star. However, such couplings are not manifested for stars with higher ejection speeds. Panel b: galactocentric radial velocity, vr, as a function of the galactocentric distance, r, for HVSs ejected radially outward with different ejection velocities, vej. The vertical dashed line of panel a and the horizontal dashed line of panel b correspond to vr = 0, while the vertical dashed line of panel b corresponds to r = 10 kpc. In both panels, HVSs are ejected in the direction $\widehat{n}(\vartheta ,\phi )=\widehat{n}(45°,45°)$ with a travel time of 160 Myr, which is the largest possible travel time for a 4 M⊙ star.

Panel a of Fig. 5 shows the relation between the magnitude of the tangential velocity, v_t = |v_ϑ|, and the radial velocity, v_r, for the above set of stars. We note that the magnitude of the tangential velocity, $v_{\mathrm{t}}={(v_{\vartheta }^2+v_{\phi }^2)}^{\frac{1}{2}}$ is equivalent to the magnitude of the latitudinal velocity |v_ϑ| because the DM halo is spherical: since the axisymmetric disk potential is the only source of nonzero v_t, the azimuthal velocity v_φ is null.

Panel b of Fig. 5 shows the relation between the radial velocity, v_r, and the distance to the Galactic Center, r, for the same set of HVSs. Both the vertical dashed line in panel a and the horizontal dashed line in panel b correspond to a null radial velocity, v_r = 0, that the star possesses at its outer and inner turnaround radii.

For the HVSs generated with ejection speeds v_ej ≳ 800 km s⁻¹ (blue and purple lines), both panels a and b show that these stars can never reach their outer turnaround radius (r = r_out) before dying out. Indeed, they never cross the line v_r = 0 from right to left in panel a, and the line v_r = 0 from top to bottom in panel b.

For these stars, panel a shows that |v_ϑ| is always less than or around a few tens of km s⁻¹, regardless of v_r: these values of v_ϑ are entirely determined by the deviation from the spherical symmetry of the gravitational potential. This deviation is due to the axisymmetric disk alone, in the case of the spherical DM halo considered here; it is due to a combination of disk and triaxial DM halo when the DM halo is nonspherical, as discussed in Sect. 3.1. At any stage of the stars’ trajectories, the distributions of |v_ϑ| for these stars can be used as an indicator of the shape of the DM halo.

On the contrary, for the HVSs with ejection speed v_ej ≲ 800 km s⁻¹ (cyan, dark green, green, orange, and red lines), both panels a and b show that these stars may undergo at least the first, outer turnaround before dying out. Indeed, these stars may cross both the line v_r = 0 from right to left in panel a, and the line v_r = 0 from top to bottom in panel b at the outer turnaround radius, r = r_out, if they are ejected with sufficiently low age. When v_ej ≲ 740 km s⁻¹, the stars may also experience the second, inner turnaround. In this case, they cross again the line v_r = 0 from left to right in panel a, and correspondingly cross the line v_r = 0 from bottom to top in panel b, at the inner turnaround radius, r = r_in. Finally, when v_ej ≲ 720 km s⁻¹, the stars may also undergo subsequent turnarounds.

For the HVSs with v_ej ≲ 800 km s⁻¹, panel a shows that, after an initial phase where the star velocity is almost purely radial, with |v_ϑ| less than or around a few tens of km s⁻¹ regardless of v_r, |v_ϑ| quickly increases after the outer turnaround and becomes very large for those stars that undergo the inner turnaround. These large values of |v_ϑ| are determined by the exchange of kinetic energy between radial and angular degrees of freedom, rather than by the shape of the gravitational potential well, especially close to the inner turnaround. Therefore, the tangential velocity of these stars could in principle be used as an indicator of the shape of the DM halo only in specific parts of the star’s trajectory.

To perform our investigation on the shape of the DM halo by means of the HVS tangential velocities, we needed to select only the stars that, at t = t_obs, were on their first outward trajectory from the Galactic Center, namely the stars that had not experienced an inner turnaround yet. Therefore, we first excluded the stars that, at t = t_obs, were on an inward trajectory (v_r < 0) toward the Galactic Center: we thus required v_r > 0 for the stars of our sample. From all the outgoing stars, we then excluded those that had already experienced the inner turnaround and may thus had uninterestingly large |v_ϑ|. To do so, we note that panel b of Fig. 5 shows that stars characterized by ejection velocities v_ej ≲ 740 km s⁻¹ may live long enough to experience at least one inner turnaround. However, panel b of Fig. 5 also shows that these stars can never go back to galactocentric distances r > 10 kpc, after having undergone the inner turnaround. We thus required a galactocentric distance r > 10 kpc for the stars of our sample. This detailed analysis of the star kinematics thus supports our choice of the two preliminary selection criteria, v_r > 0 and r > 10 kpc, adopted in Sect. 3.1.

The case of a nonspherical DM halo is slightly more complicated than the case of the spherical halo explored so far. Here, both the disk and the DM halo generate a non-null latitudinal velocity, v_ϑ; furthermore, fully triaxial or spheroidal DM halos with a symmetry axis misaligned with respect to the z axis are also sources of a non-null azimuthal velocity, v_φ, because the axial symmetry of the Galactic potential is broken.

However, we again find that both v_ϑ and v_φ are strongly coupled with v_r for the stars with v_ej ≲ 800 km s⁻¹, as it happened for v_ϑ in the case of a spherical DM halo. As an example, for the case of a (q_y, q_z) = (0.8, 1) DM halo, Fig. 6 shows the relation between the azimuthal and radial velocity components. The stars’ travel time is again 160 Myr, namely the largest possible travel time for 4 M_⊙ stars. The figure shows that for stars with v_ej ≲ 800 km s⁻¹, |v_φ| starts increasing after the outer turnaround, and reaches very large values close to the inner turnaround, irrespective of the shape of the DM halo that generated the non-null v_φ’s. To exclude all the stars that had acquired large v_φ and correspondingly large v_θ after the outer turnaround, we again selected stars with v_r > 0 and r > 10 kpc, as we did in the case of a spherical DM halo.

Fig. 6. HVS galactocentric azimuthal velocity, vφ, as a function of the radial velocity, vr, for HVSs ejected with different speeds, vej, and traveling for 160 Myr. Here, the gravitational potential of the MW is non-axisymmetric, with a qy = 0.8 and qz = 1.0 DM halo potential. Initially, the star velocity is purely radial, and vφ = 0; as time goes on, vφ becomes non-null due to the DM halo asymmetry with respect to the z axis (see Eq. (5)). Radial and angular dynamics are coupled for stars with lower ejection speeds (vej ≲ 800 km s−1); these stars start acquiring significant vφ values after the outer turnaround. However, stars with larger ejection speeds die before this increase in vφ can happen. The dashed gray line corresponds to vr = 0.

We note that these two selection criteria, illustrated for both a spherical DM halo (see Fig. 5) and a nonspherical DM halo with (q_y, q_z) = (0.8, 1) (see Fig. 6), actually hold for a DM halo with any combination of axis ratios among those investigated in this work, as illustrated in Fig. 7. This figure shows the radial velocity, v_r, as a function of the galactocentric distance, r, for a star ejected at v_ej = 740 km s⁻¹ and traveling in a Galaxy with DM halos of different shapes: a spherical DM halo, an extremely oblate or prolate DM halo that is axisymmetric about the z axis, and an extremely oblate or prolate DM halo that is axisymmetric about the y axis. Regardless of the axis ratios, stars ejected at v_ej = 740 km s⁻¹ can never reach a galactocentric distance r > 10 kpc, after having experienced the inner turnaround.

Fig. 7. Galactocentric radial velocity, vr, as a function of the galactocentric distance, r, for HVSs ejected radially outward with an ejection velocity of vej = 740 km s−1 in the direction $\widehat{n}(\vartheta ,\phi )=\widehat{n}(45°,45°)$, for different shapes of the DM halo. The triaxiality parameters (qy, qz) of the DM halos are listed in the inset. The travel time is 160 Myr. The vertical and horizontal dashed lines correspond to r = 10 kpc and vr = 0, respectively.

We emphasize that the stars’ outer and inner turnaround radii depend on the gravitational potential chosen for the computation of the trajectories (see also Chakrabarty et al. 2022). Thus, the kinematic selection criteria derived in this work would differ in different models of the Galactic gravitational potential. In addition, the kinematic selection criteria depend on the mass of the HVSs. In particular, HVSs with M < 4 M_⊙ would have longer lifetimes, and thus a higher probability of reaching the outer turnaround, undergoing the inner turnaround, and going back to large galactocentric distances before dying out; this would result in a threshold radius larger than 10 kpc.

3.3. Effect of the initial conditions on the distributions of the shape indicators

In Sect. 3.1, we identified the distribution of the magnitude of the HVS latitudinal velocity, |v_ϑ|, as the only distribution of HVS phase-space coordinates to be significantly affected by a change of shape of the DM halo from spherical to spheroidal with axis of symmetry about the z axis. Based on this finding, and extending the argument to a more general, non-axisymmetric Galactic potential, we identified the components of v_t (i.e., v_ϑ and v_φ) as indicators of the triaxiality parameters of a DM halo. Hereafter, we generally refer to ω as each of the quantities used as an indicator of the shape of the DM halo; we define the distribution of the shape indicator ω as D_ω.

In Sect. 3.1 we presented results obtained for simulated HVSs from mock catalogs A, which are unambiguously determined by the shape of the DM halo because the set of stars is ejected with the same initial conditions in halos of different shapes. If a set of real HVSs were ejected with the initial conditions used to generate our mock catalogs A, recovering the shape of the MW DM halo from the distributions D_ω of real stars would be straightforward. In reality, however, the situation is more complex. At any given time t, we can observe the phase-space distribution of a sample of HVSs that are traveling in a DM halo whose unknown shape we want to recover. Furthermore, the ejection conditions of these stars are also unknown: even assuming an ejection mechanism (e.g., the Hills mechanism, in our case), we do not know the initial conditions of the stars’ trajectories, which are subject to statistical fluctuations. Therefore, recovering the shape of the DM halo requires the comparison of the distributions D_ω of the real HVS sample with mock D_ω’s that were generated with all possible combinations of initial conditions and shape of the DM halo. In the following, we illustrate the impact of the statistical fluctuations of the initial conditions of the ejected stars on the distributions D_ω and on our ability to detect deviations from spherical of the shape of the DM halo.

When simulating the trajectories of a sample of HVSs, a different set of initial conditions for the trajectories yields a different distribution of the HVS phase-space coordinates at t = t_obs. Consequently, the results of the KS test which compares the distribution D_ω of the shape indicators against a reference set of D_ω’s of DM halos with different degrees of triaxiality, depend on the initial conditions. In Sect. 3.1, we already explored the effect of using different sets of initial conditions, each of them applied to all mock catalogs A generated for a series of spheroidal DM halos with different shapes. These different sets are responsible for the fluctuations of the p-value of the KS test within the ranges reported in Table 2 for some of the HVS phase-space coordinates. As a result, mild deviations from spherical of the shape of the DM halo (i.e., 0.05 < |q_z − 1|< 0.1) may not be recognized.

The situation becomes more complex when different sets of initial conditions are applied to each mock catalog. To investigate this case, we resorted to our HVS mock catalogs B (see Sect. 2.4): these catalogs differ from one another both for the shape of the DM halo, as catalogs A, and for the set of the stars’ initial conditions.

Depending on the combination of initial conditions and triaxiality parameters that characterizes each mock catalog B, the comparison of the resulting D_ω’s for any of these catalogs with a catalog of a spherical halo, based on the KS test, may lead to two different incorrect conclusions: (i) p < α (with α = 5%), suggesting a deviation from the spherical halo, at the significance level α, even though both star samples have traveled in the same, spherical DM halo; (ii) p > α, suggesting that both the DM halos crossed by the two star samples are spherical, even though one of the two halos is not.

Situation (i) would never occur in the comparison of two catalogs A, with the same initial conditions: the KS test for two D_ω’s from different catalogs A would always yield p = 100% for identical halo shapes. Conversely, for two mock catalogs B, situation (i) can occur as a consequence of a fundamental property of the p-value: when we compare two D_ω’s with the KS test, if the null hypothesis is true (i.e., if the two D_ω’s are drawn from the same parent distribution), the value of the KS test statistics will be at least as large as the observed value in a fraction p of the cases (e.g., Press et al. 2007). Because the p-value of the KS test is a random variable itself, and because for a true null hypothesis it is uniformly distributed in the range [0; 1] (e.g., Hung et al. 1997; Donahue 1999; Bhattacharya & Habtzghi 2002), one time out of 20 the p-value will be ≤5% for D_ω’s that are drawn from the same parent distribution. At the significance level α = 5%, this probability implies that a null hypothesis that is actually true (i.e., the shapes of the two DM halos are equal) will be rejected in 5% of the cases³.

Situation (ii) can instead occur also for catalogs A, when the shapes of the DM halos are only mildly different. However, our success in detecting actual deviations from the spherical shape with the KS test applied to the D_ω’s is lower for catalogs B because of the additional effect of the statistical fluctuations of the initial conditions. Indeed, whereas in catalogs A we cannot distinguish a spherical DM halo from a spheroidal DM halo, symmetric about the z axis, whose deviation from spherical is |q_z − 1|< 0.1, with catalogs B the range of undetectable deviations widens to |q_z − 1|≲0.2. The failure to detect the different shapes of DM halos is a failure to reject a null hypothesis that is actually false⁴.

We note that false negatives may significantly hamper the recovery of the halo shape, especially for mild deviation from spherical of the shape of a spheroidal DM halo axisymmetric about the z axis. As an example, we mention that comparing the D_ω’s obtained in a spheroidal DM halo axisymmetric about the z axis and with q_z = 0.9 (i.e., |q_z − 1| = 0.1) against those obtained in a spherical DM halo yields a rate of false negatives of ∼20%: in other words, we are not able to distinguish the shapes of the two DM halos in ∼20% of the cases. Even though the rate of false negatives decreases with an increasing deviation from the spherical shape, and it becomes null for |q_z − 1|≳0.2, it is important to reduce the rate of false negatives to improve the efficiency of the shape recovery of the DM halo.

Summarizing, when the initial conditions of the trajectories of a set of HVSs are unknown, as it happens for real HVSs, comparing the D_ω’s of real stars with mock D_ω’s with a KS test may yield results on the shape of the DM halo that are significantly affected by the initial conditions of real HVSs, especially for mild deviations of the shape of the DM halo from spherical. Therefore, it is important to design a method that optimizes the recovery of the shape of the DM halo, minimizing the effects of the statistical fluctuations of the stars’ initial conditions. We propose this optimized method in the following.

4. A new method for constraining the shape of the DM halo

Our final goal is to recover the shape of the DM halo from the distribution D_ω of the shape indicators ω of a real sample of HVSs, properly accounting for the effects of the unknown initial conditions of the stars’ trajectories. To this aim, we designed a method based on (i) the use of the KS test to compare the D_ω of the sample of real HVSs with the corresponding distributions generated with a series of DM halos of different shape; (ii) the property of the KS test’s p-value mentioned in Sect. 3.3: its uniform probability density function for a true null hypothesis which, in our case, occurs when two D_ω’s are drawn from the same parent distribution.

To implement the method and evaluate its efficiency, we resorted to our HVS mock catalogs B (see Sect. 2.4), which differ from one another for the shape of the DM halo and/or the set of the stars’ initial conditions. We also constructed a sample of synthetic HVSs, hereafter referred to as the “observed sample”: the stars of this sample are ejected from the Galactic Center, according to the Hills mechanism, with a statistically random set of initial conditions, and move across a Galaxy whose DM halo has a known shape. In our analysis, the observed sample mimics a real sample of HVSs: we used it to test the efficiency of our method in recovering the correct, known shape of the DM halo from its D_ω at t = t_obs. The distributions of the kinematic properties of one of our observed samples is highlighted in green color in Fig. A.1. We emphasize that our HVS observed samples are ideal: we include neither observational uncertainties nor observational cuts imposed by the star observability, like its magnitude or position within the Galaxy. The method success rates that we estimate in Sects. 5.2 and 6.3 are thus valid for these samples alone and not for the HVS samples that might actually be observed.

We illustrate the basic concepts of our method in Sect. 4.1, and the method implementation in Sect. 4.2.

4.1. Fundamentals of the method

In our approach, recovering the shape of the DM halo crossed by an observed sample of HVSs requires the comparison of the D_ω of the observed sample with a series of corresponding D_ω’s generated in mock catalogs characterized by a different shape of the DM halo and by a different set of initial conditions. We considered the shape of the DM halo of the mock sample that best matches the observed sample as the actual shape of the DM halo crossed by the HVS observed sample.

In principle, the comparison could be performed by means of a KS test, whose null hypothesis H₀ states that the two compared D_ω’s are drawn from the same parent distribution. At the significance level α, we would accept H₀ in those cases where the test returns p > α, and we would reject it otherwise. Accepting H₀ would correspond to considering the D_ω of the observed sample indistinguishable from the mock D_ω selected for the comparison, thus associating to the observed sample a DM halo with the same shape of the mock sample’s DM halo. However, because of the effect of the statistical fluctuations of the initial conditions of the stars’ trajectories (see Sect. 3.3), the D_ω of the observed sample may either turn out to be indistinguishable from a significant number of mock D_ω’s obtained in DM halos with different shapes (false negatives) or turn out to be significantly different from the mock D_ω’s obtained in DM halos with identical shapes (false positives).

To select the “best match” between observed and mock sample, we resorted to a property of the p-value, already mentioned in Sect. 3.3. When the null hypothesis H₀ is true, the p-value is uniformly distributed in the range [0; 1] (e.g., Hung et al. 1997; Donahue 1999); thus, its median value is p_med = 0.5. On the other hand, when the alternative hypothesis, H_a, is true, the distribution of the p-values is markedly skewed toward low p-values, because small values of p are more likely; thus, the median p-value under H_a will be p_med ≪ 0.5.

As a consequence, performing the KS test for a number n_t of D_ω pairs that are randomly drawn from the same parent distribution, namely from the ensemble of the D_ω obtained from mock catalogs B that are characterized by the same shape of the DM halo but by different initial conditions, yields a uniform distribution of p-values. Conversely, performing the KS test for a number n_t of D_ω pairs that are randomly drawn from different parent distributions, namely from different mock catalogs B, each characterized by a different shape of the DM halo and different initial conditions, yields a distribution of p-values that is markedly skewed toward small p-values. The larger is the difference in shape, the more skewed is the distribution.

The situation becomes more complex when we pick a specific distribution of ω, say ${\stackrel{\sim }{D}}_{\omega }$, as that of the observed sample, and we compare it against a series of mock distributions, performing a number n_t of KS tests. Even though both the HVS observed sample and the mock samples have traveled in the same DM halo, the distributions of the p-values is not necessarily uniform. It may be approximately uniform over the range [0; 1], skewed toward low p-values, or skewed toward high p-values, with a corresponding median p-value p_med ≃ 0.5, p_med < 0.5, and p_med > 0.5, respectively. This situation is illustrated in Fig. 8, which shows the distributions of the p-values obtained from n_t = 5000 KS test comparisons of ${\stackrel{\sim }{D}}_{\omega }$ for each of n = 3 different observed samples against all the D_ω’s of a mock catalog B generated with a spherical DM halo. In each of the three cases, the observed sample and all the mock samples of HVSs have traveled in the same, spherical DM halo; however, the p-value distributions are markedly different: the black histogram shows the case of an approximately uniform distribution, with p_med ≃ 0.5; the olive histogram shows the case of a distribution skewed toward low p-values, with p_med ≃ 0.3; the purple histogram shows the case of a distribution skewed toward high p-values, with p_med ≃ 0.7.

Fig. 8. Distributions of the p-values obtained by comparing ${\stackrel{\sim }{D}}_{\omega }$ from each of n = 3 observed samples of HVSs against the Dω’s of a series of nt = 5000 mock samples. The observed samples and the mock samples have all crossed a spherical DM halo. The black histogram shows an approximately uniform distribution of p, with a median value pmed ≃ 0.5; the olive histogram shows the case of a skewed distribution with pmed ≃ 0.3; the purple histogram shows the case of a skewed distribution with pmed ≃ 0.7. The median p-values of the three distributions are marked by the vertical dashed lines.

The fact that a uniform distribution of p-values is not guaranteed even for the same shapes of the DM halo is a consequence of the random selection of the observed sample’s ${\stackrel{\sim }{D}}_{\omega }$ that we compare against all of the mock D_ω’s. If we repeat the exercise of Fig. 8 with a series of n > 3 observed samples, we find n different distributions of p, each of them with a different median p-value. However, for sufficiently large n, the sum of all these distributions yields a uniform distribution. Moreover, the distribution of the n median p-values is markedly skewed toward high p_med’s, as the one shown in Fig. 9: here, the gray histogram is obtained assuming in turn as observed sample each of the n_t = 5000 mock catalogs B with a spherical DM halo and comparing its ${\stackrel{\sim }{D}}_{\omega }$ against all the D_ω’s of the remaining mock catalogs, for a total of n = n_t = 5000 KS tests per observed sample, and n_t corresponding p-values whose distribution has a median p_med. Repeating the exercise n = n_t = 5000 times yields a distribution of n = n_t = 5000 median p-values. In this distribution, ∼80% of the p_med’s are in the range 0.3−0.7, and the median is M₀ ≃ 0.55.

Fig. 9. Distribution of the median p-values obtained from the KS test comparison of ${\stackrel{\sim }{D}}_{\omega }$ from all possible observed samples against the Dω’s from mock catalogs B generated in the same DM halo; the DM halo is spherical (gray histogram) or spheroidal with axis of symmetry about the z axis with an oblate (qz = 0.7; green histogram) or prolate (qz = 1.4; blue histogram) shape. The solid gray, green, and blue lines are the medians of the corresponding distributions.

We note that the characteristics of the distribution of the median p-values primarily depend on the validity of H₀; indeed, the distribution is independent of the shape of the DM halo under test. As shown in Fig. 9, the p_med’s have indistinguishable distributions for extremely prolate (q_z = 1.4; blue histogram), spherical (q_z = 1; gray histogram), and extremely oblate (q_z = 0.7; green histogram) DM halos. The compatibility of the three distributions is also confirmed by a KS test. For all of the distributions, the median p_med is M₀ ≃ 0.55.

Summarizing, Figs. 8 and 9 show how the KS test comparison of the ${\stackrel{\sim }{D}}_{\omega }$ of an observed sample of HVSs with a series of n_t mock D_ω’s may yield a p-value distribution different from uniform even though the observed sample and the mock sample crossed DM halos with the same shape. This variety of p-value distributions corresponds to a variety of p_med’s. Higher p_med’s are more likely than lower p_med’s, implying that a distribution of p-values skewed toward higher p’s is more likely to occur under the null hypothesis H₀. Conversely, a distribution of p-values skewed toward lower p’s, and thus characterized by a lower p_med, is a more likely outcome under the alternative hypothesis H_a that the DM halos have different shapes.

Therefore, to select the best match between the observed sample and the mock sample, we chose the p-value distribution characterized by the largest p_med. The DM halo crossed by the HVS observed sample was then assigned the shape of the mock sample that yielded the best match.

4.2. Implementation of the method

To implement our method, we first explored the case of an axisymmetric Galactic potential, which includes the case of a spherical DM halo and that of a spheroidal DM halo axisymmetric about the z axis; we then explored the more complex case of a non-axisymmetric Galactic potential, which includes the case of a spheroidal DM halo with a symmetry axis misaligned with respect to the z axis, and that of a fully triaxial DM halo.

In both the axisymmetric and non-axisymmetric case studies, we resorted to mock catalogs B (Sect. 2.4), which we constructed as follows. We defined a series of n_s reference shapes for the DM halo, characterized by: (i) triaxiality parameters q_y = 1 and q_z that varied in steps of 0.1 within the range 0.7−1.4, for the axisymmetric Galactic potential; (ii) triaxiality parameters q_y and q_z that both varied in steps of 0.1 within the range 0.7−1.4, and q_y ≠ 1, for the non-axisymmetric Galactic potential. For each of the n_s shapes of the DM halo, we generated an ensemble of n_t = 5000 HVS mock catalogs B, each of them including the phase-space coordinates of a sample of HVSs characterized by a different random set of initial conditions. We also generated one observed sample of HVSs, randomly drawn from one of the mock catalogs B: this sample plays the role of a real HVS sample, and was used to evaluate the efficiency of the method in recovering the known shape of the DM halo crossed by the observed sample.

For each of the n_s simulated shapes of the DM halo, we performed a number n_t = 5000 of KS test comparisons of the ${\stackrel{\sim }{D}}_{\omega }$ of the observed sample against the D_ω’s of the mock samples. From these n_t KS tests, we got a distribution of n_tp-values, whose median is p_med. We repeated the procedure for all the n_s simulated shapes, and eventually obtained a set of n_sp_med’s. We selected the mock sample corresponding to the largest value of p_med as the sample that best matched the observed sample, and we associated the shape of its DM halo with the DM halo of the observed sample.

To evaluate the success rate of our method in recovering the shape of the DM halo crossed by the observed sample, we generated a series of n = n_t observed samples corresponding to the same shape of the DM halo, by randomly varying the set of the HVS initial conditions, and we computed the fraction of cases where the method correctly recovers the known shape of the DM halo crossed by the observed sample. Finally, to study the dependence of the success rate of our method on the shape of the DM halo, we repeated our analysis for the observed samples of HVSs drawn from DM halos of different shapes.

In Sect. 5 we illustrate the application of our method to a series of observed samples generated in an axisymmetric Galactic potential, where the shape indicator ω is |v_ϑ|. The test used for our analysis is the two-sample, one-dimensional KS test that we also used in the previous sections. In Sect. 6 we illustrate the application of our method to a series of observed samples generated in a non-axisymmetric Galactic potential, where the shape indicators ω are |v_ϑ| and a function of v_φ. The test used for our analysis is the two-sample, two-dimensional KS test.

A possible alternative to the KS test chosen for our method is represented by the Anderson-Darling (AD) test (Anderson & Darling 1952, 1954). However, the AD test currently exists only in its one-dimensional formulation; the use of this test would thus be limited to the case of axisymmetric Galactic potentials, when the shape of the DM halo can be recovered from the distribution of one shape indicator only, |v_ϑ|. In the comparison of two one-dimensional distributions of |v_ϑ|, the AD test is more sensitive to the tails of the distributions than the KS test. These tails are sensibly affected by the shape of the DM halo, as can be seen from Fig. 4. Therefore, the use of the AD test for our method is expected to improve our method success rate (see Sect. 5.2) in recovering the correct shape of an axisymmetric DM halo. Consequently, our choice of the KS test as the statistical test of our method represents a conservative choice in terms of method success rate.

Unlike the AD test, the KS test is available also in its two-dimensional formulation. We can thus also apply the method to non-axisymmetric Galactic potentials, where the comparison of the distributions of the two shape indicators, |v_ϑ| and a function of v_φ, is required. Adopting the same statistical test guarantees a consistent comparison of the success rate in axisymmetric and non-axisymmetric scenarios.

We remind that, from mock catalogs B and for the observed sample, we always selected subsamples of stars that fulfill the criteria defined in Sect. 3.2, namely subsamples composed of N ≃ N_obs/10 ≃ 800 stars located at r > 10 kpc and with v_r > 0.

5. Constraining the shape of the DM halo in an axisymmetric Galactic potential

If the DM halo of our Galaxy is either spherical (i.e., with triaxiality parameters q_y = q_z = 1) or spheroidal with axial symmetry about the z axis (i.e., with triaxiality parameters q_z ≠ q_y = 1), the total gravitational potential of the Galaxy (Eq. (1)) is axisymmetric about the z axis. For an axisymmetric Galactic potential, we show that the method presented in Sect. 4 can effectively recover the axial ratio q_z of the DM halo from the distribution ${\stackrel{\sim }{D}}_{\omega }$ of the shape indicator ω of an observed sample of HVSs (Sect. 5.1). We also present the evaluation of the success rate of the method (Sect. 5.2).

As demonstrated in Sect. 3.1, if the Galactic potential is axisymmetric, there is only one shape indicator ω of the DM halo: the magnitude of the latitudinal velocity of the HVSs, |v_ϑ|. We refer to its distribution as D_|vϑ|. A few examples of D_|vϑ| for DM halos with different shapes were shown in Fig. 4.

5.1. Shape recovery

We defined a series of n_s = 8 reference shapes for the DM halo by varying q_z in steps of 0.1 in the range 0.7−1.4. For each shape, we generated an ensemble of n_t = 5000 mock catalogs B, one per different set of initial conditions of the star trajectories. We thus got n_t mock samples of stars, and n_t corresponding distributions D_|vϑ| for each of the n_s shapes of the DM halo.

As a first test, we chose as the HVS observed sample one random mock sample of HVSs that have crossed a spherical DM halo (with q_z = q_y = 1). We now show that our method successfully recovers the axis ratio q_z = 1.

For each of the n_s = 8 reference shapes, we performed a set of n_t = 5000 KS test comparisons of the observed sample’s ${\stackrel{\sim }{D}}_{|v_{\vartheta }|}$ against each of the n_tD_|vϑ|’s of the mock samples associated with that reference shape. Figure 10 shows the outcome of this test. Specifically, panel b of Fig. 10 shows the test result for three of these n_s sets of comparisons: the n_t KS test comparisons of the ${\stackrel{\sim }{D}}_{|v_{\vartheta }|}$ of the observed sample against the D_|vϑ|’s of the mock star samples that traveled in the spherical DM halo yields a uniform distribution of p-values (filled, gray histogram) whose median is p_med ≃ 0.5; the comparison of ${\stackrel{\sim }{D}}_{|v_{\vartheta }|}$ against a mock sample generated either in a slightly oblate (q_z = 0.9; yellow histogram) or in a slightly prolate (q_z = 1.1; magenta histogram) DM halo yields a distribution of p-values that is markedly skewed toward very low p’s; the median p-values are p_med = 2 × 10⁻² for the observed versus oblate comparison, and p_med = 10⁻⁴ for the observed versus prolate comparison.

Fig. 10. Distributions of the p-values for an HVS observed sample that crossed a spherical DM halo. Panel a: distributions of the p-values (in logarithmic scale) obtained from the KS test comparison of the distribution ${\stackrel{\sim }{D}}_{|v_{\vartheta }|}$ of one HVS observed sample generated in a spherical DM halo, against each of the ns = 8 ensembles of nt = 5000 mock distributions, D|vϑ|, generated in DM halos with different shapes, as listed in the panel. The vertical dashed lines mark the median p-value of each distribution. Panel b: enlargement of the rightmost part of panel a with the p-value axis in linear scale. The difference in the shape of the distributions in panels a and b is due to the different size of the histogram bins in the logarithmic and linear scales.

As shown in panel a of Fig. 10, the median p-value becomes smaller and smaller for increasing departure from spherical of the shape of the DM halo of the mock catalogs used for the comparison: it reaches p_med ≃ 10⁻¹³ for the comparison against the most oblate DM halo case (q_z = 0.7; green histogram) considered in this work, and p_med ≃ 10⁻²¹ for the comparison against the most prolate DM halo case (q_z = 1.4; blue histogram). We thus confirm that the largest p_med, which we obtained in the observed versus spherical comparison, identifies the correct shape of the DM halo crossed by the observed sample (i.e., the spherical shape, with q_z = q_y = 1).

5.2. Success rate S of the method

Even though the result shown in the above test is a likely result, it is not guaranteed. Indeed, as pointed out in Sect. 4, the p-value distributions of Fig. 10 are not unique, and depend on the specific observed sample that we pick for the comparison (i.e., from the set of random initial conditions of the simulated HVSs). Therefore, it may happen that only in a fraction of cases the largest p_med is the correct indicator of the shape of the DM halo. We defined this fraction as the success rate S of our method in recovering the correct axis ratio of the DM halo. We now illustrate the evaluation of the method success rate, S, for the case of a spherical DM halo (Sect. 5.2.1), and the investigation of the dependence of S on the shape of the DM halo (Sect. 5.2.2).

5.2.1. The case of a spherical DM halo

To evaluate the success rate S of our method in recovering the correct axis ratio q_z = 1, we constructed the distributions of all possible p_med’s that could be obtained by comparing the observed sample against the mock samples corresponding to the reference shapes of the DM halo. To generate these distributions, we simulated a series of n = n_t HVS observed samples in a spherical DM halo, by randomly varying the set of the stars’ initial conditions. For each observed sample, we performed the n_t KS test comparisons against all the mock samples of a given shape, and we obtained a p_med; performing the procedure for n = n_t observed samples yields a distribution of n = n_t values of p_med for each comparison of the observed samples against the mock samples of a given shape of DM halo. Repeating this exercise for each of the n_s = 8 shapes of DM halo yields the n_s distributions of p_med’s that we show in Fig. 11.

Fig. 11. Distributions of the median p-values for an HVS observed sample that crossed a spherical DM halo. Panel a: distributions of the median p-values, pmed (in logarithmic scale), obtained from the KS test comparison of the HVS observed sample generated in a spherical halo against each of the ns = 8 ensembles of mock HVS samples generated in a DM halo with different shapes. Each distribution is the result of the KS test comparison of the ${\stackrel{\sim }{D}}_{|v_{\vartheta }|}$’s of the n = nt observed samples against the nt = 5000 mock samples generated in a DM halo with different shapes, as listed in the panel. Panel b: enlargement of the rightmost portion of panel a, with the pmed axis in linear scale. The yellow and pink distributions are the only distributions with non-null overlap with the gray distribution. The different shape of the distributions in panels a and b is due to the different size of the histogram bins in the logarithmic and linear scales.

Panel a of Fig. 11 shows that the comparison of the observed samples (generated in a spherical DM halo) against the samples generated in a spherical DM halo returns a distribution (the gray histogram) where most of the p_med’s are larger than the p_med’s of the other distributions, confirming that the method correctly recovers the shape of the DM halo crossed by the observed sample in most of the cases. However, panel a also shows that the distribution corresponding to the comparison of the observed samples against the samples generated in the spherical DM halo displays a non-null overlap with the two distributions of p_med that correspond to the comparison of (i) the observed samples against the samples generated in a slightly oblate DM halo (q_z = 0.9; yellow histogram), and (ii) the observed samples against the samples generated in a slightly prolate DM halo (q_z = 1.1; magenta histogram). This non-null overlap, which can be better appreciated in panel b, implies that, for an observed sample generated in a spherical DM halo, the shapes that mildly deviate from spherical (|Δq_z| = 0.1, in our mock catalogs) might be erroneously associated with the observed sample of HVSs, based on our method. The moderate, non-null overlap, however, does not necessarily imply that erroneous associations do occur with a rate proportional to the overlapping areas.

To illustrate this concept, we can consider two overlapping distributions of p_med, one that is right-skewed and the other that is left-skewed, with a non-null overlap (see, e.g., the gray and yellow distributions of Fig. 11). Each value of p_med of the distributions is contributed by a specific HVS observed sample, characterized by a specific random set of initial conditions and a specific ${\stackrel{\sim }{D}}_{|v_{\vartheta }|}$ (see Sect. 5.1 for details). Only when the p_med that contributes to the left-skewed distribution (p_med, ls) is higher than the p_med that contributes to the right-skewed distribution (p_med, rs) for the same random observed sample, the erroneous association does occur. Other similar pairs (p_med, ls, p_med, rs), which can be randomly drawn from the overlapping portion of the two distributions, may never occur in reality. Therefore, whereas a null overlap ensures a null rate of erroneous shape associations, the existence of a non-null overlap of two distributions is only an indication that some erroneous associations may occur, without quantifying their rate of occurrence. However, the larger is the overlap, the higher is the probability of erroneous associations.

To evaluate the rate of success of our method in recovering the shape of a spherical DM halo, for each distribution of p_med’s (as the yellow and magenta distributions in Fig. 11) that has a non-null overlap with the p_med distribution corresponding to the observed versus spherical comparison (i.e., the gray distribution in Fig. 11), we computed the rate of occurrence of pairs (p_med, ls, p_med, rs) where p_med, ls > p_med, rs, thus providing an erroneous shape recovery. We found that in 1.16% of the cases (i.e., in 58 cases out of 5000) the KS test comparison of the ${\stackrel{\sim }{D}}_{|v_{\vartheta }|}$ of the observed sample against the n_t = 5000 D_|vϑ|’s of a mildly oblate (q_z = 0.9) DM halo yields a p_med value larger than that obtained in the comparison against the n_t = 5000 D_|vϑ|’s of a spherical DM halo (i.e., p_med, ls > p_med, rs), leading to erroneously classify as mildly oblate a spherical halo; this result is equivalent to a success rate S = 98.84%. The fraction of erroneous associations drops to 0.4% for the comparison against a mildly prolate (q_z = 1.1) DM halo, yielding S = 99.6%. For |Δq_z|≥0.2, the overlap of the p_med distributions with the distribution corresponding to the spherical DM halo is null, implying a null rate of erroneous associations, and a success rate of our method S = 100%.

Overall, for an observed sample generated in a spherical DM halo, our method enables the correct axis ratio, q_z = 1, of the DM halo to be recovered in more than 98% of the cases; in other words, the method has a success rate S ≳ 98%. We stress that, in the rare cases of an erroneous shape association, the recovered axis ratio q_z is off by only |Δq_z| = 0.1.

5.2.2. Dependence of the success rate S on the shape of the DM halo of the observed sample

The success rate S of our method displays a weak, non obvious dependence on the actual shape of the DM halo crossed by the HVS observed sample. Indeed, repeating the exercise of Sect. 5.2.1 for an HVS observed sample that traveled in each of our n_s = 8 reference DM halos (with q_z = 0.7−1.4), returned a success rate that varies from 89% to 99%, as listed in Table 3. The success rate is slightly higher for spherical (q_z = 1.0) or mildly prolate (q_z = 1.1) DM halos than for oblate and markedly prolate DM halos; in other words, our method recovers the shape of a DM halo crossed by an observed sample of HVSs more easily when this DM halo is not too different from spherical. In addition, S is slightly lower for markedly prolate than for markedly oblate DM halos; this result indicates that the method recovers the shape of an oblate DM halo more easily than the shape of a prolate DM halo.

To illustrate this effect, Fig. 12 shows the analogs of the distributions of Fig. 11 for the case of an HVS observed sample that traveled in a prolate DM halo with q_z = 1.2. Panel a shows that the comparison of the HVS observed sample against the mock sample that traveled in a DM halo with the same shape returns a distribution (the pink histogram) where most of the p_med’s are larger than the p_med’s of the other distributions. This result confirms that the correct shape is recovered also for nonspherical DM halos.

Fig. 12. Distributions of the median p-values for an HVS observed sample that crossed a prolate DM halo. Panel a: distributions of the median p-values, pmed (in logarithmic scale), obtained from the KS test comparison of the HVS observed sample against each of the ns = 8 ensembles of mock HVS samples generated in a DM halo with different shapes. Each distribution is the result of the KS test comparison of the ${\stackrel{\sim }{D}}_{|{\mathit{v}}_{\vartheta }|}$’s of the n = nt observed samples obtained in a prolate DM halo with qz = 1.2 against the nt = 5000 mock samples generated in a DM halo with different shapes, as listed in the panel. Panel b: enlargement of the rightmost part of panel a with the pmed axis in linear scale. The cyan, blue, and magenta distributions are the only distributions with non-null overlap with the pink distribution. The difference in shape of the distributions in panels a and b is due to the different size of the histogram bins in the logarithmic and linear scales.

However, panel a and the corresponding enlargement in panel b also show that three distributions of p_med display a non-null overlap with the distribution corresponding to the correct shape: two of these distributions correspond to the comparison of the observed samples against the samples generated in more prolate DM halos (q_z = 1.3 and q_z = 1.4; cyan and blue histograms, respectively); the third one corresponds to the comparison against a less prolate DM halo (q_z = 1.1; magenta histogram). The number of overlapping distributions is thus larger (three instead of two) than in the case of the spherical DM halo illustrated in Fig. 11; furthermore, the overlapping area of the distributions is also larger than in Fig. 11 for the same Δq_z. Consequently, when the HVS observed sample traveled in a markedly prolate DM halo with q_z = 1.2, we get a larger number of erroneous shape associations for this halo than in the case of HVSs traveling in a spherical halo.

We note that, while the larger number of overlapping distributions is only a characteristics of the p_med distributions associated with markedly prolate (q_z ≥ 1.2) DM halos, the larger overlapping area is a property shared by all the oblate and markedly prolate DM halos. The excess in the number of overlapping distributions, where present, determines erroneous associations of the shape of the DM halo with spheroids whose q_z differs by 0.2 from the true q_z; however, it negligibly affects the success rate S of our method (≲0.02%). On the other hand, the larger overlapping area of the p_med distributions associated with DM halos whose q_z differs by 0.1 from the true q_z is the main responsible for the shape dependence of S, because the larger area increases the probability that a DM halo is assigned a shape whose q_z is off by 0.1. Summarizing, in the infrequent cases of erroneous shape associations, the recovered axis ratio q_z is typically off by |Δq_z| = 0.1, although |Δq_z| = 0.2 can seldom occur (≲0.04% of the cases) when the DM halo crossed by the HVS sample is markedly prolate.

As illustrated at the beginning of this section, the weak dependence of the success rate S on the shape of the DM halo crossed by the HVS observed sample manifests itself with a slightly higher S for spherical and mildly prolate DM halos, and a slightly lower S for markedly prolate than for markedly oblate DM halos. This weak shape dependence is a direct consequence of the following facts: (i) the difference between the D_|vϑ|’s of HVSs that traveled in a spherical or in a mildly prolate DM halo (q_z = 1.0−1.1) and the D_|vϑ|’s generated in DM halos with $q_z^{\prime }=q_z\pm 0.1$ is more pronounced than the difference between the D_|vϑ|’s generated in DM halos that are either oblate (q_z = 0.7−0.9) or markedly prolate (q_z = 1.2−1.4) and the D_|vϑ|’s produced in DM halos with $q_z^{\prime }=q_z\pm 0.1$; (ii) the differences among the D_|vϑ|’s of markedly prolate halos are milder than those among oblate halos. We ascribe these two effects to a combination of (a) the choice of a fixed resolution in q_z of our mock catalogs (i.e., Δq_z = 0.1), and (b) to the superposition of the gravitational actions of the disk and of the DM halo.

Specifically, effect (i) is mostly due to reason (a), which is the resolution in q_z of our mock catalogs. Indeed, more and more prolate (oblate) DM halos, obtained from the spherical halo (q_z = 1.0) by progressively increasing (decreasing) q_z by Δq_z = 0.1, generate a response in the D_|vϑ|’s that is stronger when q_z is closer to 1. This is true independently of the presence of the axisymmetric disk potential.

On the other hand, effect (ii) is due to a combination of reasons (a) and (b). Reason (a) causes the markedly prolate DM halos considered in our mock catalogs to have q_z’s (1.2, 1.3, and 1.4) that render their shapes more similar to one another than the q_z’s of the oblate DM halos considered in this work (0.7, 0.8, and 0.9); indeed, the DM halos whose axis ratios are the reciprocals of the axis ratios of these oblate halos would be characterized by q_z’s equal to ≃1.1, 1.25 and 1.4, respectively. Thus, the D_|vϑ|’s generated by the markedly prolate DM halos are expected to be less different from one another than the corresponding distributions obtained for the oblate DM halos.

On top of this effect, the gravitational pull of the disk makes the D_|vϑ|’s generated by the markedly prolate DM halos even less different from one another. Indeed, in a prolate DM halo the gravitational pull of the halo is opposed to the pull of the disk. For q_z ≤ 1.2, the presence of the DM halo increases the fraction of HVSs with low |v_ϑ|’s, rendering the distributions D_|vϑ|’s more left-skewed than that of a spherical DM halo. However, for q_z ≥ 1.3, the action of the DM halo not only compensates for the gravitational pull of the disk, it overcomes this pull, by bending the HVS trajectories toward the z axis: consequently the fraction of low |v_ϑ|’s drop – because v_ϑ increases in magnitudes and changes sign – and renders the D_|vϑ|’s less left-skewed and more similar to those of the less prolate DM halos.

The above effects are responsible for a slightly higher success rate S for oblate DM halos. While the effects of the choice of a fixed Δq_z to build the mock catalogs can easily be overcome, the effect of the combined gravitational actions of the disk and of the DM halo is inherent to the problem.

Overall, in an axisymmetric Galactic potential, our method recovers the axis ratio q_z of the DM halo crossed by an observed sample of HVSs with a success rate S ≳ 89%, and the erroneous shape association imply q_z that is off by ±0.1 in the overwhelming majority (≳99.96%) of the cases. In a negligible fraction (≲0.04%⁵) of the cases, and for markedly prolate DM halos only, q_z can be off by 0.2.

6. Constraining the shape of the DM halo in a non-axisymmetric Galactic potential

If our Galaxy has either a fully triaxial DM halo (i.e., triaxiality parameters q_y ≠ q_z, with q_y ≠ 1 and q_z ≠ 1) or a spheroidal DM halo with a symmetry axis misaligned with respect to the z axis (i.e., triaxiality parameters q_y ≠ q_z = 1, or q_y = q_z ≠ 1), the total gravitational potential of the Galaxy (Eq. (1)) is non-axisymmetric. For a non-axisymmetric Galactic potential, we show that the method presented in Sect. 4 can effectively recover the axial ratios q_y and q_z from the distribution ${\stackrel{\sim }{D}}_{\omega }$ of the shape indicators ω of an observed sample of HVSs (Sect. 6.2). We also present the evaluation of the success rate of the method (Sect. 6.3).

As anticipated in Sect. 3.1, in a non-axisymmetric Galactic potential both components of the tangential velocity v_t = (v_ϑ, v_φ) are affected by the halo triaxiality, and can thus be used as indicators of the shape of the DM halo. Specifically, we identify two shape indicators ω: the magnitude of the latitudinal velocity of the HVSs, |v_ϑ|, and a function ${\overline{v}}_{\phi }$ of the azimuthal velocity v_φ, which we define in Sect. 6.1. Hereafter, the distributions of the two shape indicators |v_ϑ| and ${\overline{v}}_{\phi }$ will be referred to as D_|vϑ| and $D_{{\overline{v}}_{\phi }}$, respectively. The corresponding two-dimensional distribution will be referred to as $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$.

6.1. |v_ϑ| and ${\overline{\mathit{v}}}_{\phi }$: Two indicators of the shape of the DM halo

In a Galaxy with a non-axisymmetric gravitational potential, D_|vϑ| and $D_{{\overline{v}}_{\phi }}$ are both affected by each of the two triaxiality parameters, q_y and q_z.

The behavior of D_|vϑ| in the case of a non-axisymmetric Galactic potential is similar to the behavior of D_|vϑ| in an axisymmetric Galactic potential (Sect. 3.1).

Figure 13 shows D_|vϑ| for a sample of HVSs that traveled in a Galaxy with DM halos of different shapes: the comparison of different distributions in each of the two panels shows that increasing q_z at fixed q_y leads to D_|vϑ|’s that are generally more skewed toward low values of |v_θ|. This effect was already shown in Fig. 4 for the case of spheroidal DM halos axisymmetric about the z axis (i.e., with q_y = 1): the gravitational pull of the disk, which drives the HVSs toward the x–y plane, is more and more compensated for by a DM distribution that is more and more elongated in the direction of the z axis. We note that the increase of skewness with increasing q_z depends on the value of q_y and stops when the action of the DM halo overcomes that of the disk, making v_ϑ change sign and increase in magnitude. This effect can be noticed in the left panel of Fig. 13, where the D_|vϑ|’s are comparable for q_z = 1.0 and q_z = 1.4.

Fig. 13. Distributions of the magnitude of the galactocentric tangential velocity |vϑ| for HVSs that have traveled in DM halos with qz = {1.4, 1.0, 0.7}, and with qy = 0.7 (left panel) and qy = 1.4 (right panel). The distributions were generated with the same initial conditions (mock catalogs A) to highlight the effect of the different geometries of the DM halo.

On the other hand, a comparison of the distributions corresponding to the same value of q_z (i.e., the distributions drawn with the same color) in the two panels in Fig. 13 shows the dependence of D_|vϑ| on q_y, at fixed q_z: DM halos with lower q_y at fixed q_z imply larger concentration of DM away from the x–y plane, and thus generate (i) D_|vϑ|’s that are more skewed toward lower |v_ϑ|’s, as long as q_z is not too high and the DM halo only compensates for the |v_ϑ|’s induced by the disk, (see, e.g., the green and gray histograms, corresponding to q_z = 0.7 and q_z = 1.0); (ii) D_|vϑ|’s that are less skewed toward lower |v_ϑ|’s, for large values of q_z, when the pull of the DM halo overcomes the pull of the disk, and the increase in |v_ϑ| previously discussed is enhanced by a DM distribution more concentrated about the z axis (see, e.g., the blue histograms, corresponding to q_z = 1.4).

A non-null distribution of v_φ is a distinctive characteristic of non-axisymmetric Galactic potentials: because in our model the gravitational potential is axially symmetric for the Galactic disk and spherically symmetric for the bulge, the only source of nonzero v_φ is the triaxial DM halo with q_y ≠ 1. This DM halo can be either a fully triaxial DM halo or a spheroidal DM halo with a symmetry axis misaligned with respect to the z axis (see Sect. 2.2)⁶.

Even though the distribution of v_φ is extremely sensitive to the triaxiality parameters of the DM halo, using the very value of v_φ as a shape indicator leads to a degeneracy problem: a triaxial ellipsoid with a given q_y = q_y, 1 and q_z = q_z, 1 is equivalent, in terms of geometric shape, to a triaxial ellipsoid characterized by q_y = q_y, 2 = 1/q_y, 1 and q_z = q_z, 2 = q_z, 1/q_y, 1; indeed, the role of the semimajor axes a and b is swapped in these two ellipsoids, and a rotation of 90° about the z axis would make one of the two ellipsoids coincide with the other. As a consequence, the resulting distributions of v_φ are statistically indistinguishable. This effect can be seen in the left panel of Fig. 14, where we show the case of two samples of HVSs that traveled in two triaxial gravitational potentials whose semimajor axes a and b are swapped: one potential has (q_y, 1, q_z, 1) = (0.7, 1.0) and the other has (q_y, 2, q_z, 2) = (1.4, 1.4)≃(1/q_y, 1, q_z, 1/q_y, 1).

Fig. 14. Distributions of the azimuthal velocity vφ of two samples of HVSs that have traveled in gravitational potentials whose DM halos have the same geometrical shape but differ by 90° in azimuthal orientation. Left panel: distributions of vφ indistinguishable from one another when the stars from all the quadrants are considered. Right panel: distributions of vφ manifestly different when we consider the HVSs located in the first quadrant of the x–y plane (0° < φ < 90°) only. The HVSs in the first quadrant have positive (negative) vφ when qy = 1.4 (0.7). The distributions were generated with the same initial conditions (mock catalogs A) to highlight the effect of the different geometries of the DM halo.

We note that two DM halos that are degenerate in v_φ are also degenerate in |v_ϑ|. This effect can be seen in Fig. 13, where the D_|vϑ|’s corresponding to a DM halo with (q_y, q_z) = (0.7, 1.0) (left panel, gray histogram) and to a DM halo with (q_y, q_z) = (1.4, 1.4) (right panel, blue histogram) are indistinguishable.

This degeneracy problem does not limit our understanding of the halo geometric shape in a strict sense; it rather hampers our ability of discriminating, for a DM halo with a given geometry, between two halo orientations that differ by 90° within the adopted reference frame. Breaking this degeneracy would thus enable not only the degree of triaxiality of the DM halo to be constrained, but also the orientation of the halo.

The degeneracy might be broken by considering only the HVSs that are traveling in one of the four quadrants of the x–y plane. Indeed, when q_y > 1, the mass distribution is elongated in the direction of the y axis; thus, the stars acquire a v_φ that drives them toward the y axis; conversely, when q_y < 1, the stars are attracted toward the x axis. Therefore, when the stars are located in the first quadrant (i.e., they have azimuthal coordinate 0° < φ < 90°) they have positive v_φ when q_y > 1 and negative v_φ when q_y < 1, in our sign convention. As shown in the right panel of Fig. 14, the distributions of v_φ for the HVSs located in the first quadrant are manifestly different and not overlapping: for q_y = 0.7 the v_φ’s are all negative, while they are all positive for q_y = 1.4.

Choosing only the HVSs located in one quadrant is however not the best solution to break the above mentioned degeneracy, because the size of the sample is reduced by a factor of ∼4, lowering the success rate of our method (see Sect. 7). To recover the original sample size, we defined as shape indicator the variable ${\overline{v}}_{\phi }\equiv v_{\phi }\frac{\tan \phi }{|\tan \phi |}$, where the factor $\frac{\tan \phi }{|\tan \phi |}$ is a sign plus or minus that renders the value of v_φ positive when the star is moving toward the y axis (q_y > 1), and negative when the star is moving toward the x axis (q_y < 1), independently of the quadrant where the star is located.

Figure 15 shows the distributions of ${\overline{v}}_{\phi }$ for two pairs of HVS samples that traveled through DM halos with the same geometric shape but azimuthal orientations that differ by 90°. One pair is composed of the HVS samples that traveled in DM halos with (q_y, q_z) = (0.7, 1.0) (green histogram) and (q_y, q_z) = (1.4, 1.4) (blue histogram) (i.e., the same samples investigated in Fig. 14); the other pair is composed of the HVS samples that traveled in DM halos with (q_y, q_z) = (0.9, 1.0) (yellow histogram) and (q_y, q_z) = (1.1, 1.1) (magenta histogram). The figure shows how the use of the variable ${\overline{v}}_{\phi }$ instead of v_φ enables us to easily distinguish distributions of azimuthal velocities generated in pairs of DM halos that have the semimajor axes a and b swapped. Using the two-dimensional distribution $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ thus enables us to overcome the degeneracy problem.

Fig. 15. Distributions of ${\overline{v}}_{\phi }\equiv v_{\phi }\frac{\tan \phi }{|\tan \phi |}$ for samples of HVSs that traveled in gravitational potentials whose DM halo has different triaxiality parameters. The use of the variable ${\overline{v}}_{\phi }$ instead of vφ enables us to easily distinguish distributions of azimuthal velocities generated in DM halos with the same geometry but whose semimajor axes a and b are swapped. The distributions were generated with the same initial conditions (mock catalogs A) to highlight the effect of the different geometries of the DM halo.

We note that the HVSs represented in Fig. 15 attain values of ${\overline{v}}_{\phi }$ that are of a few km s⁻¹ for mild q_y deviations (i.e., |Δq_y| = 0.1) from unity, but may reach ∼𝒪(10) km s⁻¹ for the most extreme values of q_y considered in our study (i.e., q_y = 0.7 and q_y = 1.4). This result is a consequence of the extreme flattening of the DM halo along either the x axis or the y axis for values of q_y that are extremely large or extremely small, respectively.

The distributions of |v_ϑ| and ${\overline{v}}_{\phi }$ are both very sensitive to the triaxiality parameters of the DM halo. Therefore, their combination can provide a powerful tool to detect non-sphericity of the DM halo and constrain the triaxiality parameters of the corresponding gravitational potential, as we show below.

We compared the two-dimensional distributions $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$’s of the shape indicators obtained for HVSs that traveled in DM halos of different shapes. We generated the distributions by using mock catalogs A, which were produced with HVS samples characterized by the same random set of ejection conditions, as we did in Sect. 3.1 for the case of spheroidal DM halos. We defined a series of n_s = 56 total reference shapes by varying both q_y and q_z in steps of 0.1 in the range [0.7; 1.4] and imposing q_y ≠ 1. These reference DM halos have either a fully triaxial shape or a spheroidal shape that is symmetric about the x axis or the y axis, and lead to a non-axisymmetric Galactic potential. For all these reference shapes, we generated the corresponding HVS mock catalogs.

We compared the $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$’s obtained from each of this n_s HVS samples against the $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ obtained for a spherical DM halo, where q_y = q_z = 1, and $D_{{\overline{v}}_{\phi }}$ is a distribution of null values. For this comparison, we used the two-sample, two-dimensional KS test (Press et al. 2007, and references therein), hereafter referred to as “2D KS test.” For all the 2D KS test comparisons, we found null p-values, implying that the $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$’s are significantly different from each other even for the smallest differences in triaxiality parameters with respect to the spherical DM halo. The null p-value results from the fact that, in a spherical DM halo, the HVSs do not acquire any azimuthal velocity; any non null $D_{{\overline{v}}_{\phi }}$ is thus a proof that the HVS sample has traveled in a non-axisymmetric gravitational potential. This result proves that the distribution of |v_ϑ| and ${\overline{v}}_{\phi }$ can effectively detect nonspherical shapes of the DM halo.

A 2D KS test comparison of pairs of $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ obtained for DM halos characterized by the same value of q_z but different values of q_y yielded p-values always smaller than 5% even for differences in q_y as small as Δq_y = 0.1. The high sensitivity of the $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ to small differences in q_y at fixed q_z is the result of the high sensitivity of $D_{{\overline{v}}_{\phi }}$ to q_y. In turn, this sensitivity comes from the fact that q_y ≠ 1 is the only source of non-null v_φ in triaxial DM halos.

Conversely, small differences in q_z at fixed q_y can lead to comparable $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$, according to a 2D KS test, with p-values that can exceed 5%. These larger values of p originates from the fact that differences in q_z affect v_ϑ, but v_ϑ is also affected by the gravitational pull of the disk. This effect was already discussed for the case of an axisymmetric Galactic potential, with a DM halo that is either spherical or spheroidal about the z axis, (Sect. 3.1). Overall, ${\overline{v}}_{\phi }$ is a more powerful shape indicator than |v_ϑ|, and the combination of the two shape indicators is the appropriate tool to constrain the triaxiality parameters of a DM halo with the HVSs.

We note that, in our model of the gravitational potential of the MW, we assumed the DM halo to have its principal axes aligned with the axes of our Cartesian reference frame, x, y, and z, where x indicates the direction from the Sun to the Galactic Center and z is orthogonal to the Galactic plane (see Sect. 2.2). Releasing this assumption has an effect on the shape indicator ${\overline{v}}_{\phi }$. Indeed, if we assume that one of the principal axes of the DM halo still coincides with the z axis, while the remaining two principal axes lie on the Galactic plane but are misaligned with respect to our x and y axes by an angle 0 < ϕ₀ ≤ 45°, the HVS ${\overline{v}}_{\phi }$’s would no longer be all negative or all positive. Specifically, for 0° < ϕ₀ ≲ 45°, higher ϕ₀ would correspond to higher degrees of mixing of negative and positive values, while for ϕ₀ ≃ 45°, the ${\overline{v}}_{\phi }$ distribution will be about half positive and half negative.

For very small misalignments, the very low degree of mixing of the ${\overline{v}}_{\phi }$ distribution would not prevent us from breaking the degeneracy in the orientation of the DM halo. On the other hand, for large misalignments, we would not be able to distinguish two DM halos with the same geometrical shape but with orientations that differs by 90°. However, we stress that even in those cases the degree of triaxiality of the DM halo could still be determined from the distribution of v_φ. Law et al. (2009) find that the axes of the DM halo on the plane of the Galactic disk are aligned with the x and y axes within 15°. In this case, our method would not encounter degeneracy problems.

6.2. Shape recovery

To recover the shape of the DM halo by means of the distribution ${\stackrel{\sim }{D}}_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ of the shape indicators |v_ϑ| and ${\overline{v}}_{\phi }$, we made use of the series of n_s = 56 reference shapes for a non-axisymmetric Galactic potential obtained, as described above, by varying both q_y and q_z in steps of 0.1 in the range 0.7−1.4, and imposing q_y ≠ 1. For each shape, we generated an ensemble of n_t = 5000 mock catalogs B, one per different set of initial conditions of the stars’ trajectories. From each mock catalog, we obtained one two-dimensional distribution $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ of the shape indicators, for a total of n_t distribution $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$’s per shape of the DM halo.

As an example, we chose as the HVS observed sample one random mock sample of HVSs that traveled in a DM halo with q_y = 1.2 and q_z = 0.9. We show here that our method successfully recovers the correct triaxiality parameters q_y and q_z.

Following a procedure similar to that described in Sect. 5.1, for each of the n_s = 56 reference shapes (q_y, q_z) of the DM halo we performed the 2D KS test comparisons of the observed sample’s ${\stackrel{\sim }{D}}_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ against all of the n_t = 5000 D_|vϑ|,vφ’s corresponding to each shape. Thus, for each shape of the DM halo, we obtained a distribution of n_tp-values and a corresponding p_med.

Table 4 shows the p_med’s obtained in this exercise for a selection of DM halos of different triaxiality parameters. We expect that the highest p_med points at the shape of the DM halo that best matches the shape of the DM halo actually crossed by the observed sample (see Sect. 5.1). Indeed, the highest p_med = 0.443 occurs for the DM halo with triaxiality parameters q_y = 1.2 and q_z = 0.9, namely for the DM halo actually crossed by the observed sample: this result demonstrates that our method effectively recovers the correct shape of the triaxial DM halo crossed by the observed sample.

6.3. Success rate

In this subsection we illustrate the success rate S of our method for the case of a DM halo with a specific shape (Sect. 6.3.1) and the dependence of S on the shape of the DM halo (Sect. 6.3.2).

6.3.1. The case of the triaxial DM halo with (q_y, q_z) = (1.2, 0.9)

To evaluate the success rate S of our method, we followed a procedure similar to that of Sect. 5.2. We generated a series of n = n_t HVS observed samples in a DM halo with q_y = 1.2 and q_z = 0.9, by randomly varying the set of the stars’ initial conditions. For each observed sample, we performed the n_t 2D KS test comparisons of the observed sample’s ${\stackrel{\sim }{D}}_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ against the $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$’s of all the n_t mock samples corresponding to a given reference shape of the DM halo, and we obtained a p_med. Performing the procedure for n = n_t observed samples yields a distribution of n = n_t values of p_med for each given shape of DM halo. Repeating this procedure for each of the n_s = 56 shapes of DM halo yields n_s distributions of p_med’s.

We show a representative selection of these distributions in Fig. 16. Panels a and b show that the highest p_med’s correspond to the distribution obtained for the DM halo with triaxiality parameters q_y = 1.2 and q_z = 0.9 (gray histogram), namely the halo traveled by the observed sample: this proves that our method recovers the correct shape of the triaxial halo in the large majority of the cases.

Fig. 16. Distributions of the median p-values for an HVS observed sample that crossed a triaxial DM halo. Panel a: distributions of the median p-values, pmed (in logarithmic scale), obtained from the 2D KS test comparison of the HVS observed sample against a selection of nine of the ns = 56 ensembles of mock HVS samples generated in DM halos with different triaxiality parameters. Each distribution is the result of the 2D KS test comparisons of the ${\stackrel{\sim }{D}}_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$’s of the n = nt = 5000 observed samples obtained in a DM halo with qy = 1.2 and qz = 0.9 against the nt mock samples generated in a DM halo with a different shape, as listed in the panel. Panel b: enlargement of the rightmost part of panel a, with the pmed axis in linear scale. The green and dark green distributions are the only distributions with non-null overlap with the gray distribution. The different shapes of the distributions in panels a and b are due to the different size of the histogram bins in the logarithmic and linear scales.

There are, however, cases where an erroneous shape association can occur. As shown in panel a of Fig. 16, the distributions of p_med’s, corresponding to the comparison of the ${\stackrel{\sim }{D}}_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$’s of the HVS observed sample against the $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$ of the HVS samples generated in DM halos with the same q_y = 1.2 and with q_z = 0.8 or q_z = 0.9, display a non-null overlap with the distribution corresponding to the comparison of the observed sample against the samples from DM halos with q_y = 1.2 and q_z = 0.9. This overlap can be better appreciated in the enlargement of panel b. The overlap implies that, for an observed sample generated in a DM halo with q_y = 1.2 and q_z = 0.9, our method can erroneously associate with the observed sample only the shapes characterized by the same q_y and by a q_z slightly different (|Δq_z|≤0.1) from the correct q_z.

Erroneous shape associations are however very rare, because the overlap between the distributions is very modest. In the above example, the success rate of our method is 99.98%; only in one case over 5000 the highest p_med suggests triaxiality parameters for the DM halo crossed by the observed sample that are not the correct ones, namely (q_y, q_z) = (1.2, 1.0) instead of (q_y, q_z) = (1.2, 0.9). In the analyzed case, any difference Δq_y ≥ 0.1 leads to distributions of p_med that do not overlap, namely the rate of erroneous shape associations is equal to zero; the axis ratio q_y is thus correctly recovered in 100% of cases.

6.3.2. Dependence of S on the shape of the DM halo of the observed sample

To investigate the effect of the shape of the DM halo of the observed sample on the success rate of our method in the case of a non-axisymmetric Galactic potential, we explored three additional shapes for the DM halo. Table 5 reports the axis ratios q_y and q_z of these DM halos, and the corresponding success rate, S. The axis ratio and success rate of the case explored in Sect. 6.3.1 are reported in the same table for completeness.

We obtained success rates S > 96% in all the explored cases. Thus, for a non-axisymmetric Galactic potential, the success rate of our method is less sensitive to the actual shape of the DM halo than for an axisymmetric Galactic potential. As in the case of an axisymmetric Galactic potential, illustrated in Sect. 5.2.2, for a non-axisymmetric Galactic potential we found that the DM halo with the largest q_z, namely the case (q_y, q_z) = (1.2, 1.3), yields the smallest success rate, although it is still larger than 96%.

Even though we investigated only four shapes of DM halos, we selected them to appropriately cover the axis-ratio space. We thus expect S ≳ 95% also for different combinations of axis ratios. In particular, less extreme axis ratios are expected to increase the success rate of the method: this tendency, already shown in Sect. 5.2.2 for the DM halo yielding an axisymmetric Galactic potential, is enhanced in the case of a DM halo that generates a non-axisymmetric Galactic potential, because of the presence of a non-null distribution of ${\overline{v}}_{\phi }$. Therefore, a two-dimensional distribution of shape indicators makes the constraining power of our method higher and more effective against the actual shape of the DM halo that we want to recover.

7. Sample size and method success rate

In Sects. 5 and 6 we investigated the efficiency of our method in recovering the shape of the DM halo, in both an axisymmetric and a non-axisymmetric Galactic potential, from a mock observed sample of HVSs composed of ∼800 4 M_⊙ stars. The size of the mock sample, dictated by the combination of the HVS ejection rate and the lifetime of the stars (see Sect. 2.3), as well as by our selection criteria (see Sect. 3.1), represents the optimal situation, which would occur if the HVSs were actually ejected according to the Hills mechanism with the assumed rate, and if we were able to observe all of them.

Here, we present our investigation of the dependence of the success rate S of our method on the size of the HVS observed sample. In Sect. 7.1, we show the results obtained when repeating the analysis performed in Sect. 5 on a series of mock HVS observed samples that traveled in a spherical DM halo and whose size is 50%, 25%, 10%, and 5% of the original sample size. In Sect. 7.2, we illustrate the generalization of our results to different shapes of the DM halo.

7.1. Spherical DM halos

We applied our method to an HVS observed sample that traveled in a spherical DM halo. For the case of an observed sample whose size is 50% the original size, panel a of Fig. 17 shows the distributions of the median p-values obtained by comparing each observed sample with the 5000 mock samples generated in spheroidal DM halos axisymmetric about the z axis. This figure is the analog of Fig. 11, with the exception of the reduced sample size of the observed sample and of the smaller number of p_med distributions.

Fig. 17. Distributions of the median p-values for an HVS observed sample of ≃400 stars that crossed a spherical DM halo. Panel a: distributions of the median p-values, pmed (in logarithmic scale), obtained by comparing each of the 5000 observed samples with ≃400 HVS generated in a spherical DM halo, with the 5000 mock samples generated in spheroidal DM halos axisymmetric about the z axis with a different qz, as listed in the panel. Panel b: enlargement of the rightmost part of panel a with the pmed axis in linear scale. The different shape of the distributions in panels a and b is due to the different size of the histogram bins in the logarithmic and linear scales.

Panel b is the enlargement of panel a. It shows that the distributions of the median p-values obtained from the comparison of the mock observed samples with the mock samples generated in the spherical DM halo with q_z = 1.0 (the gray histogram) is flatter than the corresponding histogram in panel b of Fig. 11. Indeed, all the distributions of Fig. 17 are flatter than the corresponding distributions of Fig. 11. As a consequence, there is a larger superposition of the gray distribution with the p_med distributions obtained from the KS test comparisons of the observed samples and the samples that crossed the mildly oblate (q_z = 0.9) and mildly prolate (q_z = 1.1) DM halos. Unlike Fig. 11, Fig. 17 also shows a non-null overlap of the gray distribution with the red distribution obtained from the comparison of the observed sample against the sample that crossed an oblate DM halo with q_z = 0.8; however, this non-null overlap does not correspond to any actual erroneous associations of the observed sample with an oblate DM halo with q_z = 0.8, because no comparison yields a p_med of the red distribution higher than the p_med of the gray distribution (see Sect. 5.2.1).

Summarizing, a 50% smaller HVS sample implies a higher probability of assigning a DM halo with an incorrect axial ratio q_z to an observed sample of HVSs. This corresponds to a lower success rate S of the method: while for the original HVS observed sample we obtain S = 98.4%, for the 50% smaller sample we obtain S = 91.8%. However, the axial ratio can be off by |Δq_z| = 0.1 only, as for the full-size sample.

The flattening of the distributions of the median p-values keeps on increasing with decreasing size of the observed sample. Figure 18 shows the change in the distribution of the median p-values obtained from the KS test comparison of the HVS observed samples against the mock samples generated in a spherical DM halo, when the size of the HVS observed samples drops from 100% (gray histogram; N ≃ 800 stars) to 50% (brown histogram; N ≃ 400 stars) and to 5% (black histogram; N ≃ 40 stars). The increasing flattening of the p_med distributions for smaller HVS samples occurs because deviations of the distribution of the p-values from the uniform distribution increase with decreasing size of the HVS sample. In turn, the difference between p_med and the expected value 0.5 increases and the distributions of p_med’s thus have higher tails.

Fig. 18. Distributions of the median p-values obtained from the KS test comparison of the ${\stackrel{\sim }{D}}_{|v_{\vartheta }|}$ of each observed sample, generated in a spherical DM halo, with the 5000 mock samples generated in a spherical DM halo. The size of the compared samples is 100% (gray histogram; N ≃ 800 stars), 50% (brown histogram; N ≃ 400 stars), and 5% (black histogram; N ≃ 40 stars), respectively, of the original HVS sample size.

For the case of a spherical DM halo, the black dots in Fig. 19 show the success rate of our method as a function of the size of the HVS observed sample. A success rate S ≳ 90% requires an HVS sample of N ≳ 400 HVSs, while S ≳ 75% is achieved with N ≳ 200 HVSs, S ≳ 55% is obtained with a sample of N ≳ 80 HVSs, and S ≳ 40% is obtained with a sample of N ≳ 40 HVSs. When the HVS sample size is N ≃ 400, the unsuccessful cases yield axial ratios that can be off by |Δq_z| = 0.1 at most, as it is for the case of the optimal sample size, N ≃ 800. On the other hand, when the sample size is N ≃ 200, the axial ratio can be off by |Δq_z|≤0.2. For the smallest sample sizes considered here (i.e., for N ≃ 80 and N ≃ 40), the axial ratio can be off by |Δq_z|≤0.4. For any sample size, the probability of associating a given erroneous shape with the DM halo crossed by the observed sample decreases with increasing difference |Δq_z| between the erroneous and the actual axis ratio.

Fig. 19. Success rate S of our method in recovering the correct shape of a DM halo as a function of the size of the HVS observed sample, for a spherical DM halo (qz = 1.0; black dots) and for two spheroidal DM halos axisymmetric about the z axis with qz = 1.1 (green dots) and qz = 1.4 (red dots). The sample size reported on the x axis is the average size of the 5000 HVS observed samples used to estimate the success rate in the spherical DM halo scenario. S is shown for sample sizes of 100% (N ≃ 800 stars), 50% (N ≃ 400 stars), 25% (N ≃ 200 stars), 10% (N ≃ 80 stars), and 5% (N ≃ 40 stars) of the original sample size.

7.2. Nonspherical DM halos

For spheroidal DM halos that are axisymmetric about the z axis, the dependence of the success rate S of our method on the sample size is comparable to the dependence found for spherical DM halos, illustrated in Sect. 7.1. However, the shape dependence of S discussed in Sects. 5.2.2 and 6.3.2 is responsible for fluctuations of the value of S for any given HVS sample size.

Figure 19 shows the dependence of the success rate S on the HVS sample size for the spheroidal DM halos axisymmetric about the z axis that yield the most extreme success rates, namely those with q_z = 1.1 (green dots) and q_z = 1.4 (red dots), superimposed to the same dependence for spherical DM halos (black dots). The different S’s between each of the two spheroidal cases and the spherical case are represented by the vertical green and red bars superimposed to the black dots. The spheroidal DM halo with q_z = 1.4 yields the lowest success rate for each sample size, and the difference in S between this case and the spherical case can be as high as ∼25%. On the other hand, the spheroidal DM halo with q_z = 1.1 yields the largest success rate for each sample size, with differences in S from the spherical case always smaller than 10%.

For fully triaxial DM halos or for spheroidal DM halos non-axisymmetric about the z axis, we refrain from performing a detailed analysis: for the optimal sample size of ∼800 HVSs, our method recovers the correct shape of the DM halo with a success rate S = 96%−100% for the explored scenarios, depending on the actual shape of the DM halo (see Sect. 6.3.2). This success rate is always higher than the success rate obtained, for the optimal sample size, for a spheroidal DM halo axisymmetric about the z axis and with q_z = 1.4; for some shapes of the non-axisymmetric DM halo, S can also be higher than that obtained for a spheroidal DM halo axisymmetric about the z axis and with q_z = 1.1. We thus expect that, for the case of non-axisymmetric DM halos, S will be generally larger than for the axisymmetric cases with the same sample size.

8. Discussion

The HVSs are ejected from the Galactic Center and cross a large range of distances during their journey across the Galaxy. As shown in Sect. 3.1, their use as tracers of the DM halo of the MW is appropriate in any region where the DM halo is expected to dominate their kinematics, namely at galactocentric distances r ≳ 10 kpc. However, the phase-space coordinates of the HVSs at any radius r stores the information on the triaxiality parameters of the dark halo within r. Therefore, the HVSs appear to be a promising probe of the triaxiality of the DM halo of the MW over a large range of spatial scales.

Conversely, other observational probes can constrain the shape of the DM halo over more limited spatial scales (see Sect. 1). Our HVS-based method is thus a powerful tool to complement the currently available constraints on the shape of the Galactic DM halo provided by different tracers on different, limited scales.

Clearly, the applicability of the presented method depends on a few working hypotheses, which we discuss below. In addition, the size of the available HVS samples with measured velocity is still an issue both for the success rate of the method and for the magnitude of the offset between the recovered and the actual axis ratio. However, even though our method was developed in the framework of a specific model of HVS production and with a set of assumptions, our working hypotheses do not limit the validity of the method, as we illustrate in the following.

As for the HVS production, we considered here the Hills mechanism, because of its unique ability in generating both a large number of unbound main-sequence HVSs and B-type stars in close orbit about SgrA^⋆ (see Sect. 1). However, our method can be applied to any mechanism that can eject stars from the Galactic Center on a purely radial direction with a velocity v_ej ≳ 730 km s⁻¹, thus enabling them to reach r > 10 kpc with non-null outward radial velocity.

All our simulated HVSs have a mass of 4 M_⊙. The large majority of the HVS candidates are B-type stars with mass in the range ∼2−4 M_⊙, while other candidates are classified as A, F, G, or K type stars. The generation of a realistic sample of ejected stars would require sampling the masses of the binary stars that encounter the SMBH from a mass distribution. However, under the simplifying assumption that the population of binary stars were entirely composed of equal-mass binaries, the distribution of the ejection velocities of a sample of HVSs with different masses would be the sum of the independent distributions of the ejection velocities of equal-mass HVSs. Thus, our simulated mock HVS sample of 4 M_⊙ stars has the same fractional distribution of ejection velocities that would be obtained for the subsample of 4 M_⊙ stars in a simulated extended spectrum of masses for equal-mass binaries. As a consequence, our mock observed sample would mirror a real sample limited to 4 M_⊙ HVSs.

In fact, in a realistic scenario, unequal mass binaries are likely. For hyperbolic encounters, such as those simulated in this work, the primary member of each binary is usually ejected as an HVS; its ejection velocity depends on the mass m₂ of the secondary star. Thus, in a distribution of binaries whose primary member is a 4 M_⊙ star, the velocity of the ejected 4 M_⊙ HVS depends on the mass of the companion star, m₂ ≤ 4 M_⊙. Nevertheless, Bromley et al. (2006) show that the fractional velocity distribution of 4 M_⊙ HVSs located at r = 10−120 kpc is insensitive to m₂ when m₂= 0.5, 1, 2, and 4 M_⊙. Our simulations confirm this result for binaries with m₂ uniformly sampled in the range 0.1−4 M_⊙. The insensitivity of the velocity distribution to m₂ follows from the fact that the stars that reach r > 10 kpc populate the high velocity tail of the distribution of ejection velocities (v_ej ≳ 730 km s⁻¹); the high-speed tail of the ejection velocity distribution obtained in the case of 4 + 4 M_⊙ binaries is statistically indistinguishable from the high-speed tail of the corresponding distribution obtained when m₂ is sampled in the range 0.1−4 M_⊙. Thus, limiting our mock HVS sample to 4 M_⊙ stars ejected from equal mass binaries returns, for stars in the region of interest (r > 10 kpc), final velocity distributions that are statistically consistent with those we would obtain by fixing the mass of the primary star to 4 M_⊙ and by varying the mass of the secondary star in the range 0.1−4 M_⊙. However, the normalization of the latter distributions is smaller than the normalization we obtain by simulating 4 + 4 M_⊙ binaries. Indeed, when m₂ = 0.1−4 M_⊙, the ejection velocity distribution peaks at lower values, implying that fewer HVSs can reach r > 10 kpc.

The size of the HVS sample also decreases with decreasing ejection rate R. By adopting R = 10⁻⁴ yr⁻¹ (see Sect. 2.3), we generated an optimal sample of ∼800 4 M_⊙ stars. We investigated the effect of different sizes of the HVS samples on the method success rate in Sect. 7, and showed that the success rate decreases for smaller sample sizes. This analysis is suggestive of the effect of a lower ejection rate on the success rate of our method. Indeed, the ejection rate depends on a number of assumptions and it is still poorly constrained (see, e.g., Hills 1988; Yu & Tremaine 2003; Bromley et al. 2012; Zhang et al. 2013; Brown 2015). Bromley et al. (2012) assume continuous star formation in the Galactic Center and estimate R ≈ 1−2 × 10⁻³ yr⁻¹ when integrating over all the mass spectrum of the ejected stars considered. Our adopted R = 10⁻⁴ yr⁻¹ for 4 M_⊙ stars appears roughly consistent with this analysis of Bromley et al. (2012). However, by dropping the assumption of continuous star formation, Bromley et al. (2012) find substantial smaller values: R ∼ 2−8 × 10⁻⁵ yr⁻¹. This range partly overlaps with the range ∼10⁻⁵− a few 10⁻⁴ yr⁻¹ found by Zhang et al. (2013) who consider different origins of the injected binaries and different models of the initial mass function of the primary stars (see also Yu & Tremaine 2003).

The decrease of the success rate with the decreasing size of the HVS sample illustrated in Sect. 7 also mimics the effect of the reduction of the HVS sample because of observational limitations, like the star magnitude or the star position within the Galaxy.

The dependence of the size of the HVS sample on the mass of the stars is more intricate: the larger number of stars with M < 4 M_⊙ and lifetime τ_L > 160 Myr would determine a larger number of HVSs that are alive at the observation time compared to our 4 M_⊙ star sample. However, longer-lived, lower-mass stars can reach larger galactocentric distances after experiencing the inner turnaround (Sect. 3.2); therefore, the lower limit r = 10 kpc we adopted here for the 4 M_⊙ stars will be larger for lower-mass stars, thus potentially reducing the number of HVSs suitable for our method. Opposite effects would be determined by the moderate number of stars with M > 4 M_⊙, and lifetime τ_L < 160 Myr. We plan to investigate in future work how these effects combine to determine the optimal sample size and, in turn, the success rate of the method.

Our HVS mock samples include both unbound and bound HVSs, which are HVSs whose ejection velocity does not exceed the escape velocity of the MW. The bound HVSs that satisfy the selection criteria that we defined in Sect. 3.2 are indicators of the shape of the DM halo as good as the unbound HVSs. Therefore, the observation of bound HVSs is of fundamental importance to increase the HVS sample size and the success rate of the method.

As highlighted by Marchetti (2021), bound HVSs may be the best candidate stars ejected from the Galactic Center that can be observed in the Gaia catalogs, while the majority of the unbound HVS population is expected to be too far from the Sun’s position (Kenyon et al. 2014) for its radial velocity to be measured by the Gaia Radial Velocity Spectrometer (RVS). In the current sample of HVS candidates, ∼70% of the stars are located at a galactocentric distance r > 10 kpc, at the 3σ level. The possibility to measure the radial velocities of fainter objects in the outer halo will be of fundamental importance to identify new bound and unbound HVS candidates that satisfy our selection criteria. For example, the forthcoming 4-metre Multi-Object Spectroscopic Telescope (4MOST; de Jong et al. 2019) will be able to increase the volume of the spectroscopic sample provided by Gaia; specifically, it will measure the radial velocities of Gaia photometric sources with magnitude G < 20.5, while the Gaia RVS will provide radial velocities for sources with G_RVS ≤ 16.2 mag, corresponding to G ≤ 15.9 for B0V stars and G ≤ 17.4 mag for K4V stars (Jordi et al. 2010; Fitzgerald 1970).

If our sample of 4 M_⊙ HVSs reached the optimal size of ∼800 stars, in the ideal case of null uncertainties, our method would be able to recover the correct axis ratios of the DM halo in ≳90% of cases, while being off by 0.1 in the remaining ≲10%. The offset of 0.1 found in the minority of unsuccessful cases is set by the resolution in triaxiality parameters that we choose in our simulations, Δq_y, z = 0.1.

Here, we did not investigate the effect of any observational limitations on our HVS observed sample, nor the effect of the uncertainties on the HVS distances and galactocentric tangential velocities on the success rate of our method. This investigation will be pursued elsewhere.

The main contribution to the uncertainty on the tangential velocity of distant HVSs comes from their proper motion. The proper motion measurements currently available in the Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2016b, 2021) are affected by uncertainties that can lead to relative uncertainties on the tangential velocities larger than 100%. On the other hand, for nearby HVS candidates, the largest contribution to the uncertainties on the tangential velocities comes from the star distances inferred from Gaia parallaxes: in the Gaia EDR3, the uncertainties on the parallaxes range from 0.02−0.03 mas, for sources with G < 15 mag, to 1.3 mas, for sources with G = 21 mag (Gaia Collaboration 2021). The relative uncertainties on the parallax-inferred distances of an HVS candidate at 10 kpc from the Sun can thus vary from 20−30% for the brightest stars to 1300% for the faintest stars. To preserve the high success rate of our method, it clearly is of utmost importance to reduce the uncertainties on both the star distances and proper motions.

A future Theia-like mission (Malbet et al. 2016, 2019, 2021; The Theia Collaboration 2017), designed for unprecedented high precision astrometry, may achieve an end-of-mission uncertainty on proper motions of a few micro-arcseconds per year (i.e., ∼100 times smaller than the uncertainty of Gaia), opening up the possibility for significantly constraining the shape of the DM halo. The availability of more precise measurements of proper motions will also make it possible to better constrain the birthplace of the current HVS candidates.

The confirmation of the galactocentric origin of the trajectories of HVS candidates and their use to constrain the Galactic gravitational potential requires a final cautionary note. This research program is far from obvious, and can suffer from a circularity problem. The galactocentric origin of an HVS candidate can be assessed through a backtracking of the star trajectory (e.g., Brown et al. 2014, 2018; Marchetti et al. 2019; Irrgang et al. 2018, 2021; Koposov et al. 2020; Kreuzer et al. 2020); in turn, the backtracking requires an assumption on the gravitational potential that one aims to constrain. In addition to the matter distribution of the Galaxy, the gravitational potential is set by the distribution of the satellites of the Galaxy (e.g., Kenyon et al. 2018). Moreover, the gravitational potential of the Galaxy is time dependent because of the interaction of the Galaxy with its satellites (Boubert et al. 2020). On a wider perspective, the galactocentric origin of an HVS candidate also depends on the theory of gravity (Chakrabarty et al. 2022). Therefore, a solid confirmation of the galactocentric origin is currently limited to high-speed stars that are close to the Galactic Center, where the effects mentioned above are negligible (e.g., Koposov et al. 2020). An appropriate method for constraining the gravitational potential with slower and more distant HVS candidates would require, for example, an iterative approach. In the absence of a self-consistent procedure, constraints on the Galaxy mass model with HVS candidates remain questionable.

9. Conclusions

We have presented a new method for inferring the shape of the DM halo of the MW from the distribution of the components of the galactocentric tangential velocities of a sample of HVSs. We applied our method to an ideal optimal sample of ∼800 4 M_⊙ simulated HVSs. We refer to this sample as the observed sample. We illustrated the method for both axisymmetric and non-axisymmetric Galactic potentials.

In the axisymmetric scenario, we recovered the axial ratio of the DM halo from the one-dimensional distribution of one shape indicator only: the magnitude of the latitudinal velocity, |v_ϑ|, of the HVSs of the observed sample. In the non-axisymmetric scenario, we recovered the axial ratios from the two-dimensional distribution of two shape indicators: |v_ϑ| and a function ${\overline{v}}_{\phi }$ of the azimuthal velocity, v_φ, of the HVSs of the observed sample.

The method is based on the use of the one- or two-dimensional KS test to compare the distribution of the shape indicator(s) of the observed sample, D_|vϑ| or $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$, against the corresponding distributions of HVS mock samples that traveled in DM halos with different shapes. The resolution in axial ratios of our ensemble of mock catalogs is Δq_y, z = 0.1.

For each shape of the DM halo, we compared the observed sample with a set of 5000 D_|vϑ|’s or $D_{|v_{\vartheta }|,{\overline{v}}_{\phi }}$’s, which account for the different ejection initial conditions of the HVSs. A distribution of 5000 p-values was thus obtained for each shape. We used the median of these p-value distributions, p_med, to identify the shape of the DM halo that best matched the shape of the DM halo crossed by the observed sample: the highest p_med value comes from the p-value distribution associated with the correct shape of the DM halo.

In our ideal case of galactocentric velocities with null uncertainties and no observational limitations, the method has a success rate S ≳ 89% in recovering the correct shape of the DM halo of an axisymmetric Galactic potential, and a success rate S > 96% in recovering the correct shape of the DM halo of the non-axisymmetric Galactic potentials that we explored in this work.

The higher success rate of our method for a non-axisymmetric Galactic potential is due to (i) the availability of two shape indicators, |v_ϑ| and ${\overline{v}}_{\phi }$, compared to one shape indicator for an axisymmetric Galactic potential; and (ii) the sensitivity of the azimuthal velocity, v_φ, to the shape of the DM halo. In the small fraction of unsuccessful cases, an erroneous shape association occurs; however, the discrepancy from the correct axial ratio is small. Indeed, in the axisymmetric Galactic potential scenario, the erroneous DM halo axis ratio, q_z, typically differs from the correct ratio by |Δq_z| = 0.1; an offset |Δq_z| = 0.2 is very rare (≲0.04% of cases). In the case of a non-axisymmetric Galactic potential, the incorrect shape associations are expected to be even rarer because our method has a higher constraining power for a non-axisymmetric than for an axisymmetric Galactic potential.

The success rate of our method depends on the size of the HVS sample and is higher for larger HVS samples. In the case of an axisymmetric Galactic potential, a decrease in the sample size corresponds to a decrease in the success rate that depends on the actual shape of the DM halo: for a spherical halo, a decrease in the sample size from 800 to 40 mock observed HVSs implies a drop in S from ∼98% to ∼38%; for a spheroidal halo axisymmetric about the z axis with q_z = 1.1, which yields the highest S, S drops from ∼99% to ∼41% for the same decrease in the sample size; for a spheroidal halo with q_z = 1.4, which yields the lowest S, S drops from ∼89% to ∼32%. On the other hand, for any non-axisymmetric Galactic potential, and for any given sample size, S is expected to always be larger than the rate obtained in the most unfavorable axisymmetric case (q_z = 1.4) and typically also larger than the rate obtained in the most favorable axisymmetric case (q_z = 1.1). In addition to the increase in the success rate, S, increasing the size of the HVS sample decreases the discrepancy between the inferred shape and the correct shape.

Accurate estimates of the success rate of our method when applied to real data require the generation of more realistic mock observed HVS samples that account for (i) the uncertainties on the distances and velocities of the HVSs, (ii) the observational limitations of the HVS sample, and (iii) an appropriate mass distribution of the ejected stars. Nevertheless, our analysis suggests that a robust determination of the shape of the DM potential requires measuring the galactocentric velocities of a few hundred HVSs whose galactocentric origin is robustly confirmed.

We will assess elsewhere the sensitivity that is required for our method and that might be reached with future astrometric space missions (Malbet et al. 2016, 2019, 2021; The Theia Collaboration 2017). Similarly, we will investigate the success rate of our method with expected realistic HVS samples that might come from upcoming spectroscopic surveys (de Jong et al. 2019).

Acknowledgments

We sincerely thank Scott Kenyon, Margaret Geller, and Warren Brown for numerous inspiring discussions that stimulated our initial interest in the dynamics of HVSs. We are in debt with the referee whose sharp comments contributed to clarify the presentation of our work. The Departments of Excellence grant L.232/2016 of the Italian Ministry of Education, University and Research (MIUR) fully supported the PhD fellowship of AG and partly supported the fellowship of SSC. SSC is also supported by the grant “The Milky Way and Dwarf Weights with Space Scales” funded by University of Torino and Compagnia di S. Paolo (UniTO-CSP), by the grant no. IDROL 70541 IDRF 2020.0756 funded by Fondazione CRT, and by INFN. 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 research has made use of NASA’s Astrophysics Data System Bibliographic Services, of the routines of the NUMERICAL RECIPES (Press et al. 2007), of TOPCAT (Taylor 2005) and of the NUMPY (Harris et al. 2020) and MATPLOTLIB (Hunter 2007) PYTHON modules.

References

Appendix A: Mock catalogs and observed sample

Figure A.1 shows the distributions of the distances to the Galactic Center, and of the galactocentric radial, latitudinal, and azimuthal velocities of our simulated HVSs at the steady state (i.e., at the observation time t_obs; see Sect. 2.3) for a mock catalog (see Sect. 2.4) generated in a gravitational potential where the DM halo is triaxial with axis ratios (q_y, q_z) = (1.2, 0.9). The gray histograms correspond to the full sample of simulated HVSs, while the green histograms represent our observed sample of ∼800 HVSs, namely the ideal optimal sample of HVSs that we derived from the full sample by applying the selection criteria r > 10 kpc and v_r > 0 (see Sects. 3.1, 3.2, and 4). The distributions of the kinematic properties of the stars in the other Galactic potentials we investigated in this work are qualitatively similar to those shown in Fig. A.1.

Fig. A.1. Kinematic properties of the HVSs of one of our mock catalogs, corresponding to a Galaxy whose DM halo has a triaxial gravitational potential with (qy, qz) = (1.2, 0.9). The gray histograms represent the full sample of HVSs; the green histograms represent our observed sample, derived from the full sample by requiring r > 10 kpc and vr > 0. Top left panel: Distribution of galactocentric distances; the vertical dashed line corresponds to our threshold radius r > 10 kpc. Top right panel: Distribution of the radial velocities; the vertical dashed line corresponds to vr = 0. Bottom left panel: Distribution of the latitudinal velocities. Bottom right panel: Distribution of the azimuthal velocities.

Appendix B: Transformation of coordinates and velocities

The Galactic heliocentric position of a star is (d, ℓ, b), with d the heliocentric distance to the star, ℓ its Galactic longitude, and b its Galactic latitude; the Galactic heliocentric velocity is (v_d, v_ℓ, v_b), with v_d the radial velocity of the star, and v_ℓ and v_b the longitudinal and latitudinal components of the heliocentric tangential velocity. The Galactic heliocentric position and velocity components are

$$\begin{array}{c}\hfill d=\sqrt{{(x+R_{\odot })}^2+y^2+z^2},\end{array}$$(B.1)

$$\begin{array}{c}\hfill l=\arctan \left(\frac{x\tan \phi }{x+R_{\odot }}\right),\end{array}$$(B.2)

$$\begin{array}{c}\hfill b=\arcsin \left(\frac{z}{d}\right),\end{array}$$(B.3)

and

$$\begin{array}{cc}\hfill v_d=& (v_x-U_{\odot })\cos l\cos b+[v_y-(V_{\odot }+{\mathrm{\Theta }}_0)]\sin l\cos b\hfill \\ \hfill & +(v_z-W_{\odot })\sin b,\hfill \end{array}$$(B.4)

$$\begin{array}{c}\hfill v_l=[v_y-(V_{\odot }+{\mathrm{\Theta }}_0)]\cos l-(v_x-U_{\odot })\sin l,\end{array}$$(B.5)

$$\begin{array}{cc}\hfill v_b=& (U_{\odot }-v_x)\cos l\sin b-[v_y-(V_{\odot }+{\mathrm{\Theta }}_0)]\sin l\sin b\hfill \\ \hfill & +(v_z-W_{\odot })\cos b.\hfill \end{array}$$(B.6)

The longitudinal component of the proper motion μ is μ_ℓ = v_ℓ/(d cos b), and its latitudinal component is μ_b = v_b/d.

We transformed the Galactic heliocentric star position and proper motion to the equatorial system. In this system, the star position is (d, α, δ), with α the right ascension and δ the declination; the star velocity is (v_d, v_α, v_δ), with v_d still the radial velocity of the star, and v_α and v_δ the components of the star tangential velocity along the right ascension and declination, respectively; the star proper motion is μ = (μ_α, μ_δ), with μ_α = v_α/(d cos δ), and μ_δ = v_δ/d.

For the star position, we adopted the transformation equations (Duffett-Smith 1979)

$$\begin{array}{c}\hfill \alpha =\arctan \left[\frac{\cos b\cos (\ell -{\ell }_{\mathrm{asc}})}{\sin b\cos {\delta }_G-\cos b\sin {\delta }_G\sin (\ell -{\ell }_{\mathrm{asc}})}\right]+{\alpha }_G,\end{array}$$(B.7)

$$\begin{array}{c}\hfill \delta =\arcsin [\cos b\cos {\delta }_G\sin (\ell -{\ell }_{\mathrm{asc}})+\sin b\sin {\delta }_G],\end{array}$$(B.8)

where α_G and δ_G are the equatorial coordinates of the North Galactic Pole, and ℓ_asc is the Galactic longitude of the ascending node of the Galactic plane, related to the Galactic longitude of the North Celestial Pole by ℓ_asc = ℓ_NCP − 90°. We chose J2000.0 as a reference point for the coordinate conversion. This choice yields α_G = 192.85948123°, δ_G = 27.12825120° (Cox 2000), and ℓ_asc = 32.93192° (Poleski 2013).

For the components of the proper motion μ, we adopted the transformation equations (Poleski 2013)

$$\begin{array}{c}\hfill {\mu }_{\delta }=\left(\frac{C_2\cos b}{C_1^2+C_2^2}\right)({{\mu }_{\ell }}^{\ast }+\frac{C_1}{C_2}{\mu }_b),\end{array}$$(B.9)

$$\begin{array}{c}\hfill {\mu }_{\alpha }^{\ast }=-\frac{1}{C_1}(C_2{\mu }_{\delta }-{\mu }_{\ell }\cos b),\end{array}$$(B.10)

with

$$\begin{array}{cc}\hfill C_1& =\sin {\delta }_G\cos \delta -\cos {\delta }_G\sin \delta \cos (\alpha -{\alpha }_G),\hfill \end{array}$$(B.11)

$$\begin{array}{cc}\hfill C_2& =\cos {\delta }_G\sin (\alpha -{\alpha }_G),\hfill \end{array}$$(B.12)

and ${\mu }_{\mathcal{l}}^{*}={\mu }_{\mathcal{l}}\text{cos}b$, ${\mu }_{\alpha }^{*}={\mu }_{\alpha }\text{cos}\delta$, and ${\mu }_{\alpha }^{\ast 2}+{\mu }_{\delta }^2={\mu }_{\ell }^{\ast 2}+{\mu }_b^2={\mu }^2.$

All Tables