On measuring the Galactic dark matter halo with hypervelocity stars

Hypervelocity stars (HVSs) travel from the Galactic Centre to the dark matter halo of the Milky Way, where they are observed with velocities in excess of the Galactic escape speed. Their trajectories make them a unique probe of the still poorly constrained dark matter component of the Galaxy. In this paper, we present a new method to constrain the Galactic potential with HVSs. The likelihood is constructed by efficiently calculating the local HVS density at any point of the Galaxy by back-propagating the phase space position and quantifying the ejection probability along the orbit. This method is particularly suited to the data from the ESA mission {\it Gaia}. Therefore, to showcase our method, we applied it to a simulated \Gaia sample of $\sim200$ stars in Galactic potentials with three different dark matter components, parametrized by a spheroidal NFW profile. We find that individual HVSs exhibit a degeneracy in the scale mass-scale radius plane ($M_s-r_s$) and are able to measure only the combination $\alpha = M_s/r_s^2$, likewise a degeneracy is also present between $\alpha$ and the spheroidal axis-ratio $q$. When the whole sample is considered, both parameters are nailed down with {\it sub-percentage} precision (about $1\%$ and $0.1\%$ for $\alpha$ and $q$ respectively) and no systematic bias is observed. This remarkable power to constrain deviations from a symmetric halo is a consequence of the Galactocentric origin of HVSs. To compare our results with other probes, we break the degeneracy in the scale parameters and impose a mass-concentration relation. The result is a competitive precision on the virial mass $M_{200}$ of about $10\%$.


INTRODUCTION
In the concordance ΛCDM model of cosmology, galaxies are embedded inside larger structures known as haloes. These are made of a dissipationless fluid called dark matter, visible only through its gravitational effects. Over cosmic time, haloes grow in mass and size through hierarchical clustering, starting from the initial perturbations of a slightly inhomogeneous matter density field. Despite its central role in structure formation, the nature of dark matter and its microscopic physics are still unknown (see, e.g., Garrett & Duda 2011).
There is a number of theoretical predictions associated to the shape and mass of dark matter structures. Pure cold dark matter simulations suggest that collapsed haloes acquire a triaxial ellipsoid shape, but more recently it has been found that the inclusion of baryonic matter results in rounder shapes (e.g., Debattista et al. 2008). Similarly, self E-mail: contigiani@strw.leidenuniv.nl interacting dark matter is also expected to induce spherical haloes in the innermost regions (Peter et al. 2013). Furthermore, measurements of the Milky Way's halo, together with observations of surrounding dwarf galaxies, can be used as a test for the concordance model (Moore et al. 1999;Klypin et al. 1999). For example, a total mass of the Milky Way lower than 10 12 M can align the observed number of satellite galaxies with what is predicted in simulations (Wang et al. 2011).
Gravitational lensing is the most common technique used to measure the dark matter distributions of statistical samples of distant galaxies and galaxy groups (e.g., Hoekstra et al. 2013;Mandelbaum 2014). In the case of the Milky Way, our privileged position mandates the use of a different set of techniques and dynamical tracers are employed to measure the structure of its dark matter halo. Objects traveling through the halo act as test particles subjected to its gravitational potential and their trajectories in phase space can be traced to constrain a parametric model for the density profile. This procedure usually requires assumptions about the initial conditions or the steady-state configuration of the system.
In the Galactic bulge and disc, where baryons dominate the matter density, established techniques based on the kinematics of field stars or HI emission are used (Portail et al. 2017;Reid & Dame 2016). Unfortunately, the scarcity of these tracers outside the Galactic disc limits their constraining power where the dark matter halo dominates (e.g., Huang et al. 2016). Since its discovery (Newberg et al. 2002), the Sagittarius stellar stream has proved to be a valuable dynamical tracer (e.g., Law et al. 2009;Deg & Widrow 2013;Gibbons et al. 2014). The Sagittarius dwarf galaxy is one of the closest satellites to the Milky Way and it is in the process of being tidally disrupted. The strong tidal forces give rise to a long stream of tidal debris which orbits the Milky Way. Other tidal streams have been discovered over the years: some of them connected to globular clusters (e.g. Palomar 5, Odenkirchen et al. 2001) and some others represent the last remnants of now defunct dwarf galaxies (e.g. Virgo, Duffau et al. 2014).
Despite the existence of these multiple tracers, there is no consensus in the literature on the mass of the Milky Way halo (Wang et al. 2015): measurements differ up to a factor 5 and relative precisions range from below 10% to roughly 100%. The situation is no different when, instead of its mass, the halo shape is considered. While the Milky Way's dark matter halo is often measured to be a spheroid with two of its axes being equal and aligned with the disc galaxy within (e.g., Bovy et al. 2015;Pearson et al. 2015), conflicting measurements are present in the literature and triaxial shapes have also been suggested (e.g. Law & Majewski 2010). The halo shape could also be a function of radius, spheroidal in the centre and triaxial in the outer region (Vera-Ciro & Helmi 2013). In the case of a pure spheroid, the ratio between the third axis and one of the others is usually referred to as c/a or, like in this paper, just q. A ratio q = 1 corresponds to a sphere, while the conditions q > 1 and q < 1 correspond respectively to a prolate or an oblate spheroid. Reports range from a spherical halo (e.g., Bovy et al. 2016b) to oblate (e.g., Loebman et al. 2014) or prolate (e.g., Bowden et al. 2016;Posti & Helmi 2018). It is clear that when previous endeavours to measure the Galactic halo are put together, the tensions between different probes imply the existence of systematic biases.
In future years, hypervelocity stars (HVSs) are expected to be introduced to this landscape as a powerful probe. For the purposes of this work we will refer to HVSs as high velocity objects (Galactocentric velocity > 450 km/s) travelling from the Galactic Centre (GC) along quasi-radial orbits. In 2005 the first HVS was discovered (Brown et al. 2005): a Btype main sequence star with radial velocity in the Galactic rest frame of about 700 km/s. Subsequent observations have measured its distance from the GC, found to be of the order of 100 kpc (Brown et al. 2014). Given its high velocity, the object was measured to be unbound form the Galaxy. Over the years, objects with similar stellar properties have been found and to date the largest and most studied sample is composed of the 21 HVS candidates reported by Brown et al. (2014), a survey targeting B-type stars in the outer halo. In the near future, the high-quality sample of HVSs predicted to be observed by the satellite mission Gaia (Gaia Collaboration 2016; Gaia Collaboration et al. 2018) by the early 2020s is expected to contain several hundred objects (Marchetti et al. 2018) and will offer a new diffuse dynamical tracer for the Galactic potential.
The goal of this paper is to introduce a new method to exploit this tracer. Gnedin et al. (2005) already showed that a few HVSs can be a powerful tool to constrain the shape and orientation of the Galactic halo and a precision of about 10% can be reached if accurate proper motion and Galactocentric distances are known. Later, Yu & Madau (2007) have shown how the triaxiality of an ellipsoidal halo can be estimated directly from observed HVS positions and velocities under a specific halo model. Other similar attempts include Perets et al. (2009), who explored how asymmetries in the radial velocity distribution of halo stars due to HVSs depend on the Milky Way mass, and Fragione & Loeb (2017), which is an application of such method. In other cases, inferences about the Galactic gravitational potential behind the deceleration of HVSs assume a certain class of ejection velocity distributions (e.g., Sesana et al. 2007;Rossi et al. 2017).
We expand on previous works by developing a new versatile technique that can be adapted with minimal assumptions to a variety of models for the ejection mechanism and Galactic potential. This is of the uttermost importance to produce unbiased joint constraints in combination with other probes (see Rossi et al. 2017, where two of us have shown the power of this approach). Our method is based on a reconstruction of the HVS orbital history and it has the advantage of not requiring simulations of the entire population for every potential/ejection model studied. 1 To validate our method we will focus here on HVSs ejected through one realization of the Hills mechanism (Hills 1988;Yu & Tremaine 2003;Sari et al. 2010). According to this mechanism, the three-body interaction between a binary system and a massive black hole (MBH) results in one star orbiting closely around the black hole and the other one being ejected at high velocity. The aforementioned observations of high-velocity stars in the Galactic halo are consistent with the existence of such mechanism and at present, it is still considered the leading explanation (Brown 2015;Brown et al. 2018). Note also that HVSs are expected to be an observational consequence of the massive black hole located in the GC (Ghez et al. 2003).
In Sec. 2 we construct mock populations of this sample, based on previous work (Rossi et al. 2014Marchetti et al. 2018). Our mock catalogues are based on the expected astrometry and photometry of the final Gaia data release. Afterwards, we lay the foundations of our technique and we arrive in Section 3 at an integral formula for the phase space distribution of these objects, which allows us to write down a likelihood function for an observed sample of HVSs. In the same section, we also discuss the advantages and limitations of the method. In Sec. 4 we then test our approach and try to recover the dark matter halo inside which the simulated sample was propagated. In the same section, we also address issues related to the practical implementation of our technique for a Gaia -like sample.

SIMULATED HVS CATALOGUES
The first step to verify how and if HVSs can constrain the dark matter halo of the Milky Way is to produce an observational mock catalogue of HVSs. To produce such sample we need to specify three important ingredients: 1) the ejection distribution that determines how the velocities, positions and masses of our stars are distributed at the moment of ejection from the GC; 2) a survival function that dictates the fraction of HVSs alive after a time t post-ejection, and 3) a gravitational potential under the influence of which the stars trace their orbits. In the next three subsections we present an implementation of these quantities and we conclude, in the last subsection, with the details of our numerical simulation.

Ejection rate distribution
We aim to parametrize the distribution of velocities, positions and masses at ejection for HVSs generated through the Hills mechanism (Hills 1988) by writing down an explicit expression for an ejection rate distribution R(w ), which has the units of a configuration space density per unit time. We call w our configuration space coordinate, w = (x, v, m), where (x, v ) is the usual phase space coordinate (position, velocity) and m is the stellar mass. We follow the set up first described in Rossi et al. (2017) and then implemented by Marchetti et al. (2018).
In a reference system centred on the massive black hole (or equivalently, the GC) we can write: where we have introduced the ejection rate per unit time Γ and the δ terms are Dirac deltas. In this work we will not assume any value for Γ and we will normalize all of the other functions appearing in this expression to unity. The main prediction of the Hills mechanism quantifies the asymptotic velocity of the ejected objects at an infinite distance from the massive black hole. In practice, this distance can be modelled as the radius of the gravitational sphere of influence of the black holer, defined as the radius of the sphere centred on the black hole and containing twice its mass. For distances larger than its radius the potential of the black hole becomes a negligible fraction of the total Galactic potential. We pick the valuer = 3 pc (Genzel et al. 2010) and impose this to be the ejection radius through the Dirac delta function δ (|x |).
The term R H (|v |, m) quantifies the relative probabilities of different initial velocities and masses of HVSs. It can be computed using Monte Carlo (MC) simulations as done in Rossi et al. (2014Rossi et al. ( , 2017; Marchetti et al. (2018). In the first paper it is also shown that the resulting distributions can be easily fitted with analytic functions. By fitting the Hills mechanism MC catalogue in Marchetti et al. (2018) to the functional form suggested by Rossi et al. (2014) we obtain (2) v 0 (m) = 1530 (M /m) 0.65 km/s.
Notice that the velocity distribution for a fixed value of m has a high velocity tail starting from the value v 0 (m).
The last term in eq. 1 is a Dirac delta function imposing zero angular momentum. This condition must be satisfied at any ejection distance |x | >r since every HVS is a product of a close encounter of the progenitor binary with the back hole at a much closer distance. Assuming this distance to be tidal disruption radius r bt , for a massive black hole of mass M bh = 10 6 M and a binary with semi-major axis a ∼ 1 R and total mass m * ∼ 1 M we get r bt = a(M bh /m * ) 1/3 3 pc.

Survival function
If there is no preferred time of ejection, the flight time t f of an HVS of mass m is sampled according to where 1 , 2 are random variables uniformly distributed between 0 and 1, and t L (m) is the stellar lifetime (Marchetti et al. 2018). In our implementation this is taken to be equal to the main sequence lifetime, modelled according to Hurley et al. (2000).
The probability density function of the variable t f is found to be equal to Note that the average value of t f /t L is then expected to be 0.25, i.e. on average HVS fly for a quarter of their lifetime. The function g(t f , m) is then the corresponding survival function: for t f < t L (m).

Galactic potential
We model the Milky Way gravitational potential as the sum of four components: central black hole, bulge, disc and dark matter halo. Depending on the symmetry, we use Cartesian (x, y, z), spherical (r, θ, φ) or cylindrical coordinates (R, ϕ, z).
In all three cases we position the GC at the origin and the z axis perpendicular to the Galactic disc. The first component is a simple Keplerian potential and it is meant to describe the massive black hole at the centre of the Galaxy with a mass of M bh = 4 × 10 6 M (Eisenhauer et al. 2005;Ghez et al. 2008), The second and third components are an Hernquist spheroid (Hernquist 1990) and a Miyamoto-Nagai disc (Miyamoto & Nagai 1975) respectively. This form and model parameters are chosen because they are commonly used to parametrize the baryonic components of the Galactic potential in similar studies (e.g., Johnston et al. 1995   The last component of the potential is a spheroidal NFW density profile (Navarro et al. 1997) and it models the dark matter halo of the Milky Way, Notice that for the sake of simplicity, in our parametrization q corresponds to the dimensionless axis ratio of our spheroid. The potential associated to this matter density is found by solving Poisson's equation, In this study we will focus on three different fiducial Galactic haloes (see Table 1). For model A the chosen values for M s , r s are the best fit parameters to the rotation curve of the Milky Way for a spherical halo . In model B we consider an oblate spheroid as variation of this model, and in model C we consider a spherical halo with a significantly different scale radius and mass.

Mock catalogues
We follow a procedure similar to the one detailed in Marchetti et al. (2018) to generate mock HVS Galactic populations ejected through the Hills mechanism and simulate the effect of the Gaia selection function. In the same paper, we estimated the current Galactic population of HVSs produced through the Hills mechanism to include 10 5 members. An ejection sample of this size is therefore generated by sampling the distribution R(w ) in eq. 1 using a Markov chain Monte Carlo method implemented through the python library emcee  (2010)). The sample is then propagated numerically by the software galpy (Bovy et al. 2015), through the three fiducial models of the Galactic potential presented in Sec. 2.3, using a time step δt = 0.01 Myr. The integration time, i.e. the flight time, for each star is dictated by the formulas presented in section 2.2. At the end, the photometric properties of the stars are simulated using the stellar models provided by Hurley et al. (2000), the BaSeL SED Library 3.1 (Westera et al. 2002) and a map of the dust reddening in the Milky Way presented in Bovy et al. (2016a) (a combination of Drimmel et al. 2003;Marshall et al. 2006;Green et al. 2015). The magnitude in the Gaia band GRVS is then computed using the polynomial fitting functions provided by Jordi et al. (2010).
From this catalogue we define a golden sample by imposing two conditions. First, only stars brighter than the 16th magnitude in the G RV S band are selected. This cut filters objects for which Gaia is expected to measure the line of sight velocity (Cacciari et al. 2016;Katz et al. 2018). The second condition that we impose is related to the velocity of the objects appearing in the sample: we impose a total velocity at present time in the Galactic reference frame higher than 450 km/s. This threshold filters objects which will be clearly recognizable as high velocity -i.e. faster than three times the one dimensional Galactic velocity dispersion (Battaglia et al. 2005;Brown et al. 2010;King III et al. 2015). At the end of this selection process our golden samples contain 195, 192, 211 objects for halos A, B, C respectively. 2 Additional information about the catalogue is available in Marchetti et al. (2018).

DISTRIBUTION FUNCTION
We study the distribution of HVSs in the configuration space is the usual phase space coordinate (position, velocity) and m is the stellar mass. We then introduce the density function of HVSs in this space at a time t: In this expression, dN(t) represents the number of HVSs in the volume d 3 v d 3 x dm. We now aim to write down this distribution as a combination of two other functions: the ejection rate distribution R(w ), which parametrizes the density of HVSs ejected at a given position of the configuration space per unit time, and the survival function g(t, m), which quantifies the fraction of stars of mass m which survives for at least a time t after ejection. In Sec. 2.4 we have provided two example of how these functions might be defined. Notice that we assume a stationary process for the creation of HVSs, meaning that R(w ) is not a function of time.
These definitions allow us to write down the total number of HVSs present in the Galaxy at a time t after the formation of the Milky Way or, equivalently, when the first HVS was ejected: In this expression we integrate R(w ) over the entire configuration space and over every possible ejection time t . In the last integral, the weight function g(t − t , m) accounts for the fact that not all stars ejected at a time t will still be alive after a time t − t . From the expression for N(t) we can derive the density function by applying a Dirac delta function in configuration space: In this expression we introduced the solution of the equations of motion in configuration space, W (w , t ; t), which maps the initial condition w at a time t to the phase space position W at a time t > t . The delta function imposes that objects in the position w at a time t must have been generated inside the appropriate orbit W (w , t ; t) at the appropriate time. Note that if we assume that the stellar mass is not a function of time, Liouville's theorem ensures that the map (w , t ) ↔ (W , t) is bijective and conserves the volume d 7 w. Because of this, applying the Dirac delta in the integral over d 7 w does not introduce a Jacobian term despite the argument of the Dirac delta not being a trivial function of w . Furthermore, because Hamilton's equations for a single HVS are time invariant, we can write In conclusion, we derive the following: In this expression we have introduced the trajectory w (w ; t − t ) which is a solution of the argument of the delta function in eq. (14) and it can be found by integrating numerically back in time the equations of motion from the starting point w .
Notice that this final result is completely general: it does not discriminate between bound or unbound objects and can be applied to a variety of ejection mechanisms and lifetime models.
In this analysis we are interested in exploring how the distribution in eq. (15) is affected by the Galactic potential. The dependence on the dynamics is not made clear from the expression itself, but it is hidden in the backwards trajectory w (w ; t − t ) . If we model the Galactic potential using a set of parameters θ, we can write down the parametric configuration space distribution as: We can then assign for every value of this parameter vector a likelihood to the observation of N HVS HVSs in the configuration space points {w 1 , . . . , w N HVS } at a time t: While the likelihood function formally depends on the observations, in order to simplify the notation our expression does not make this dependence of L on w i explicit.
Our implementation is strictly a forward-fitting algorithm, meaning that it does not produce model-independent results, but it can be used to constrain any parametric model. The first obvious advantage of this technique is that it allows us to parametrize (hence fit) any aspect of the HVS population. For example, we could easily use an observed sample to constrain a parametric version of R(w ). In this case, we would write the dependence on model parameters explicitly into its expression. Notice however that, in order to compare different ejection mechanisms, the rates Γ should be fixed or at least be left as free parameters. Second, we stress that the technique described here can be implemented for unbound and bound trajectories alike. The periodicity is not an issue thanks to the explicit time dependence of g(t f , m). The presence of this function in the integral also means that the time integration should be performed only between now and a time t L (m) in the past, since g(t f , m) is zero by design after this point. Third, because the stars are tested individually and not as a sample, a single one is able to rule out any Galactic potential not consistent with Galactocentric origin or any ejection model unable to reproduce the range of allowed initial velocities.
In practice, the evaluation of L(θ) is performed by integrating numerically back in time the HVS orbits from the observed positions w i under the influence of the potential specified by θ. For consistency, our set-up matches the one employed for the creation of the mock catalogue. Given the orbit w (t −t ) as a function of the backwards time coordinate t we can then evaluate the integral in eq. 15 in the configuration space volume where R(w ) is non-zero. Because of the presence of the Dirac deltas this volume is formally a 4−d space embedded in the 7-d configuration space. To perform the integral and account for numerical errors we swap the Dirac deltas with Gaussian kernels calibrated against the numerical precision of the orbit back-propagation code. This introduces two smoothing parameters which correspond to σ r = 10 pc and σ L = 10 pc×(km/s).

LIKELIHOOD FUNCTION
To test our method, we study the likelihood L(θ) for the golden sample of HVSs simulated in Sec. 2 as input data and the halo potential parameters as variable θ. In Sec. 4.1 we explore the parameter space θ = (M s , r s ), while keeping q fixed at the fiducial value; and in Sec. 4.2 we assume θ = (q) and freeze M s , r s . We then discuss in Sec. 4.3 the implications for the full parameter space θ = (M s , r s , q). We use the subscript 0 (e.g. q 0 ) to indicate the fiducial value of our halo parameters.
Our analysis also helps us identifying which stars are particularly suited to measure the Galactic halo. In Sec. 4.4 we characterize this sample and discuss the impact of observational errors on the orbital reconstruction method we apply to it.

Likelihood in M s − r s plane
For the three haloes A, B, C we evaluate the likelihood, eq. (17), in the space r s , M s using a coarse grid of size 27 × 27. Based on the results, we define three classes of HVSs in our golden sample: strong, average and poor constrainers. This classification is based on the number of points on our grid with non-zero likelihood. Fig. 1 and Fig. 2   significantly different trend of the latter two classes for halo A. The strong constrainers (not shown) are stars for which no particular trend in the likelihood was identified and have non-zero likelihood only in the fiducial model.
For every star we call n the number of non-zero likelihood points associated to it: strong constrainers have n = 1 and average constrainers have 1 ≤ n ≤ 300. The value 300 is picked from visual inspection of the individual likelihoods.
From Fig. 1 we infer that HVSs are sensitive exclusively to the parameter in the M s − r s plane. Following the established notation, we call α 0 the fiducial value of this parameter. Notice that, for s . The top histogram shows the distribution of n for the stars in our golden sample. The relation saturates for n 30, when grid effects start to hinder the estimate of σ α . Therefore, the relation σ α (n) is calibrated using only points outside the shaded area. a spherical NFW potential, this degeneracy is natural and every α corresponds to a value of the local force at small radii: where we expanded the integral around r/r s = 0. A simple physical interpretation of this degeneracy is that the innermost region of the halo is responsible for the majority of the deceleration experienced by these HVSs. For any single star we can interpolate the likelihood in our coarse grid and obtain an estimate of the 1σ error on α associated to it, which we call σ α . Fig. 3 shows how n and σ α are related to each other. The scaling σ α ∝ √ n is indicative of the fact that a constant α represents a 1d curve in the 2d M s − r s plane. After confirming the absence of bias in the measurement of α for the individual stars, we estimate the 1σ error of the stacked likelihood by assuming normality and using the geometric mean of the individual variances: where the proportionality constant for σ α (n) is fitted independently for our three halo models and, for halo A, it is presented in Fig. 3. The variable n i represents the n corresponding to the i−th star. The proportionality constant in σ α (n) is found to be equal to (1.7, 1.2, 3.3) × 10 7 M k pc −2 for halo A, B, C respectively. Table 2 reports the number of sources in each class for our three fiducial haloes and the effective precision in α expected from the strong and average constrainers using this method. Notice that the value of the combined errobarσ α,q 0 is dominated by low-n stars, meaning that we are extremely susceptible not only to the inferred σ α (n) but also to changes in the distribution of the variable n. To mitigate this effect, theσ α /α 0 mentioned in the table does not use an extrapo-  lated σ α (n) to values n < 30, but assumes the constant value σ(n = 30).
In Sec. 4.4 we discuss how various orbital properties strongly correlate with the likelihood classification and how this information can be used to guide future detections of HVSs.

Likelihood in q
Similarly to what we did for the parameter α, we evaluate the likelihood in eq. (17) by varying the parameter q, while fixing the values of M s and r s to their fiducial values (Table 1). We develop again a classification based on the number of non-zero likelihood points, which is shown in Fig. 4. Notice that while the poor constrainers might prefer the fiducial model, we confirm that their individual likelihoods are either extremely broad or significantly biasedsometimes excluding the fiducial value at the 3σ level.
We stress that the labels we attach to the HVSs (either poor, average or strong constrainer) are independent statements for the two parameters α and q. We find, however, significant overlap between them: for halo A, among the 142 strong constrainers for q, 126 are in the average category for α and 15 are in the strong one. In fact, the performance of every star for q is always equal or better than for α. This implies that HVSs are more sensitive to one parameter than the other in our scheme. This is expected; while the parameters M s and r s set the deceleration, a incorrect parameter q can disrupt the ejection point of quasi-radial orbits by introducing additional torque.
For each star, we can relate the number n of non-zero likelihood points for the parameter q to the expected confidence interval σ q . Fig. 5 shows how the two are related and provides the distribution of the values of n for the strong and average constrainers for our halo A model (the same trend is observed in all models). As discussed in Sec. 4.1 for the parameter α, this scatter plot also fixes the proportionality constant for the stacked uncertainty: Notice that this time σ q (n) ∝ n, since our mesh for this section is constructed on the space of the parameter q directly. As before, because we are extremely susceptible to our reconstruction of σ q (n), we do not extrapolate σ q (n) below the value n < 3, but we assume a constant value. Notice how bias notwithstanding, some of the poor constrainers can still be used to calibrate σ q (n). The proportionality constant in σ q (n) is found to be equal to (1.8, 2.3, 2.3) × 10 −3 for halo A, B, C respectively.
Our results are summarized in Table 2, where we present the estimated precisionσ q for the combination of our average and strong constrainers.

Correlation between α and q
In Fig. 6 we show how many stars allow a certain halo model parametrized by α and q. To generate this figure, we have explored the whole parameter space θ = (M s , r s , q) only for 15 stars in the average category for both α and q, as defined in the previous two subsections. We consider only this subset because these stars have a broad likelihood in both projections and are particularly suited to show the presence of correlation.
A correlation is clearly visible for every star, but while in the plane M s − r s they all constrain the same combination α, the same is not true in the plane q − α. As an example of this, in Fig. 6 we also show the degeneracy stripe for 2 stars. Because of this, we expect both direction and size of the combined constraints to depend on the particular selection bias of our sample.
To give an estimate of the impact of this correlation on Table 2. Number of HVSs with different constraining power for the fiducial haloes considered in this work and predictions for the combined relative errorsσ α /α 0 andσ c /α q on the NFW effective parameters α = M s /r 2 s , q. The lower bound on the error σ q is found in Sec. 4.2 (4.1) by fixing M s , r s (q) to the fiducial values and exploring only the direction q (plane M s − r s ). The upper bound is found in Sec. 4.3 after estimating the correlation coefficient between the two parameters.  Number of HVSs Figure 6. Number of HVSs propagated in halo A for which a given halo, parametrized by the effective scale parameter α = M s /r 2 s and shape parameter q, is allowed. The peak at (1, 1) marks the fiducial values for halo A. The figure was created using stars for which a visible spread in the likelihood is present in our grid and the result proves that there is correlation between the two parameters. The contours are created by linearly interpolating the values found on a grid. To illustrate the origin of this degeneracy, the dashed and solid lines delimit regions allowed by two particular stars.
our reconstructed errors we assume the combined likelihood to be a bivariate normal distribution. Notice that in the previous two sections we verified that the two one-directional log-likelihoods log L(α, q 0 ) and log L(α 0 , q) are both normal. However, the 1σ error barsσ α,q 0 andσ q,α 0 we found in Sec. 4.1, 4.2 do not correspond to the standard deviations of the full log L(α, q) in the presence of correlation. A bivariate log-likelihood up to constant terms can be written as: Whereσ α ,σ q are the standard deviation for the two parameters and −1 < ρ < 1 is the correlation coefficient. From this expression it is clear that the standard deviations found in Sec. 4.1, 4.2 are an underestimate of the real error bars in the full α, q parameter space and should be multiplied by a factor (1 − ρ 2 ) −1/2 ≥ 1. An estimate of the correlation coefficient ρ can be found by fitting the function in Fig. 6 by assuming that the number of stars with non-zero likelihood trace the underlying likelihood contours. By doing this, we obtain ρ = −0.74, which corresponds to a factor 1.5 for the uncertainties. This multiplication provides us with upper limits for the 1σ errors, as reported in Table 2. Notice that we consider this to be an overestimate of the real uncertainties because the individual contours are in reality nonnormal, non-linear and have slightly orthogonal constraints among each others.
The quoted precisions for α, q in our summary table are remarkable. This is a by-product of the extremely stringent condition that all HVS orbits should be radial and cross the ejection region near the GC, which represents a limited volume of the Galactic phase-space. In our numerical implementation, this volume is determined by the hyperparameters σ r and σ L , which set the maximum distance from the GC, r, and the maximum angular momentum, L, allowed inside the ejection region. In our testing, relaxing the condition on the angular momentum worsens the constraints in α, q considerably, meaning that the zero-angular momentum condition is the dominant factor that allows HVSs to constrain the NFW profile. Fig. 7 show the orbital characteristics of the strong, poor and average constrainers for the parameters α, q.

Observational prospects
From the scatter plots, it is clear that there is a correlation between how constraining stars are and how much time they have spent being affected by the gravitational potential (see flight time panel). The most powerful stars in our golden sample are therefore tightly bound and have spent hundreds or thousands of Myr orbiting around the Galaxy. Unfortunately part of these stars spend most of their time in a region where the Galactic Disc dominates the gravitational potential and while we have assumed perfect knowledge of this component, in reality this will hinder the halo reconstruction. In addition, the identification of these HVSs is difficult because of their low Galactocentric velocities.
On the other hand, we identify a useful sample made from average constrainers for α and strong constrainers for q, located at Galactocentric distances above 2 kpc. Because around half of this sample is moving along unbound trajectories, the identification of these HVSs is feasible.
The distributions in Galactocentric velocity and position presented here set clear targets for observations aimed The results shown are for model A for the Milky Way's dark matter halo, but identical trends are found in model B and C. Overall, these plots show the presence or absence of correlation between the kinematic properties of HVSs and their ability to constrain the parameter α (left column) or q (right column), parametrized by the expected individual relative errors σ α /α 0 and σ q /q 0 . In the ejection velocity plots, the two vertical dashed lines mark, from left to right, the minimum velocity necessary to reach a Galactocentric distance equal to the scale radius r s and 250 kpc respectively. at measuring the Galactic halo with HVSs. Notice in particular that in order to produce an average constraint for α, a HVS needs an ejection velocity sufficiently high to reach Galactocentric distances equal to the scale radius r s , but it is not required to be there when observed. This is not surprising, since α is the effective parameter measured in the M s − r s plane. Note also that while not shown, the results for halo B and C follow the same trends.
Another important factor to consider regarding future observations is the presence of uncertainties in the astrometry of our stellar sample. In Appendix A we show that for true HVSs, Gaia-like observations should not affect our halo reconstruction, but might still cause complications in the identification of HVSs. In particular, we expect a certain amount of contamination due to stars which, by chance, follow HVS-compatible orbits.
While quantifying the impact of this contamination and correcting for its effect in the inferred halo parameters is not the goal of this paper, we can still quantify its expected magnitude using simple arguments. Robin et al. (2012) estimated the number of halo stars in the final Gaia catalogue to be equal to 10 7 . Of these, around 10 4 will have velocities higher than our golden sample threshold of 450 km/s. Assuming an isotropic velocity distribution and the fiducial star with Gaia error bars (see Appendix A), the fraction of stars with a proper motion vector consistent with the radial direction is then ∼ 10 −3 . According to this estimate, the number of halo stars polluting our sample should be ∼ 10, close to the ∼ 100 real HVSs that we expect in our average constrainer class. Notice however that this bound is particularly conservative, since we have neglected additional information -like metallicity -which correlate greatly with being a HVS.

DISCUSSION AND CONCLUSIONS
Hypervelocity stars are remarkable objects. According to the leading model, they are ejected from the GC with high velocity (around 10 3 km/s) and travel along orbits spanning at least tens of kiloparsecs. This allows them to probe the gravitational potential of the Milky Way where the dark matter halo is dominant.
In this work we have developed a technique to extract information about the Galactic potential and the ejection mechanism from the observed HVS distribution in mass, velocity and position. Our method predicts the density of HVSs for a given stellar mass and phase-space position by back-propagating the observed location to the ejection point. The orbit is therefore required to cross the GC within a stellar lifetime to result in a non-zero distribution function. This is the basis of our likelihood pipeline, used to produce model constraints. To test our method we have applied it to mock HVS populations, designed to mimic what the European Space Agency's mission Gaia will observe in the next few years. In our simulations, HVSs are propagated in three fiducial axisymmetric potentials and then used to reconstruct the dark matter components, modelled using a spheroidal NFW potential defined by a scale radius, scale mass and axis ratio (r s , M s , q).
The results of our analysis are very promising, we find that ∼ 200 HVSs are able to provide an unbiased measurement of the NFW potential parameters with sub-percent uncertainties, thanks mainly to the strict constraints we impose on the ejection location and angular momentum at that instant. While promising, it should be kept in mind that our results were obtained in an idealized scenario. We assumed perfect knowledge of the baryonic potential and the parametric form of dark matter halo, not accounting for modelling errors.
We test the robustness of our results for both spherical and oblate geometries, and for two different values of the fiducial scale parameters. In all cases we also observe a natural degeneracy, whereby only the combination α = M s /r 2 s can be constrained. We also identify two special classes of stars, named "poor" and "strong" constrainers. The first class contains 15−30% of our sample and the stars in it are not be able to produce likelihood contours because, of the model explored, only the fiducial one produces a non-zero likelihood. The second class, similarly sized, contains stars which are unable to tell the majority of the potentials in our grid from each other. If we neglect the poor constrainers, we identify a useful sample of ∼ 150 stars which can individually measure α with a precision of ∼ 20% and q with a precision of ∼ 5%.
In Fig. 8 we summarize the final results of this paper by showing how many stars in our sample we expect to allow a given halo potential. It should be noted that the correlation visible in the rightmost panel of this figure is not a characteristic of individual HVSs, but arises because single HVSs constrain different combinations of q and α. To account for this in our estimates, in Sec. 4.3 we have estimated the covariance between the two variables in the stacked likelihood. Our final estimate of the 1σ relative uncertainties is < 1% for α and < 0.1% for q (see Table 2).
We also show that the constraining power of HVSs correlates with some observational quantities. In particular, we identify two essential properties characterizing a useful sample of HVSs: their orbits should be able to reach the NFW scale radius and their flight-time should be as long as possible. This roughly translates into distances from the GC between 2 and 20 kpc and Galactocentric velocities 900 km/s. We then discuss how our method can be applied to real data and examine the impact of observational challenges. We find that if the origin of a HVS is assumed to be the GC, Gaia-like uncertainties have no impact on the reconstructed radial orbit. Furthermore, the number of contaminants moving along quasi-radial orbits by chance should be negligible.
At last, while a detailed comparison between our forecast and actual measurements of the Milky Way halo using other probes is not straightforward, we find it useful to report the result for our primary model (halo A) in a standard format. Notice that even if we assume a spherical halo, the degeneracy in the M s −r s plane does not allow us to constrain the virial mass M 200 or the virial radius R 200 . 3 This degeneracy can be broken if we assume that the Milky Way concentration parameter c = R 200 /r s is related to the virial mass through a mass-concentration relation, as seen in ΛCDM numerical simulations (see, e.g. Navarro et al. 1997 Figure 8. Summary of the simulated constraints obtained in this paper using 197 HVSs propagated in a Galaxy with dark matter halo A. Assuming a spheroidal NFW potential, HVSs are able to measure the axis ratio q and the effective scale parameter α = M s /r 2 s . The plots on the left show number of stars allowing a certain value for one of the parameters (q or α) when freezing the second one to its fiducial expectation (q 0 for one and α 0 for the other). By looking at the likelihood evaluated along to the two directions marked by the dashed lines in the α − q plane, and taking into account covariance, we are able to estimate the marginalized 3σ error bars visible on the top-most and right-most side of the contour plot (see Sec. 4.1 and 4.2). An estimate for the correlation between q and α is found in Sec. 4.3. The plotted ellipses represent the 5 and 10 sigma contours for the inferred bivariate distribution. The figure shows the robustness of our result to changes in the fiducial values for α and q (halo A, B, C, see Table 1) and it illustrates how the width of the peaks in the histograms translates, through the combination of multiple stars, in tight constraints for the halo parameters.
assuming a spherical halo, we use the relation from Dutton & Macciò (2014) and the latest Planck 2015 cosmology Planck Collaboration (2016) to translate our precision in α into a virial mass of M 200 = (1.5 ± 0.2) × 10 12 M , corresponding to ∼ 10% precision. Notice that, although our reconstruction of α is not affected by bias, the recovered virial mass is 2.3σ away from the true M 200 = 1.04 × 10 12 M of the fiducial halo A. As observed before in Wang et al. (2015), this is a perfect example of how assumptions, like imposing a massconcentration relation, can affect the results obtained with dynamical tracers. A direct comparison of our forecast with the constraints of other probes provided by the same paper (their figure 1) also suggests that our technique is able to achieve competitive results.
Our conclusions paint an optimistic picture for the introduction of HVSs as a new dynamical tracer of the Galactic potential, especially when combined with the prospects of HVS detections in the final release of Gaia (Marchetti et al. 2018). The wealth of data that will become available in the next few years will allow measurements of the dark matter distribution in the Milky Way of unprecedented precision. However, in order to produce accurate results and combine the information provided by multiple tracers, particular care should be taken and modelling biases be carefully considered.
we focus on initial velocity and time because the likelihood in eq. (15) is sensitive to these parameters. With our test, we verify that a change of the order of 5 km/s and 0.1 Myr in the initial conditions is enough to obtain a proper motion vector inconsistent with the one obtained from the real v 0 , t 0 . These deviations should be compared to the typical scale functions of the individual likelihood function. The typical time-scale of the survival function g(t, m) is the lifetime t L (m), which respects ∆t 0 t L (m) for our ranges of stellar mass m, while the typical time-scale for the velocity is dictated by R H and is equal to 1000 km/s ∆v 0 . This confirms that the freedom of initial conditions allowed by the observational errors ∆w does not affect the inferred parameters of the gravitational potential, because all orbits within the bounds provide an identical constraint. This test is promising, but there is still the possibility that sets of initial parameters disjoint from the real w 0 , t 0 might result in an orbit crossing the configuration space volume delimited by the error bar ∆w . Said otherwise, it is possible that very different initial conditions might result in orbits which eventually become close enough to be confused with each other. We argue that the chances of this happening are slim because HVS trajectories, originated from a small region of phase space, are rare in the Galaxy. Notice, in particular, that the chances of this happening are manifestly zero for unbound stars, which travel along mostly radial orbits. In conclusion: for our purposes we expect only one HVS orbit to be within the Gaia-like error bars.
Unfortunately, exactly because HVSs' orbits are rare, one of the main complications of HVS searches is confirming if an observed position w ± ∆w is consistent with being a HVS (see, e.g., Brown et al. 2015;Marchetti et al. 2017). Note that previous attempts have neglected the almost zero angular momentum condition which allows HVSs to constrain the halo parameters with high precision. As noted in Hattori et al. (2018), this requirement is particularly restrictive: when transforming from Equatorial to Galactocentric coordinates, this condition can also be used to constrain the Sun's motion relative to the GC.
To identify the allowed phase space regions consistent with true HVS orbits a fine sampling within the error bars ∆w is required. As an example, consider a star with angular momentum at ejection equal to our smoothing value σ L = 10 pc × km/s and let us assume conservation of angular momentum. When the star reaches ∼ 10 kpc from the GC, the non-radial component of the velocity must then be known with ∼ 10 −3 km/s precision. Note that this value is at least three orders of magnitude smaller than the expected velocity error-bars mentioned at the start of this Section.
To expedite this process, we propose to start by artificially enlarging the phase space volume associated to the GC by modifying the pericentre and angular momentum smoothing lengths σ r , σ L . This has the intended effect of allowing a larger region within the observed w ± ∆w to be consistent with a trajectory crossing the GC. An iterative search can then be performed, where the parameters σ r , σ L are slowly lowered and the observed phase space error ∆w shrunk until a HVS orbit is found. In our case, for every step, we zoom in around the maximum likelihood orbit.
To test the feasibility of this approach we apply it to a few average constrainers in our sample. We begin with an observed position w which we shift by a factor 0.5∆w compared to the real one to simulate noisy data, and for each iteration step we back-propagate 10 3 positions within the estimated uncertainties in velocity and distance from the sun. Starting from Gaia-like uncertainties, we obtain in < 10 iterations a phase space region 10 −8 times smaller, which contains orbits passing through the GC (according to our original strict σ r , σ L ) and with the same initial condition as the starting w .
This procedure we suggest effectively fits the position of the observed star itself while fitting the potential parameters under the assumption that it is a HVS. As discussed in Sec. 4.4 this exposes us to contamination from halo stars which might be observationally consistent with radial orbits. This contamination is found however not to be dominant.