# The escape speed curve of the Galaxy obtained from Gaia DR2 implies a heavy Milky Way

G. Monari, Benoit Famaey, I. Carrillo, T. Piffl, M. Steinmetz, R. F. G. Wyse, F. Anders, C. Chiappini, K. Janssen

► To cite this version:

G. Monari, Benoit Famaey, I. Carrillo, T. Piffl, M. Steinmetz, et al.. The escape speed curve of the Galaxy obtained from Gaia DR2 implies a heavy Milky Way. Astronomy and Astrophysics - A&A, 2018, 616, pp.L9. 10.1051/0004-6361/201833748. hal-02410681

HAL Id: hal-02410681

https://hal.science/hal-02410681

Submitted on 3 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LETTER TO THE EDITOR

# The escape speed curve of the Galaxy obtained from Gaia DR2 implies a heavy Milky Way

G. Monari<fnref n="1" />, B. Famaey<fnref n="2" />, I. Carrillo<fnref n="1" />, T. Piffl<fnref n="1" />, M. Steinmetz<fnref n="1" />, R. F. G. Wyse<fnref n="3" />,<fnref n="4" />, F. Anders<fnref n="1" />,
C. Chiappini<fnref n="1" />, and K. Janßen<fnref n="1" />

Received 29 June 2018 / Accepted 5 July 2018

ABSTRACT

We measure the escape speed curve of the Milky Way based on the analysis of the velocity distribution of ~2850 counter-rotating halo stars from the Gaia Data Release 2. The distances were estimated through the StarHorse code, and only stars with distance errors smaller than 10% were used in the study. The escape speed curve is measured at Galactocentric radii ranging from ~5 kpc to ~10.5 kpc. The local Galactic escape at the Sun's position is estimated to be $v_c(r_\odot) = 580 \pm 63$ km s$^{-1}$, and it rises towards the Galactic centre. Defined as the minimum speed required to reach three virial radii, our estimate of the escape speed as a function of radius implies for a Navarro–Frenk–White profile and local circular velocity of 240 km s$^{-1}$ a dark matter mass $M_{200} = 1.28_{-0.50}^{+0.68} \times 10^{12} M_\odot$ and a high concentration $c_{200} = 11.09_{-1.79}^{+2.94}$. Assuming the mass-concentration relation of ΛCDM, we obtain $M_{200} = 1.55_{-0.51}^{+0.64} \times 10^{12} M_\odot$ and $c_{200} = 7.93_{-0.27}^{+0.33}$ for a local circular velocity of 228 km s$^{-1}$.

**Key words.** Galaxy: kinematics and dynamics – Galaxy: fundamental parameters

1. Introduction

While much progress has been made in recent years regard-
ing 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 valid-
ity 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 exhaus-
tive model, we can already take advantage of the accuracy, qual-
ity, 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. Sev-
eral 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¹<fnref n="2" /> 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 sim-
ilar to those of the Galactic halo are bound to the Galaxy and

showed that the virial mass should be higher than 1.4 × 10¹<fnref n="2" /> 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
1.67±0.79 × 10¹<fnref n="2" /> M☉. Finally, Posti & Helmi (2018) also used the
kinematics of globular clusters to fit a two-component distribu-
tion 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 dif-
ference 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 previ-
ous 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

<footnote n="1">¹ Leibniz – Institut fuer Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
e-mail: gmonari@aip.de</footnote>

<footnote n="2">² Université de Strasbourg, Observatoire astronomique de Strasbourg, CNRS UMR 7550, 11 rue de l'Université, 67000 Strasbourg, France</footnote>

<footnote n="3">³ Department of Physics & Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA</footnote>

<footnote n="4">⁴ Leverhulme Trust Visiting Professor, University of Edinburgh, Edinburgh, UK</footnote>
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;

- VISIBILITY_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 AllWISE) 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 r = √R² + z² and the speed v = √v_R² + v_φ² + v_z². 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),

$$f(v|v_e, k) = \begin{cases} (k+1)(v_e - v)^k / (v_e - v_{\text{cut}}), & v \le v_e, \\ 0, & v > v_e, \end{cases} \quad (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

$$P(v_e, k|v_{i=1,...,N}) = \frac{P(v_e)P(k)\Pi_{i=1}^N f(v_i|v_e, k)}{\int \int P(v_e)P(k)\Pi_{i=1}^N f(v_i|v_e, k)dv_e dk}, \quad (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,

$$P(v_e|v_{i=1,...,N}) = \int_{2.3}^{3.7} P(v_e, k|v_{i=1,...,N})dk, \quad (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 v_e^i and its uncertainty σ_e^i = max(v_e^i - v_{16}^i, v_{84}^i - v_e^i), where v_{16}^i and v_{84}^i 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 v_e^i estimates of each iteration, and its uncertainty as the sum in quadrature of the square root of the mean σ_e^i² and of the standard deviation of the values v_e^i. 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).

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
<figure><img src="image_2.png" /><figcaption>Fig. 1. Histograms in different r-bins for one realisation of r and v of the sample, compared with $f(v|v_e, k)$ (blue line), where $v_e$ is derived for that particular bin and $k=3$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

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

$$v_e(r_{\odot}) = 580 \pm 63 \text{ km s}^{-1}. \qquad (4)$$

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_{\odot}) = 593 \pm 76 \text{ km s}^{-1}$. 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 $\Phi(x)$ is the potential of the

Galaxy at some position x, and $\Phi \to 0$ at infinity, then

$$v_{e,\infty}(x) = \sqrt{2|\Phi(x)|}. \qquad (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_{340}$, where $r_{340}$ is the spherical radius within which the average density of the whole Galaxy is equal to 340 times the critical density at redshift 0,

$$\rho_c = 3H^2/8\pi G, \qquad (6)$$

where we take $H = 73 \text{ km s}^{-1} \text{ Mpc}^{-1}$. The distance $3r_{340}$ 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<fnref n="1" />, so that it depends on r alone, that is, $\Phi(r)$. The escape speed can in this case be modelled as $v_{e,\infty}(r)$ or as

$$v_{e,340}(r) = \sqrt{2|\Phi(r) - \Phi(3r_{340})|}, \qquad (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_{200}$ and its concentration $c_{200}$, both estimated at $r_{200}$, the radius where the average density of the dark halo reaches 200 times the critical density $\rho_c$. As an additional constraint, we imposed that the circular velocity at the Sun is $v_c(r_{\odot}) = 240 \text{ km s}^{-1}$, as in the assumptions that we used to transform l.o.s. velocities and proper motions to v. We obtain a fit using a $\chi^2$ minimisation, with the points weighted according to their error bars. The best-fit model corresponds to the orange line and orange $1\sigma$ error bands in the left panel of Fig. 3. The error bands are obtained using the $\chi^2$ statistics. The best-fit model has $M_{200} = 0.97_{-0.36}^{+0.47} \times 10^{12} M_{\odot}$ and $c_{200} = 12.63_{-2.09}^{+3.41}$ if the fit is made assuming the $v_{e,\infty}$ interpretation of the escape speed, and $M_{200} = 1.28_{-0.50}^{+0.68} \times 10^{12} M_{\odot}$ and $c_{200} = 11.09_{-1.79}^{+2.94}$ if the fit is made with the more realistic $v_{e,340}$ (shown in the figure).

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

$$\log_{10}(c_{200}) = 0.905 - 0.101 \log_{10}(M_{200}/[10^{12} \text{ h}^{-1} M_{\odot}]), \qquad (8)$$

<footnote n="1">¹ The escape speed estimated in the Galactic plane or towards the Galactic pole at a fixed r can differ by a few tens of km s⁻¹, which is included in the error bars of our estimate of the escape speed.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Left panel: fit to the escape velocity points assuming $v_{e,340}$, and obtained using the disc and bulge potential of model III of Irgang et al. (2013), and a dark halo of $M_{200} = 1.28 \times 10^{12} M_{\odot}$ and $c_{200} = 11.09$ (orange line and uncertainty bands). Right panel: same as in the left panel, but for $M_{200} = 1.55 \times 10^{12} M_{\odot}$ and $c_{200} = 7.93$ (blue line and uncertainty bands) when the ΛCDM relation between $M_{200}$ and $c_{200}$ is assumed. The orange and blue bands represent the $1\sigma$ uncertainties of the models.</figcaption></figure>

where $h = H/(100 \text{ km s}^{-1} \text{ Mpc}^{-1})$. We repeated the fit using this relation. In this case, the relation between $M_{200}$ and $c_{200}$ no longer guarantees that $v_c = 240 \text{ km s}^{-1}$ at $R_{\odot}$. 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_{\odot}$ (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\sigma$ error bands. The best-fit model in this case is $M_{200} = 1.16^{+0.47}_{-0.38} \times 10^{12} M_{\odot}$ and $c_{200} = 8.17^{+0.33}_{-0.28}$ in the $v_{e,\infty}$ interpretation, and $M_{200} = 1.55^{+0.64}_{-0.51} \times 10^{12} M_{\odot}$, $c_{200} = 7.93^{+0.33}_{-0.27}$ 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_{\odot}) = 228 \text{ km s}^{-1}$. Such a high mass for the Milky Way dark matter halo is expected from abundance matching. Behroozi et al. (2013) estimated that $\log(M_{200}) \approx 12.25$ should correspond to a stellar mass of between $3 \times 10^{10} M_{\odot}$ and $5.5 \times 10^{10} M_{\odot}$.

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

$$ \Phi_h = \frac{v_0^2}{2} \log(r^2 + r_c^2). \qquad (9) $$

The fit was obtained using $v_{e,340}$ ($\Phi_h$ diverges at infinity) and imposing that $v_c(r_{\odot}) = 240 \text{ km s}^{-1}$. The best fit corresponds to $v_0 = 176^{+22}_{-22} \text{ km s}^{-1}$ and $r_c = 4.7^{+2.1}_{-3.8} \text{ kpc}$ (the error bars are not shown in Fig. 4). This fit corresponds to a high halo mass $M_{200} = 1.74^{+0.72}_{-0.58} \times 10^{12} M_{\odot}$, as expected for a pseudo-isothermal sphere with a density decreasing slowly as $r^{-2}$ 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_0 = 235 \text{ km s}^{-1}$ and $r_c = 10.7 \text{ kpc}$, which gives a too high a value for the escape speed if no cut-off on the $r^{-2}$ 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.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Same as in Fig. 2, but for the isothermal best fit model, with $v_0 = 176 \text{ km s}^{-1}$ and $r_c = 4.7 \text{ kpc}$ (blue line), compared with the model of Englmaier & Gerhard with $v_0 = 235 \text{ km s}^{-1}$ and $r_c = 10.7 \text{ kpc}$ (orange line).</figcaption></figure>

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,\infty}$ 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_{\odot}) = 580 \pm 63 \text{ km s}^{-1}$ (in $1\sigma$ agreement with Williams et al. 2017). The escape speed varies from ~650 km s$^{-1}$ at ~5 kpc to ~550 km s$^{-1}$ at 10.5 kpc. The uncertainty 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_{340}$ (as in P14). In the latter way, and fixing the circular velocity of the Galaxy at 240 km s⁻¹ at the Sun, $M_{200} = 1.28_{-0.50}^{+0.68} \times 10^{12} M_⊙$, a concentration $c_{200} = 7.93_{-0.27}^{+0.33}$ (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_{340}$, and we find $M_{200} = 1.74 \times 10^{12} 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).

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 ($M_{200} = 1.55_{-0.51}^{+0.64} \times 10^{12} M_⊙$), a concentration $c_{200} = 7.93_{-0.27}^{+0.33}$ (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_{340}$, and we find $M_{200} = 1.74 \times 10^{12} 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 ~2.64$_{-0.38}^{+0.42} \times 10^{12} M_⊙$. 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_{200} > 10^{12} M_⊙$, especially when considering the mass-concentration 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.

*Acknowledgements.* 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

Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17

Battaglia, G., Helmi, A., Morrison, H., et al. 2005, MNRAS, 364, 433

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bovy, J., Bird, J. C., García Pérez, A. E., et al. 2015, ApJ, 800, 83

Casagrande, L., & VandenBerg, D. A. 2018, MNRAS, 479, L102

Cole, D. R., & Binney, J. 2017, MNRAS, 465, 798

Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359

Englmaier, P., & Gerhard, O. 2006, Celestial Mech. Dyn. Astron., 94, 369

Famaey, B., Bruneton, J.-P., & Zhao, H. 2007, MNRAS, 377, L79

Fritz, T. K., Carrera, R., & Battaglia, G. 2018, A&A, submitted

[arXiv:1805.07350]

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Hattori, K., Valluri, M., Bell, E. F., & Roederer, I. U. 2018, ApJ, submitted

[arXiv:1805.03194]

Hawkins, K., & Wyse, R. F. G. 2018, ArXiv e-prints [arXiv:1806.07907]

Irrgang, A., Wilcox, B., Tucker, E., & Schiefelbein, L. 2013, A&A, 549, A137

Kafle, P. R., Sharma, S., Lewis, G. F., Robotham, A. S. G., & Driver, S. P. 2018, MNRAS, 475, 4043

Kochanek, C. S. 1996, ApJ, 457, 228

Kruijssen, J. M. D., Pfeffer, J. L., Reina-Campos, M., Crain, R. A., & Bastian, N. 2018, MNRAS, in press [arXiv:1806.05680]

Leonard, P. J. T., & Tremaine, S. 1990, ApJ, 353, 486

Lindegren, L., Hernandez, J., Bombrun, A., et al. 2018, A&A, 616, A2

Marchetti, T., Rossi, E. M., & Brown, A. G. A. 2018, MNRAS, submitted [arXiv:1804.10607]

McMillan, P. J. 2018, Res. Notes AAS, 2, 51

Peñarrubia, J., Gómez, F. A., Besla, G., Erkal, D., & Ma, Y.-Z. 2016, MNRAS, 456, L54

Piffl, T., Scannapieco, C., Binney, J., et al. 2014, A&A, 562, A91

Posti, L., & Helmi, A. 2018, A&A, submitted [arXiv:1805.01408]

Queiroz, A. B. A., Anders, F., Santiago, B. X., et al. 2018, MNRAS, 476, 2556

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130

Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829

Smith, M. C., Ruchti, G. R., Helmi, A., et al. 2007, MNRAS, 379, 755

Watkins, L. L., van der Marel, R. P., Sohn, S. T., & Evans, N. W. 2018, AAS J., submitted [arXiv:1804.11348]

Williams, A. A., Belokurov, V., Casey, A. R., & Evans, N. W. 2017, MNRAS, 468, 2359