© ESO 2018

1. Introduction

While much progress has been made in recent years regarding our understanding of the formation of galaxies, the very basic fundamental parameters of our own Galaxy are still very poorly known. Of these fundamental parameters, the mass of our Galaxy is perhaps the most important, and it is still unknown within a factor of four (Bland-Hawthorn & Gerhard 2016). This parameter is of particular interest because it provides a test for the current Λ cold dark matter (ΛCDM) paradigm, and the validity of recipes such as abundance matching (e.g. Behroozi et al. 2013).

The primary objective of the cornerstone Gaia mission (Gaia Collaboration 2018) is to allow us to build a detailed dynamical model of the Galaxy, including all of its components, and to provide us insight into its structure, its formation, and its evolutionary history. Before building such a detailed and exhaustive model, we can already take advantage of the accuracy, quality, and extent of the data provided by the Gaia Data Release 2 (DR2) to determine the fundamental parameters of the Galaxy such as its virial mass more accurately than ever before. Several recent studies have tried to estimate this parameter through various methods: Fritz et al. (2018) found that a virial mass of at least 1.6×10¹²M_⊙ would be needed to keep a reasonable number of satellite dwarf galaxies bound. Hattori et al. (2018) assumed that extremely high velocity stars with chemical properties similar to those of the Galactic halo are bound to the Galaxy and showed that the virial mass should be higher than 1.4 × 10¹²M_⊙. Hawkins & Wyse (2018) also studied a sample of hyper-velocity stars and found that their chemistry is compatible with that of the stellar halo. Watkins et al. (2018) used a mass estimator based on the kinematics of halo globular clusters to find a virial mass of Finally, Posti & Helmi (2018) also used the kinematics of globular clusters to fit a two-component distribution function in action space as well as the Galactic potential including a prolate Navarro–Frenk–White (NFW) halo, to find a slightly lower value for the virial mass.

Here, our new estimate of the Milky Way mass will be made through measuring the escape speed curve of the Galaxy at Galactocentric radii ranging from ~5 kpc to ~10.5 kpc based on Gaia DR2. Since the escape speed directly measures the difference between the local potential and that in the outskirts of the Galaxy, it allows us to constrain its total mass based on local measurements. In doing this, we follow up on the previous studies by Leonard & Tremaine (1990), Kochanek (1996), Smith et al. (2007), and Piffl et al.(2014; hereafter P14).

The general idea is to use a simple model for the tail of the velocity distribution of halo stars at a given radius, based on a truncated power law (Leonard & Tremaine 1990) motivated by simulations (Smith et al. 2007; P14). Based on a model for the baryonic distribution of the Milky Way, the virial mass of the Galaxy can be adjusted to best fit the escape speed curve. In Sect. 2 we present our selection of Gaia DR2 stars. The method for estimating the escape speed is then presented in Sect. 3, and the results are summarized in Sect. 4. The corresponding virial mass of the Galaxy dark matter halo is presented in Sect. 5, and we conclude in Sect. 6.

2. Data

We used stars from the Gaia DR2 with line-of-sight (l.o.s.) velocity information. Following Marchetti et al. (2018), we selected only stars with the following Gaia flags:

• ASTROMETRIC_GOF_AL < 3;
  

• ASTROMETRIC_EXCESS_NOISE_SIG ≤ 2;
  

• −0.23 ≤ MEAN_VARPI_FACTOR_AL ≤ 0.32;
  

• VISIILITY_PERIODS_USED > 8;
  

• – RV_NB_TRANSITS > 5.

We used new distance estimates (Anders et al. in prep.) obtained with the StarHorse code (Queiroz et al. 2018), which combines multiband photometric information (APASS, 2MASS, and All-WISE) and the Gaia astrometric information with a Bayesian approach, accounting for the global Gaia DR2 parallax zero-point shift of −0.029 mas (Lindegren et al. 2018; Arenou et al. 2018). Based on these distances and on the Gaia DR2 astrometry and l.o.s. velocities, we computed the Galactocentric cylindrical positions (R, z, ϕ) and velocities (v_R, v_ϕ, v_z) of these stars. For this, we assumed a distance of the Sun from the centre of the Galaxy of r_⊙ = 8.34 kpc, a circular velocity at the Sun of v_c(r_⊙) = 240 km s⁻¹ (Reid et al. 2014), and a peculiar solar motion of (U_⊙, V_⊙, W_⊙) = (11.1, 12.24, 7.25) km s⁻¹, as estimated by Schönrich et al. (2010). Of these stars we selected only those that counter-rotate, that is, stars with v_ϕ < 0, in order to ensure that the kinematic tail is representative of stellar halo-type stars. We checked that reflecting v_ϕ for these stars around v_ϕ = 0, we obtain the kind of kinematics typical of the stellar halo, with velocity dispersions of ~140 km s⁻¹ (e.g. Battaglia et al. 2005).

Of these stars, we selected only those with distance errors smaller than 10 percent. We then computed their Galactocentric spherical distance and the speed . However, we did not simply combine the observable quantities, but used a Monte Carlo technique (explained below), which allowed us to achieve better estimates of the r and v uncertainties than a simple propagation of errors.

3. Methods

As in P14, we assumed that the velocities v of stars in the kinematic tail of the halo are distributed according to a power law, with a truncation at the escape velocity v_e (Leonard & Tremaine 1990),(1)

where v_cut corresponds to a lower cut in v that defines the fast stars that should be distributed as in Eq. (1). We note that v_e and k depend in general on the position in the Galaxy. Both v_cut and typical values of k can be read off simulations. In particular, P14 have shown from simulations that the typical range for k is 2.3 < k < 3.7 and that stars follow this power law for v_cut ~ 250 km s⁻¹.

The Bayes theorem can be used to show that the probability of the model parameters (v_e, k), given N stars with velocity determination v_i (for the i-th star) is(2)

where the uncertainty on the determination of v_i is neglected for the moment, and the denominator is the same for all the points in the (v_e, k) space, so that we do not need to determine it. For the a priori probabilities we chose P(v_e) ∝ 1/v_e and P(k) uniform, as in Leonard & Tremaine (1990).

Since we are interested mostly in v_e, we finally marginalized along k to obtain the probability density of v_e,(3)

where the range of marginalisation is obtained from P14, as explained above.

However, at this stage, Eq. (2) does not consider the uncertainty on the individual stellar parameters (r and v).

Hence, in order to obtain a robust result, we used a Monte Carlo or bootstrap technique, repeating the scheme detailed below 100 times.

• For each star in the sample, we drew a random distance, proper motion, and l.o.s. velocity from Gaussians with means and dispersions (the quantity and its uncertainty) given by the StarHorse code and the Gaia catalogue; from these, we computed r and v assuming for the transformation the same Galactic parameters as were used in Sect. 2, and negligible uncertainties on the positions on the sky of the stars.

• We used only stars with drawn random distance smaller than 5 kpc in order to have a sample of stars with Bayesian distance estimates dominated by the Gaia parallax information, and not by the photometry; considering the parallax shift, 99% of the stars in Gaia DR2 with l.o.s. velocities and parallax errors better than 10% are found inside ~5 kpc, and inside the same sphere, 90% of the stars have parallax error smaller than ~10%; the typical size of the samples obtained in this way is ~2850 stars.

• We resampled the data set obtained in this way with replacement, so that each time some stars are randomly excluded and some stars are repeated more times in the sample.

• We binned the data into bins of size 0.6 kpc, with the innermost bin centred at r = 5.34 kpc and the outermost at r = 10.14 kpc.

• For the i-th iteration (and resampling), we calculated for each bin the escape velocity and its uncertainty , where and are the 16th and 84th percentiles of the v_e distribution.

We then estimated v_e for each r-bin as the median of the estimates of each iteration, and its uncertainty as the sum in quadrature of the square root of the mean and of the standard deviation of the values . The error bars obtained in this way contain a systematic component because the top of the error bar for each bin always depends critically on the fastest star that can be included in the sample in one of the 100 Monte Carlo realisations.

4. Escape speed curve

In Fig. 1 we show the histograms in different r-bins for one realisation of r and v of the sample in order to check that a power law is a reasonable description of the velocity distribution. We overplot f (v|v_e, k), where v_e is derived for those particular bins, and k = 3 (the mean between k = 2.3 and k = 3.7, used for the marginalisation to obtain the v_e probability distribution function).

Fig. 1. Histograms in different r-bins for one realisation of r and v of the sample, compared with f (v|ve, k) (blue line), where ve is derived for that particular bin and k = 3.

In the left panel of Fig. 2, we show the results of the analysis in different bins (black dots with error bars). Only the bins that contain at least ten stars in each resampling were considered. The escape velocity curve decreases with r, as expected in any reasonable Galactic potential (see below). The value of the escape velocity at the Sun is(4)

Fig. 2. Escape speed at different radii r in the Milky Way obtained with the StarHorse code (black points), and McMillan distances (red points) and their error bars obtained using the procedure explained in Sect. 3.

The error bar increases for the farthest bins, especially the outermost bins that are less densely populated.

We then compared the results obtained with the StarHorse code with the distance estimates by McMillan (2018). We used the same counter-rotating stars as selected in Sect. 2, but with the velocities v and radii r, computed with the new distances. The results in different bins are also shown in Fig. 2 with red dots and associated error bars. The results for all bins are very similar to the results obtained with the StarHorse distances. The escape velocity at the Sun is estimated in this case to be v_e(r_⊙) = 593 ± 76 km s⁻¹. However, the farthest bins do not contain enough stars to perform the analysis. With the McMillan distances it is not possible to go very far from the Sun because the distance quality degrades rapidly and the bins in our sample rapidly become empty, while the multiband photometric information used in StarHorse allows us to achieve accurate distances deeper inside the Galaxy.

5. Dark matter halo mass and concentration

The escape velocity is in principle defined as the velocity necessary to bring a star to infinity. If Φ(x) is the potential of the Galaxy at some position x, and Φ → 0 at infinity, then(5)

However, it is unphysical to think that stars faster than v_e derived in the previous section can reach infinity; they will instead turn around beyond some very large distance. On the basis of simulations, P14 chose this distance to be 3r₃₄₀, where r₃₄₀ is the spherical radius within which the average density of the whole Galaxy is equal to 340 times the critical density at redshift 0,(6)

where we take H = 73 km s⁻¹ Mpc⁻¹. The distance 3r₃₄₀ approximates the virial radius of the Galaxy.

To model the escape speed, we used the Milky Way potential model III of Irrgang et al. (2013), which fits a number of observables of the Milky Way. As our estimates of the escape speed depend only on the spherical radius, we computed the Irrgang potential only within the Galactic plane¹, so that it depends on r alone, that is, Φ(r). The escape speed can in this case be modelled as v_e,∞(r) or as(7)

the latter being more meaningful in a cosmological context. We kept the Irrgang model III disc and bulge potentials fixed and varied the parameters of the NFW halo, which we chose to describe with its “virial” mass M₂₀₀ and its concentration c₂₀₀, both estimated at r₂₀₀, the radius where the average density of the dark halo reaches 200 times the critical density ρ_c. As an additional constraint, we imposed that the circular velocity at the Sun is v_c(r_⊙) = 240 km s⁻¹, as in the assumptions that we used to transform l.o.s. velocities and proper motions to v. We obtain a fit using a χ² minimisation, with the points weighted according to their error bars. The best-fit model corresponds to the orange line and orange 1σ error bands in the left panel of Fig. 3. The error bands are obtained using the χ² statistics. The best-fit model has and if the fit is made assuming the v_e,∞ interpretation of the escape speed, and and if the fit is made with the more realistic v_e,340 (shown in the figure).

Fig. 3. Left panel: fit to the escape velocity points assuming ve,340, and obtained using the disc and bulge potential of model III of Irrgang et al. (2013), and a dark halo of M200 = 1.28 × 1012M⊙ and c200 = 11.09 (orange line and uncertainty bands). Right panel: same as in the left panel, but for M200 = 1.55 × 1012M⊙ and c200 = 7.93 (blue line and uncertainty bands) when the ΛCDM relation between M200 and c200 is assumed. The orange and blue bands represent the 1σ uncertainties of the models.

However, N-body simulations in ΛCDM cosmology display a relation between the halo concentrations and their mass (Dutton & Macciò 2014),(8)

where h = H/(100 km s⁻¹ Mpc⁻¹). We repeated the fit using this relation. In this case, the relation between M₂₀₀ and c₂₀₀ no longer guarantees that v_c = 240 km s⁻¹ at R_⊙. We note that the value of the circular velocity at the solar circle is degenerate with the actual peculiar velocity of the Sun, so that a lower value of the circular velocity and a higher value of V_⊙ (e.g. Bovy et al. 2015) are still compatible with our assumptions. The fit is shown in the right panel of Fig. 3 with the blue line and blue 1σ error bands. The best-fit model in this case is and in the v_e,∞ interpretation, and in the more realistic v_e,340 interpretation. In the former case, the circular speed is still well behaved in the Galaxy and corresponds to v_c(r_⊙) = 228 km s⁻¹. Such a high mass for the Milky Way dark matter halo is expected from abundance matching. Behroozi et al. (2013) estimated that log(M₂₀₀) ≃ 12.25 should correspond to a stellar mass of between 3 × 10¹⁰M_⊙ and 5.5 × 10¹⁰M_⊙.

Finally, in Fig. 4 we show with a blue line the fit to the escape velocity points that we obtained with a quasi-isothermal potential for the dark halo of the form(9)

Fig. 4. Same as in Fig. 2, but for the isothermal best fit model, with v0 = 176 km s−1 and rc = 4.7kpc (blue line), compared with the model of Englmaier & Gerhard with v0 = 235 km s−1 and rc = 10.7 kpc (orange line).

The fit was obtained using v_e,340 (Φ_h diverges at infinity) and imposing that v_c(r_⊙) = 240 km s⁻¹. The best fit corresponds to and (the error bars are not shown in Fig. 4). This fit corresponds to a high halo mass , as expected for a pseudo-isothermal sphere with a density decreasing slowly as r⁻² in the outer parts. We also compared our escape speed curve with the v_e,340 values obtained from the halo model fitting the Milky Way gas dynamics by Englmaier & Gerhard (2006), with v₀ = 235 km s⁻¹ and r_c = 10.7 kpc, which gives a too high a value for the escape speed if no cut-off on the r⁻² is applied in the outer halo.

The analysis we performed in Sect. 4 can also be performed by keeping the exponent k fixed. In this case, a larger v_e is obtained for larger k because of the correlation between these two parameters, as has been shown by Leonard & Tremaine (1990). To check that our estimates are consistent with such an alternative method, we derived v_e for a fixed value k = 2.3 and k = 3.7, the extremes of the k-range used in our marginalisation. The v_e points obtained in this way are either shifted at lower or higher values, for k = 2.3 and 3.7 respectively, leaving the global shape of the curve unchanged. Fitting v_e,∞ and v_e,340 to these points, we then find models that are safely included within the error bands on mass and concentration given above.

6. Discussion and conclusions

We have calculated the escape speed curve of the Milky Way in a wide range of Galactocentric spherical radii using stars from the second data release of Gaia with line-of-sight velocities. The distances of these stars were computed using the StarHorse pipeline, which allows determining accurate and precise distances even for stars very far from the Sun thanks to its treatment of the extinction, especially in the central regions of the Galaxy. We used stars with distance estimates better than 10% at a distance of less than 6 kpc that are counter-rotating in order to have a stellar halo sample that is not contaminated by the disc.

In the solar neighbourhood we estimated the escape speed v_e(r_⊙) = 580 ± 63 km s⁻¹ (in 1σ agreement with Williams et al. 2017). The escape speed varies from ~650 km s⁻¹ at ~5 kpc to ~550 km s⁻¹ at 10.5 kpc. The uncertainly in the determination of v_e(r_⊙) is quite large despite the relatively large stellar sample, mainly because this uncertainty takes into account the errors in position and velocity of the stars in the analysis, and treats the possibility that they might be bound or not in the Galactic potential with a bootstrap technique without excluding stars on arbitrary criteria. Moreover, our determination of the escape speed at the Sun and at other radii in the Galaxy is “local” (i.e. performed in small radial bins), meaning that the estimates do not require any velocity rescaling that assumes a Galactic potential a priori (see P14); this introduces a bias in the result. Our determination of v_e(r_⊙) is slightly larger than but agrees well with what P14 determined based on stars from the RAVE survey, within the error bars. Our result also agrees with the finding by Hattori et al. (2018), who reported that the escape speed in the solar neighbourhood is expected to be of the order of ~600 km s⁻¹.

We have modelled the escape speed data points obtained using simple forms for the Galactic potential, including a fixed disc and bulge, by fitting a number of observables and a dark halo specified by the free parameters for the fit. The escape speed can be interpreted either as the velocity necessary to bring a star to infinity from a location in the Galaxy, or, more realistically, to a typical radius very far from the Galactic centre, motivated by cosmological simulations. The latter is chosen to be 3r₃₄₀ (as in P14). In the latter way, and fixing the circular velocity of the Galaxy at 240 km s⁻¹ at the Sun, .

With this fit imposing the local circular velocity at 240 km s⁻¹, we find a slightly more concentrated halo than in ΛCDM simulations. However, when we performed the fit imposing the ΛCDM relation between halo masses and concentrations, we find a higher halo mass (), a concentration (e.g. Kruijssen et al. 2018), and still a reasonable circular velocity curve. We also fit models consisting of the same disc and bulge potentials, but including cored pseudo-isothermal dark halos. In this case, the only interpretation for the escape velocity is the velocity to bring a star at 3r₃₄₀, and we find M₂₀₀ = 1.74 × 10¹²M_⊙, with a large core radius r_c = 4.7 kpc. In the future, we aim to determine whether such a measurement is compatible with alternative models (e.g. Famaey et al. 2007)

The high estimated mass of the Milky Way can be compared with recent estimates of the total mass of the Local Group based on the timing argument (Peñarrubia et al. 2016), which sets the total mass to . The heavy Milky Way implied by our work seems in accordance with recent estimates of the mass of M 31, yielding a lower mass than previously thought (0.8 ± 0.1 × 10¹²M_⊙, Kafle et al. 2018), which leaves room for an increased Milky Way mass.

A caveat of this work is that the distance determination for stars using Gaia stellar parameters is still an ongoing process, and the distances and velocities of the stars could undergo modifications, especially at large distances (e.g. because of the magnitude-dependent offset in the Gaia G magnitudes, see Casagrande & VandenBerg 2018). A dynamical caveat is that the analysis with which we estimated the escape velocities is based on distribution functions and assumes phase-mixing for the kinematic tail of the stellar halo, which is not guaranteed. However, the choice of the velocity distribution is motivated by simulations, and the comparison of the derived distribution functions with the histograms of the velocity distribution of the halo stars that we performed shows that this description is reasonable. Finally, the models of the Milky Way potential that we presented here are rather simple and do not take into account a number of issues, such as the fact that the dark halo should self-consistently adapt its shape to the bulge and disc potential (e.g. Cole & Binney 2017) and that the fit could be performed leaving the disc and bulge parameters free. Such detailed modelling is beyond the scope of this paper, however. Finally, spectroscopic follow-up of the stars we used in this analysis might be used to study whether their properties are compatible with the chemical properties of the stellar halo (e.g. Hawkins & Wyse 2018).

In conclusion, the uncertainties in the measurements of the escape velocity (and of the corresponding dark halo masses) are quite large, mostly because we were realistic and considered the uncertainties on both distance and velocity of the stars, and because we lack information on whether high-velocity stars are bound or unbound in the Galactic potential. Nevertheless, our estimate of the escape speed curve for a wide range of Galactocentric radii in the Milky Way conclusively implies a massive Galaxy M₂₀₀ > 10¹²M_⊙, especially when considering the massconcentration relation of ΛCDM, in agreement with other recent estimates using independent methods, both dynamical and using chemical tagging (e.g. Hawkins & Wyse 2018). It also agrees with typical abundance-matching expectations.

Acknowledgments

We thank Paul McMillan for useful discussions. RFGW thanks her sister, Katherine Barber, for support, and the Leverhulme Trust for a Visiting Professorship at the University of Edinburgh, held while this work was being completed. We thank the E-Science and Supercomputing Group at Leibniz Institute for Astrophysics Potsdam (AIP) for their support with running the StarHorse code on AIP cluster resources. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

All Figures

	Fig. 1. Histograms in different r-bins for one realisation of r and v of the sample, compared with f (v|ve, k) (blue line), where ve is derived for that particular bin and k = 3.
In the text

	Fig. 2. Escape speed at different radii r in the Milky Way obtained with the StarHorse code (black points), and McMillan distances (red points) and their error bars obtained using the procedure explained in Sect. 3.
In the text

Fig. 3. Left panel: fit to the escape velocity points assuming ve,340, and obtained using the disc and bulge potential of model III of Irrgang et al. (2013), and a dark halo of M200 = 1.28 × 1012M⊙ and c200 = 11.09 (orange line and uncertainty bands). Right panel: same as in the left panel, but for M200 = 1.55 × 1012M⊙ and c200 = 7.93 (blue line and uncertainty bands) when the ΛCDM relation between M200 and c200 is assumed. The orange and blue bands represent the 1σ uncertainties of the models.

In the text

	Fig. 4. Same as in Fig. 2, but for the isothermal best fit model, with v0 = 176 km s−1 and rc = 4.7kpc (blue line), compared with the model of Englmaier & Gerhard with v0 = 235 km s−1 and rc = 10.7 kpc (orange line).
In the text