© ESO, 2014

1. Introduction

In recent years quite a large number of studies concerning the mass of our Galaxy have been published. This parameter is of particular interest, because it provides a test of the current cold dark matter paradigm. There is now convincing evidence (e.g., Smith et al. 2007) that the Milky Way (MW) exhibits a similar discrepancy between luminous and dynamical mass estimates as already found in the 1970’s for other galaxies. A robust measurement of this parameter is needed to place the Milky Way in a cosmological framework. Furthermore, detailed knowledge of the mass and the mass profile of the Galaxy is crucial for understanding and modeling the dynamic evolution of the MW satellite galaxies (e.g., Kallivayalil et al. 2013, for the Magellanic clouds) and the Local Group (van der Marel et al. 2012b,a).

Generally, it can be observed that mass estimates based on stellar kinematics yield low values ≲10¹² M_⊙ (Smith et al. 2007; Xue et al. 2008; Kafle et al. 2012; Deason et al. 2012; Bovy et al. 2012), while methods exploiting the kinematics of satellite galaxies or statistics of large cosmological dark matter simulations find higher values (Wilkinson & Evans 1999; Li & White 2008; Boylan-Kolchin et al. 2011; Busha et al. 2011; Boylan-Kolchin et al. 2013). There are some exceptions, however. For example, Przybilla et al. (2010) find a rather high value of 1.7 × 10¹² M_⊙ when taking the star J1539+0239 into account. This is hyper-velocity star approaching the MW. Gnedin et al. (2010) find a similar value using Jeans modeling of a stellar population in the outer halo. On the other hand, Vera-Ciro et al. (2013) estimate a most likely MW mass of 0.8 × 10¹² M_⊙ by analyzing the Aquarius simulations (Springel et al. 2008) in combination with semi-analytic models of galaxy formation. Watkins et al. (2010) report an only slightly higher value based on the line-of-sight velocities of satellite galaxies (see also Sales et al. 2007), but when they include proper motion estimates they again find a higher mass of 1.4 × 10¹² M_⊙. Using a mixture of stars and satellite galaxies, Battaglia et al. (2005, 2006) also find a low mass below 10¹² M_⊙. McMillan (2011) find an intermediate mass of 1.3 × 10¹² M_⊙, including constraints from photometric data. A further complication comes from the definition of the total mass of the Galaxy which is different for different authors and so a direct comparison of the quoted values has to be made with care. Finally, there is an independent, strong upper limit for the Milky Way mass coming from Local Group timing arguments that estimate the total mass of the combined mass of the Milky Way and Andromeda to 3.2 ± 0.6 × 10¹² M_⊙ (van der Marel et al. 2012b).

In this work we attempt to estimate the mass of the MW through measuring the escape speed at several Galactocentric radii. In this we follow up on the studies by Leonard & Tremaine (1990), Kochanek (1996), and Smith et al. (2007, S07, hereafter). The latter work made use of an early version of the RAdial Velocity Experiment (RAVE; Steinmetz et al. 2006), a massive spectroscopic stellar survey that finished its observational phase in April 2013, and the almost complete set of data will soon be publicly available in the fourth data release (Kordopatis et al. 2013). This tremendous data set forms the foundation of our study.

The escape speed measures the depth of the potential well of the Milky Way and therefore contains information about the mass distribution outside the radius for which it is estimated. It thus constitutes a local measurement connected to the very outskirts of our Galaxy. In the absence of dark matter and a purely Newtonian gravity law, we would expect a local escape speed of $\sqrt{\mathrm{2}}{\mathit{V}}_{\mathrm{LSR}}\mathrm{=}\mathrm{311}$ km s^-1, assuming the local standard of rest V_LSR to be 220 km s^-1 and neglecting the small fraction of visible mass outside the solar circle (Fich & Tremaine 1991). However, the estimates in the literature are much greater than this value, starting with a minimum value of 400 km s^-1 (Alexander 1982) to the currently most precise measurement by S07, who find [498, 608] km s^-1 as a 90% confidence range.

The paper is structured as follows. In Sect. 2 we introduce the basic principles of our analysis. Then we go on (Sect. 3) to describe how we use cosmological simulations to obtain a prior for our maximum likelihood analysis and thereby calibrate our method. After presenting our data and the selection process in Sect. 4, we obtain estimates on the Galactic escape speed in Sect. 5. The results are extensively discussed in Sect. 6, and mass estimates for our Galaxy are obtained and compared to previous measurements. Finally, we conclude and summarize in Sect. 7.

2. Methodology

The basic analysis strategy applied in this study was initially introduced by Leonard & Tremaine (1990) and later extended by S07. They assumed that the stellar system could be described by an ergodic distribution function (DF) f(E) that satisfied f → 0 as E → Φ, the local value of the gravitational potential Φ(r). Then the density of stars in velocity space will be a function n(ν) of speed ν and tend to zero as ν→ ν_esc = (2Φ)^1/2. Leonard & Tremaine (1990) proposed that the asymptotic behavior of n(ν) could be modeled as $\mathit{n}\mathrm{\left(}\nu \mathrm{\right)}\mathrm{\propto }\mathrm{(}{\nu }_{\mathrm{esc}}\mathrm{-}\nu {\mathrm{)}}^{\mathit{k}}\mathit{,}$(1)for ν< ν_esc, where k is a parameter. We should therefore be able to obtain an estimate of ν_esc from a local sample of stellar velocities. S07 used a slightly different functional form, $\mathit{n}\mathrm{\left(}\nu \mathrm{\right)}\mathrm{\propto }\mathrm{(}{\nu }_{\mathrm{esc}}^{\mathrm{2}}\mathrm{-}{\nu }^{\mathrm{2}}{\mathrm{)}}^{\mathit{k}}\mathrm{=}\mathrm{(}{\nu }_{\mathrm{esc}}\mathrm{-}\nu {\mathrm{)}}^{\mathit{k}}\mathrm{(}{\nu }_{\mathrm{esc}}\mathrm{+}\nu {\mathrm{)}}^{\mathit{k}}\mathit{,}$(2)that can be derived if f(E) ∝ E^k is assumed, but, as we see in Sect. 3, results from cosmological simulations are approximated better by Eq. (1).

Currently, the most accurate velocity measurements are line-of-sight velocities, ν_los, obtained from spectroscopy via the Doppler effect. These measurements have typical uncertainties of a few km s^-1, which is an order of magnitude smaller than the typical uncertainties on tangential velocities obtained from the proper motions currently available. Leonard & Tremaine (1990) have already showed with simulated data that because of this, estimates from radial velocities alone are as accurate as estimates that use proper motions as well. The measured velocities ν_los have to be corrected for the solar motion to enter a Galactocentric rest frame. These corrected velocities we denote with ν_∥.

Following Leonard & Tremaine (1990), we can infer the distribution of ν_∥ by integrating over all perpendicular directions: $\begin{array}{ccc}{\mathit{n}}_{\mathrm{\parallel }}\mathrm{(}{\nu }_{\mathrm{\parallel }}\mathrm{|}{r}\mathit{,k}\mathrm{)}& \mathrm{\propto }& \int \mathrm{d}{v}\mathit{n}\mathrm{(}{v}\mathrm{|}{r}\mathit{,k}\mathrm{\left)}\mathit{\delta }\mathrm{\right(}{\nu }_{\mathrm{\parallel }}\mathrm{-}{v}\mathrm{·}{m̂}\mathrm{)}\\ & \mathrm{\propto }& {\left({\nu }_{\mathrm{esc}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{-}\mathrm{\left|}{\nu }_{\mathrm{\parallel }}\mathrm{\right|}\right)}^{\mathit{k}\mathrm{+}\mathrm{1}}\end{array}$(3)again for | ν_∥| < ν_esc. Here δ denotes the Dirac delta function and represents a unit vector along the line of sight.

We do not expect that our approximation for the velocity DF is valid over the whole range of velocities, but only at the high velocity tail of the distribution. We therefore impose a lower limit ν_min for the stellar velocities. A further important requirement is that the stellar velocities come from a population that is not rotationally supported, because such a population is clearly not described by an ergodic DF. In the case of stars in the Galaxy, this means that we have to select for stars of the Galactic stellar halo component.

We now apply the following approach to the estimation of ν_esc. We adopt the likelihood function $\mathit{L}\mathrm{\left(}{\nu }_{\mathrm{\parallel }}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}{\nu }_{\mathrm{esc}}\mathrm{-}\mathrm{|}{\nu }_{\mathrm{\parallel }}\mathrm{|}{\mathrm{)}}^{\mathit{k}\mathrm{+}\mathrm{1}}}{{\mathrm{\int }}_{{\nu }_{\min }}^{{\nu }_{\mathrm{esc}}}\mathrm{d}\nu \hspace{0.17em}\mathrm{(}{\nu }_{\mathrm{esc}}\mathrm{-}\mathrm{|}{\nu }_{\mathrm{\parallel }}\mathrm{|}{\mathrm{)}}^{\mathit{k}\mathrm{+}\mathrm{1}}}\mathrm{=}\mathrm{(}\mathit{k}\mathrm{+}\mathrm{2}\mathrm{)}\frac{\mathrm{(}{\nu }_{\mathrm{esc}}\mathrm{-}\mathrm{|}{\nu }_{\mathrm{\parallel }}\mathrm{|}{\mathrm{)}}^{\mathit{k}\mathrm{+}\mathrm{1}}}{\mathrm{(}{\nu }_{\mathrm{esc}}\mathrm{-}{\nu }_{\min }{\mathrm{)}}^{\mathit{k}\mathrm{+}\mathrm{2}}}$(4)and determine the likelihood of our catalog of stars that have | ν_∥| > ν_min for various choices of ν_esc and k, then we marginalize the likelihood over the nuisance parameter k and determine the true value of ν_esc as the speed that maximizes the marginalized likelihood.

2.1. Non-local modeling

Leonard & Tremaine (1990) (and in a similar form also S07) used Eq. (3) and the maximum likelihood method to obtain constraints on ν_esc and k in the solar neighborhood. This rests on the assumptions that the stars whose velocities are used are confined to a volume that is small compared to the size of the Galaxy and thus that ν_esc is approximately constant in this volume.

In this study we go a step further and take the individual positions of the stars into account. We do this in two slightly different ways: (1) one can sort the data into Galactocentric radial distance bins and analyze them independently. (2) Alternatively, all velocities in the sample are rescaled to the escape speed at the Sun’s position, ${\nu }_{\mathrm{\parallel }\mathit{,i}}^{\mathrm{\prime }}\mathrm{=}{\nu }_{\mathrm{\parallel }\mathit{,i}}\left(\frac{{\nu }_{\mathrm{esc}}\mathrm{\left(}{{r}}_{\mathrm{0}}\mathrm{\right)}}{{\nu }_{\mathrm{esc}}\mathrm{\left(}{{r}}_{\mathit{i}}\mathrm{\right)}}\right)\mathrm{=}{\nu }_{\mathrm{\parallel }\mathit{,i}}\sqrt{\frac{\mathrm{\left|}\mathrm{\Phi }\mathrm{\right(}{{r}}_{\mathrm{0}}\mathrm{\left)}\mathrm{\right|}}{\mathrm{\left|}\mathrm{\Phi }\mathrm{\right(}{{r}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right|}}}\mathit{,}$(5)where r₀ is the position vector of the Sun. For the gravitational potential, Φ(r), model assumptions have to be made. This approach makes use of the full capabilities of the maximum likelihood method to deal with unbinned data and thereby exploit all the information available.

We compare the two approaches using the same mass model: a Miyamoto & Nagai (1975) disk and a Hernquist (1990) bulge for the baryonic components, and for the dark matter halo an original or an adiabatically contracted NFW profile (Navarro et al. 1996; Mo et al. 1998). As structural parameters of the disk and the bulge, we use common values that were also used by S07 and Xue et al. (2008) and are given in Table 1. Apart from its virial mass the NFW profile has the (initial) concentration parameter c as a free parameter. In most cases we fix c by requiring the circular speed at the solar radius, ν_circ(R₀), to be equal to the local standard of rest, V_LSR (after a possible contraction of the halo). As a result our simple model only has one free parameter, its virial mass. For our results from the first approach using Galactocentric bins, we alternatively apply a prior for c taken from the literature to reduce our dependence on the somewhat uncertain value of the local standard of rest.

2.2. General behavior of the method

To learn more about the general reliability of our analysis strategy, we created random velocity samples drawn from a distribution according to Eq. (3) with ν_esc = 550 km s^-1 and k = 4.3. For each sample we computed the maximum likelihood values for ν_esc and k. Figure 1 shows the resulting parameter distributions for three different sample sizes: 30, 100, and 1000 stars. Five thousand samples were created for each value. One immediately recognizes a strong degeneracy between ν_esc and k, and the method tends to find parameter pairs with too low a escape speed. This behavior is easy to understand if one considers the asymmetric shape of the velocity distribution. The position of the maximum likelihood pair strongly depends on the highest velocity in the sample – if the highest velocity is relatively low, the method will favor too low a escape speed. This demonstrates the need for additional knowledge about the power index k as already noticed by S07.

3. Constraints for k from cosmological simulations

Almost all of the recent estimates of the Milky Way mass have made use of cosmological simulations (e.g., Smith et al. 2007; Xue et al. 2008; Busha et al. 2011; Boylan-Kolchin et al. 2013). In particular, those estimates that rely on stellar kinematics (Smith et al. 2007; Xue et al. 2008) make use of the realistically complex stellar velocity distributions provided by numerical experiments. In this study we also follow this approach. S07 used simulations to show that the velocity distributions indeed reach all the way up to the escape speed, but more importantly, from the simulated stellar kinematics, they derived priors on the power-law index k. This was fundamental for their study on account of the strong degeneracy between k and the escape speed shown in Fig. 1 because there were not enough data to break this degeneracy. As we show later, despite our larger data set, we still face the same problem. However, with the advanced numerical simulations available today, we can do a much more detailed analysis.

In this study we make use of the simulations by Scannapieco et al. (2009). This suite of eight simulations comprises resimulations of the extensively studied Aquarius halos (Springel et al. 2008), including gas particles using a modified version of the Gadget-3 code including star formation, supernova feedback, metal-line cooling, and the chemical evolution of the interstellar medium. The initial conditions for the eight simulations were randomly selected from a dark-matter-only simulation of a much larger volume. The only selection criteria were a final halo mass similar to what is measured for the mass of the Milky Way and no other massive galaxy in the vicinity of the halo at redshift zero. We adopted the naming convention for the simulation runs (A – H) from Scannapieco et al. (2009). The initial conditions of simulation C were also used in the Aquila comparison project (Scannapieco et al. 2012). The galaxies have virial masses between 0.7–1.6 × 10¹² M_⊙ and span a wide range of morphologies, from galaxies with a significant disk component (e.g., simulations C and G) to pure elliptical galaxies (simulation F). The mass resolution is 0.22–0.56 × 10⁶ M_⊙. For a detailed description of the simulations, we refer the reader to Scannapieco et al. (2009, 2010, 2011). Details regarding the simulation code can be found in Scannapieco et al. (2005, 2006) and also in Springel (2005).

An important aspect of the Scannapieco et al. (2009) sample is that the eight simulated galaxies have a broad variety of merger and accretion histories, thereby providing a more or less representative sample of Milky Way-mass galaxies formed in a ΛCDM universe (Scannapieco et al. 2011). Our set of simulations is thus useful for the present study, since it gives us information on the evolution of various galaxies, including all the necessary cosmological processes acting during the formation of galaxies and at relatively high resolution.

Also, we note that the same code has been successfully applied to the study of dwarf galaxies (Sawala et al. 2011, 2012), using the same set of input parameters. Despite a mismatch in the baryon fraction (common to almost all simulations of this kind), the resulting galaxies exhibited structures and stellar populations consistent with observations, proving that the code is able to reproduce the formation of galaxies of different masses in a consistent way. Considering that the outer stellar halo of massive galaxies form from smaller accreted galaxies, that we do not need to fine-tune the code differently for different masses once more proves the reliability of the simulation code and its results.

To allow a better comparison to the Milky Way we rescale the simulations to have a circular speed at the solar radius, R₀ = 8.28 kpc, of 220 km s^-1 by the following transformation: $\begin{array}{c}\\ {{r}}_{\mathit{i}}^{\mathrm{\prime }}& \mathrm{=}& {{r}}_{\mathit{i}}\mathit{/}\mathit{f,}& {{v}}_{\mathit{i}}^{\mathrm{\prime }}& \mathrm{=}& {{v}}_{\mathit{i}}\mathit{/}\mathit{f,}\\ {\mathit{m}}_{\mathit{i}}^{\mathrm{\prime }}& \mathrm{=}& {\mathit{m}}_{\mathit{i}}\mathit{/}{\mathit{f}}^{\mathrm{3}}\mathit{,}& {\mathrm{\Phi }}_{\mathit{i}}^{\mathrm{\prime }}& \mathrm{=}& {\mathrm{\Phi }}_{\mathit{i}}\mathit{/}{\mathit{f}}^{\mathrm{2}}\end{array}$(6)where m_i and Φ_i are the mass and the gravitational potential energy of the ith star particle in the simulations. The resulting masses, M₃₄₀, radii, R₃₄₀, and velocities, V₃₄₀, as well as the scaling factors, are given in Table 2. Throughout this work we use a Hubble constant H = 73 km s^-1 Mpc^-1 and define the virial radius to contain a mean matter density 340   ρ_crit, where ρ_crit = 3H²/8πG is the critical energy density for a closed universe¹. The above transformations do not alter the simulation results since they preserve the numerical value of the gravitational constant G governing the stellar motions, the mass density field ρ(r) that governs the gas motions, and the numerical star formation recipe. Only the supernova feedback recipe does not scale in the same way, but since our scaling factors f are close to unity, this is not a major concern.

Since the galaxies in the simulations are not isolated systems, we have to define a limiting distance above which we consider a particle to have escaped its host system. We set this distance to 3R₃₄₀ and set the potential to zero at this radius, which results in distances between 430 and 530 kpc in the simulations. This choice is an educated guess, and our results are not sensitive to small changes, because the gravitational potential changes only weakly with radius at these distances, and in addition, the resulting escape speed is only proportional to the square root of the potential. However, we must not choose too low a value, because otherwise we underestimate the escape speed encoded in the stellar velocity field. On the other hand, we must set the cut-off in a regime where the potential is not yet dominated by neighboring (clusters of) galaxies. Our choice is, in addition, close to half of the distance of the Milky Way and its nearest massive neighbor, the Andromeda galaxy. We test our choice further below. With this definition of the cut-off radius, we obtain local escape speeds at R₀ from the center between 475 and 550 km s^-1.

Now we select a population of star particles belonging to the stellar halo component. In many numerical studies, the separation of the particles into disk and bulge/halo populations is done using a circularity parameter that is defined as the ratio between the particle’s angular momentum in the z-direction² and the angular momentum of a circular orbit either at the particle’s current position (Scannapieco et al. 2009, 2011) or at the particle’s orbital energy (Abadi et al. 2003). A threshold value is then defined to divide disk and bulge/halo particles. We opt for the very conservative value of 0, which means that we only take counter-rotating particles. Practically speaking, this is equivalent to selecting all particles with a positive tangential velocity with respect to the Galactic center. This choice allows us to do exactly the same selection as we will do later with the real observational data, for which we have to use a very conservative value because of the larger uncertainties in the proper motion measurements.

For similar reasons, we also keep only particles in our sample that have Galactocentric distances between 4 and 12 kpc, which reflects the range of values of the stars in the RAVE survey that we use for this study. This further ensures that we exclude particles belonging to the bulge component.

Finally, we set the distance R₀ of the observer from the Galactic center to be 8.28 kpc, choose an azimuthal position φ₀, and compute the line-of-sight velocity ν_∥ ,i for each particle in the sample. We furthermore know the exact potential energy Φ_i of each particle and therefore their local escape speed ν_esc,i is easily computed.

We do this for four different azimuthal positions separated from each other by 90°. The positions were chosen such that the inclination angle with respect to a possible bar is 45°. The corresponding samples are analyzed individually and also combined. These samples are practically statistically independent even though a particle could enter two or more samples. However, because we restrict ourselves the line-of-sight component of the velocities, only in the unlikely case that a particle is located exactly on the line-of-sight between two observer positions, it would gain an incorrect double weight in the combined statistical analysis.

Figure 2 shows the velocity-space density of star particles as a function of 1−ν_∥/ν_esc, and we see that, remarkably, these plots have a reasonably straight section at the highest speeds, just as Leonard & Tremaine (1990) hypothesized. The slopes of these rectilinear sections scatter around k = 3 as we see later.

We also considered the functional form proposed by S07 for the velocity DF; that is, $\mathit{n}\mathrm{\left(}\nu \mathrm{\right)}\mathrm{\propto }\mathrm{(}{\nu }_{\mathrm{esc}}^{\mathrm{2}}\mathrm{-}{\nu }^{\mathrm{2}}{\mathrm{)}}^{\mathit{k}}$. Figure 3 tests this DF with the simulation data. The curvature implies that this DF does not represent the simulation data as well as the formula proposed by Leonard & Tremaine (1990).

If we fit Eq. (3) to the velocity distributions while fixing k to 3, we recover the escape speeds within 6%. This confirms our choice of the cut-off radius for the gravitational potential, 3R₃₄₀, that was used during the definition of the escape speeds.

3.1. The velocity threshold

We now try to find the best value for the lower threshold velocity ν_min. S07 had to use a high threshold value for their radial velocities of 300 km s^-1, because the threshold had an additional purpose, namely to select stars from the non-rotating halo component. If one can identify these stars by other means, the velocity threshold can be lowered significantly. This adds more stars to the sample, thereby putting our analysis on a broader basis. If the stellar halo had the shape of an isotropic Plummer (1911) sphere, the threshold could be set to zero, because for this model the S07 version of our approximated velocity distribution function would be exact. However, for other DFs we need to choose a higher value to avoid regions where our approximation breaks down. Again, we use the simulations to select an appropriate value.

We compute the likelihood distribution of k in each simulation using different velocity thresholds using the likelihood estimator ${\mathit{L}}_{\mathrm{tot}}\mathrm{(}\mathit{k}\mathrm{|}{\nu }_{\min }\mathrm{)}\mathrm{=}\underset{\mathit{i}}{\mathrm{􏽙}}\mathit{L}\mathrm{(}{\nu }_{\mathrm{\parallel }\mathit{,i}}\mathrm{)}\mathit{.}$(7)Figure 4 plots the median values of the likelihood distributions as a function of the threshold velocity. We see a trend toward increasing k for ν_min ≲ 150 km s^-1 and roughly random behavior above. For low values of ν_min, simulation G does not follow the general trend. This simulation is the only one in the sample that has a dominating bar in its center (Scannapieco & Athanassoula 2012), which could contain counter-rotating stars. Given this fact, a likely explanation for its peculiar behavior is that with a low value for the velocity threshold, bar particles start entering the sample, thereby altering the velocity distribution.

Simulation E exhibits a dip around ν_min ≃ 300 km s^-1. A spatially dispersed stellar stream of significant mass is counter-orbiting the galaxy and is entering the sample at one of the observer positions. This is also clearly visible in Fig. 2 as a bump in one of the velocity distributions between 0.2 and 0.3. Furthermore, this galaxy has a rapidly rotating spheroidal component (Scannapieco et al. 2009).

The galaxy in Simulation C has a satellite galaxy very close by. We exclude all star particles in a radius of 3 kpc around the satellite center from our analysis, but there will still be particles entering our samples that originate in this companion and which do not follow the general velocity DF.

All three cases are unlikely to apply for our Milky Way. Our galaxy hosts a much shorter bar, and up to now no signatures of a massive stellar stream were found in the RAVE data (Seabroke et al. 2008; Williams et al. 2011; Antoja et al. 2012). However, it is very interesting to see how our method performs in these rather extreme cases.

We adopt threshold velocities ν_min = 200 km s^-1 and 300 km s^-1. Both are far enough from the regime where we see systematic evolution in the k values (ν_min ≤ 150 km s^-1). For the higher value, we can drop the criterion for the particles to be counter-rotating because we can expect the contamination by disk stars to be negligible (S07) and thus partly compensate for the reduced sample size.

3.1.1. An optimal prior for k

From Fig. 4 it seems clear that the different simulated galaxies do not share exactly the same k, but cover a considerable range of values. Thus in the analysis of the real data we have to consider this whole range. We fix the extent of this range by requiring that it delivers optimal results for all four observer positions in all eight simulated galaxies, so we applied our analysis to the simulated data by computing the posterior probability distribution $\mathit{p}\mathrm{\left(}{\nu }_{\mathrm{esc}}\mathrm{\right)}\mathrm{\propto }{\mathrm{\int }}_{{\mathit{k}}_{\min }}^{{\mathit{k}}_{\max }}\mathrm{d}\mathit{k}\underset{\mathit{i}}{\mathrm{􏽙}}\mathit{L}\mathrm{(}{\nu }_{\mathrm{\parallel }\mathit{,i}}^{\mathrm{\prime }}\mathrm{|}{\nu }_{\mathrm{esc}}\mathit{,k}\mathrm{)}\mathit{,}$(8)where L was defined in Eq. (4), and ${\nu }_{\mathrm{\parallel }\mathit{,i}}^{\mathrm{\prime }}$ is the ith rescaled line-of-sight velocity as defined in Eq. (5). We define the median of p(ν_esc), ${\begin{array}{c}{\mathrm{˜}}^{}\\ \nu \end{array}}_{\mathrm{esc}}$, as the best estimate. For a comparison of the estimates between different simulations, we consider the normalized estimate $\hat{\nu }{}_{\mathrm{esc}}\mathrm{=}{\begin{array}{c}{}^{\mathrm{˜}}\\ \nu \end{array}}_{\mathrm{esc}}\mathit{/}{\nu }_{\mathrm{esc}\mathit{,}\mathrm{true}}$ with ν_esc,true as the true local escape speed in the simulation. By varying k_min and k_max, we identify those values that minimize the scatter in the sample of 32 $\hat{\nu }{}_{\mathrm{esc}}$ values and at the same time leave the median of the sample close to unity. We find very similar intervals for both threshold velocities and adopt the interval $\mathrm{2.3}\mathit{<}\mathit{k}\mathit{<}\mathrm{3.7.}$(9)Reassuringly, this is very close to the lower part of the interval found by S07 (2.7–4.7) using a different set of simulations. The scatter of the $\hat{\nu }{}_{\mathrm{esc}}$ values is less than 3.5% (1σ) for both velocity thresholds. This scatter cannot be completely explained by the statistical uncertainties of the estimates, so there seems to be an additional uncertainty intrinsic to our analysis technique itself. We try to quantify this in the next section.

3.1.2. Realistic tests

One important test for our method is whether it still yields correct results if we have imperfect data and a non-isotropic distribution of lines of sights. To simulate typical RAVE measurement errors we attached random Gaussian errors on the parallaxes (distance^-1), radial velocities, and the two proper motion values with standard deviations of 30%, 3 km s^-1 and 2 mas, respectively. We computed the angular positions of each particle (for a given observer position) and selected only those particles that fell into the approximate survey geometry of the RAVE survey. The latter we define by declination δ < 0° and Galactic latitude | b| > 15°.

Figure 5 shows the resulting distributions of $\hat{\nu }{}_{\mathrm{esc}}$ for the two velocity thresholds. Again, the width of the distributions cannot be solely explained by the statistical uncertainties computed from the likelihood distribution, but an additional uncertainty of ≃4% is required to explain the data in a Gaussian approximation. The distribution for ν_min = 300 km s^-1 in addition exhibits a shift to higher values by ≃3%. Due to the low number statistics the significance of the shift is unclear (~3σ). As we see in Sect. 5, compared to the statistical uncertainties arising when we analyze the real data it would presents a minor contribution to the overall uncertainty and we neglect the shift for this study.

We can go a step further and try to recover the masses of the simulated galaxies using the escape speed estimates. To do this we use the original mass profile of the baryonic components of the galaxies to model our knowledge about the visual parts of the Galaxy and impose an analytic expression for the dark matter halo. As we will do for the real analysis, we try two models: an unaltered and an adiabatically contracted NFW sphere. We adjust the halo parameters, the virial mass M₃₄₀ and the concentration c, to match both boundary conditions, the circular speed, and the escape speed at the solar radius. Figure 6 plots the ratios of the estimated masses and the real virial masses taken from the simulations directly. The adiabatically contracted halo on average overestimates the virial mass by 25%, while the pure NFW halo systematically understates the mass by about 15%. For both halo models we find examples that obtain a very good match with the real mass (e.g., simulation B for the contracted halo and simulation H for the pure NFW halo). However, the cases where the contracted halo yields better results coincide with those cases where the escape speed was underestimated. The colored symbols in Fig. 6 mark the mass estimates obtained using the exact escape speed computed from the gravitational potential in the simulation directly. This reveals that the mass estimates from the two halo models effectively bracket the real mass as expected. We also recover the masses of the three simulations C, E, and G that show peculiarities in their velocity distributions. Only for simulation E and one azimuthal observer position do we completely fail to recover the mass. In this case there is a prominent stellar stream moving along the line of sight.

4. Data

4.1. The RAVE survey

The major observational data for this study comes from the fourth data release (DR4) of RAVE, a massive spectroscopic stellar survey conducted using the 6dF multi-object spectrograph on the 1.2-m UK Schmidt Telescope at the Siding Springs Observatory (Australia). A general description of the project can be found in the data release papers: Steinmetz et al. (2006), Zwitter et al. (2008), Siebert et al. (2011), and Kordopatis et al. (2013). The spectra are measured in the Ca ii triplet region with a resolution of R = 7000. To provide an unbiased velocity sample, the survey selection function was kept as simple as possible: it is magnitude limited (9 < I < 12) and has a weak color cut of J − K_s > 0.5 for stars near the Galactic disk and the bulge.

In addition to the very precise line-of-sight velocities, ν_los, several other stellar properties could be derived from the spectra. Effective temperature T_eff, surface gravity log   g, and metallicity [M/H] of the stars were estimated several times using different analysis techniques (Zwitter et al. 2008; Siebert et al. 2011; Kordopatis et al. 2011). Breddels et al. (2010), Zwitter et al. (2010), Burnett et al. (2011), and Binney et al. (2013) independently used these estimates to derive spectro-photometric distances for a large fraction of the stars in the survey. Matijevič et al. (2012) performed a morphological classification of the spectra and in this way identified binaries and other peculiar stars. Finally Boeche et al. (2011) developed a pipeline to derive individual chemical abundances from the spectra.

The DR4 contains information for nearly 500 000 spectra of more than 420 000 individual stars. The target catalog was also cross-matched with other databases to be augmented with additional information, such as apparent magnitudes and proper motions. For this study we adopted the distances provided by Binney et al. (2013)³ and the proper motions from the UCAC4 catalog (Zacharias et al. 2013).

4.2. Sample selection

The wealth of information in the RAVE survey presents an ideal foundation for our study. Since S07 the amount of available spectra has grown by a factor of 10, and stellar parameters have become available. The number of high-velocity stars has unfortunately not increased by the same factor, which is most likely because RAVE concentrated more on lower Galactic latitudes where the relative abundance of halo stars – which can have these high velocities – is much lower.

We use only high-quality observations by selecting only stars that fulfill the following criteria:

• the stars must be classified as “normal” according to the classification by Matijevič et al. (2012),

• the Tonry-Davis correlation coefficient computed by the RAVE pipeline measuring the quality of the spectral fit (Steinmetz et al. 2006) must be larger than 10,

• the radial velocity correction due to calibration issues (cf. Steinmetz et al. 2006) must be smaller than 10 km s^-1,

• the signal-to-noise ratio (S/N) must be over 25,

• the stars must have a distance estimate by Binney et al. (2013),

• the star must not be associated with a stellar cluster.

The first requirement ensures that the star’s spectrum can be well fitted with a synthetic spectral library and excludes, among other things, spectral binaries. The last criterion removes the giant star (RAVE-ID J101742.6-462715), in particular, from the globular cluster NGC 3201 that would have otherwise entered our high-velocity samples. Stars in gravitationally self-bound structures like globular clusters are clearly not covered by our smooth approximation of the velocity distribution of the stellar halo. We further excluded two stars (RAVE-IDs J175802.0-462351 and J142103.5-374549) because of their peculiar location in the Hertzsprung-Russell diagram⁴ (cf. Fig. 9, upper panel).

In some cases RAVE observed the same target multiple times. In this case we adopt the measurements with the highest S/N, except for the line-of-sight velocities, ν_los, where we use the mean value. The median S/N of the high-velocity stars used in the analysis is 56.

We then convert the precisely measured ν_los into the Galactic rest frame using the following formula: ${\nu }_{\mathrm{\parallel }\mathit{,i}}\mathrm{=}{\nu }_{\mathrm{los}\mathit{,i}}\mathrm{+}\mathrm{(}{\mathit{U}}_{\mathrm{\odot }}\cos {\mathit{l}}_{\mathit{i}}\mathrm{+}\mathrm{(}{\mathit{V}}_{\mathrm{\odot }}\mathrm{+}{\mathit{V}}_{\mathrm{LSR}}\mathrm{\left)}\sin {\mathit{l}}_{\mathit{i}}\mathrm{\right)}\cos {\mathit{b}}_{\mathit{i}}\mathrm{+}{\mathit{W}}_{\mathrm{\odot }}\sin {\mathit{b}}_{\mathit{i}}\mathit{.}$(10)We define the local standard of rest, V_LSR, to be 220 km s^-1, and for the peculiar motion of the Sun, we adopt the values given by Schönrich et al. (2010): U_⊙ = 11.1 km s^-1, V_⊙ = 12.24 km s^-1, and W_⊙ = 7.25 km s^-1.

As mentioned in Sect. 2 we need to construct a halo sample and we do this in the same way as for the simulation data. We compute the Galactocentric tangential velocities, ν_φ, of all stars in a Galactocentric cylindrical polar coordinate system using the line-of-sight velocities, proper motions, distances, and the angular coordinates of the stars. For the distance between the Sun and the Galactic center we use the value R₀ = 8.28 kpc (Gillessen et al. 2009). We perform a full uncertainty propagation using the Monte-Carlo technique with 2000 resamplings per star to obtain the uncertainties in ν_φ. As already done for the simulations we discard all stars with positive median estimate of ν_φ and also those for which the upper end of the 95% confidence interval of ν_φ reaches above 100 km s^-1 to obtain a pure stellar halo sample. This is important because a contamination by stars from the rapidly rotating disk component(s) would invalidate our assumptions made in Sect. 2. Only for this step do we make use of proper motions.

We use the measurements from the UCAC4 catalog (Zacharias et al. 2013) and avoid entries that are flagged as (projected) double star in UCAC4 itself or in one of the additional source catalogs that are used for the proper motion estimate. In such cases we perform the Monte-Carlo analysis with a flat distribution of proper motions between –50 and 50 mas yr^-1, both in right ascension, α and declination, δ.

In principle, we could also use a metallicity criterion to select halo stars. There are several reasons why we did not opt for this. First, we want to be able to reproduce our selection in the simulations. Unfortunately, the simulated galaxies are all too metal-poor compared to the Milky Way (Tissera et al. 2012) and are thus not very reliable in this aspect. This is particularly important in the context of the findings by Schuster et al. (2012), who identified correlations between kinematics and metal abundances in the stellar halo that might be related to different origins of the stars (in-situ formation or accretion). However, despite the unrealistic metal abundances, the formation of the stellar halo is modeled realistically in the simulations including all aspects of accretion and in-situ star formation. In the simulated velocity distributions (Fig. 2), we do not detect any characteristic features that would indicate that the duality of the stellar halo as found by Schuster et al. (2012) is relevant for our study. Second, we would have to apply a very conservative metallicity threshold in order to avoid contamination by metal-poor disk stars. Because of this our sample size would not significantly increase using a metallicity criterion instead of a kinematic one.

It is worth mentioning that the star with the highest ν_∥ = −448.8 km s^-1 in the sample used by S07 (RAVE-ID: J151919.7-191359) did not enter our samples, because it was classified as having problems with the continuum fitting by Matijevič et al. (2012). S07 showed via re-observations that the velocity measurement is reliable, however, the star did not get a distance estimate from Binney et al. (2013). Zwitter et al. (2010) estimate a distance of 9.4 kpc, which owing to its angular position (l,b) = (344.6°,31.4°), would place the star behind and above the Galactic center. The star thus clearly violates the assumption by S07 to deal with a locally confined stellar sample and potentially leads to an overestimate of the escape speed. For the sake of a homogeneous data set, we ignored the alternative distance estimate by Zwitter et al. (2010) and discarded the star.

Figure 7 depicts the velocities ${\nu }_{\mathrm{\parallel }}^{\mathrm{\prime }}$ of all RAVE stars as a function of Galactic longitude l and the two velocity thresholds ν_min = 200 and 300 km s^-1. By selecting for a counter-rotating (halo) population we automatically select against the general sinusoidal trend of the RAVE stars in this diagram. Figure 8 illustrates the spatial distribution of our high-velocity sample. As a result of RAVE avoiding the low Galactic latitudes, stars with small Galactocentric radii are high above the Galactic plane. Furthermore, because RAVE is a southern hemisphere survey, the stars in the catalog are not symmetrically distributed around the Sun. The stars in our high-velocity sample are mostly giant stars with a metallicity distribution centered at −1.25 dex as can be seen in Fig. 9.

4.3. Including other literature data

To increase our sample sizes we have also considered other publicly available and kinematically unbiased data sets. We used the sample of metal-poor dwarf stars collected by Beers et al. (2000, B00 hereafter). The authors also provided the full 6D phase space information including photometric parallaxes. We updated the proper motions by cross-matching with the UCAC4 catalog (Zacharias et al. 2013). We found new values for 2011 stars using the closest counterparts within a search radius of 5 arcsec. For ten stars we found two sources in the UCAC4 catalog closer than 5 arcsec, so discarded these stars. There were a further five cases where two stars in the B00 catalog have the same closest neighbor in the UCAC4 catalog. All these ten stars were discarded as well. Finally, we kept only those stars with uncertainties in the line-of-sight velocity measurement below 15 km s^-1.

There is a small overlap of 123 stars with RAVE, 68 of which have a parallax estimate, ϖ, by Binney et al. (2013) with σ(ϖ) < ϖ. By chance two of these stars entered our high-velocity samples. This, on first glance, very unlikely event is not so surprising if we consider our selection for halo stars, the strong bias towards metal-poor halo stars of the B00 catalog, and the significant completeness of the RAVE survey >50% in the brighter magnitude bins (Kordopatis et al. 2013).

To compare the two distance estimates, we convert all distances, d, into distance moduli, μ = 5log (d/10  pc), because both estimates are based on photometry, so the error distribution should be approximately⁵ symmetric in this quantity. We find that σ_Beers should be about 1.3 mag for the weighted differences (Fig. 10, upper panel) to have a standard deviation of unity. B00 quote an uncertainty of 20% on their photometric parallax estimates, while our estimate corresponds to roughly 60%. We adopt our more conservative value and emphasize that this uncertainty is only used during the selection of counter-rotating halo stars.

We further find a systematic shift by a factor f_dist = 1.5 (δμ = 0.9 mag) between the two distance estimates, in the sense that the B00 distances are greater. Since more information was taken into account to derive the RAVE distances we consider them more reliable. To have consistent distances, we decrease all B00 distances by ${\mathit{f}}_{\mathrm{dist}}^{-1}$ and use these calibrated values in our further analysis.

The data set with the currently most accurately estimated 6D phase space coordinates is the Geneva-Copenhagen survey (Nordström et al. 2004) providing Hipparcos distances and proper motions, as well as precise radial velocity measurements. However, this survey is confined to a very small volume around the Sun and is therefore even more strongly dominated by disk stars than the RAVE survey. We find only two counter-rotating stars in this sample with | ν_∥| > 200 km s^-1, as well as two (co-rotating) stars with | ν_∥| > 300 km s^-1. For the sake of homogeneity in our sample we neglect these measurements.

5. Results

5.1. Comparison to Smith et al. (2007)

As a first check we did an exact repetition of the analysis applied by S07 to see whether we get a consistent result. This is interesting because strong deviations could point to possible biases in the data due to, say, the slightly increased survey footprint of the sky. RAVE contains 76 stars fulfilling the criteria, which is an increase by a factor 5 (3 if we take the 19 stars from the B00⁶ catalog into account). The median values of the distributions are effectively the same (537 km s^-1 instead of 544 km s^-1), and the uncertainties resulting from the 90% confidence interval ([504,574]) are reduced by a factor 0.6 (0.7) for the upper (lower) margin. If we assume that the precision is proportional to the square root of the sample size, we expect a decrease in the uncertainties of a factor ${\mathrm{3}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\simeq }\mathrm{0.6}$.

With the distance estimates available now, we know that this analysis rests on the incorrect assumption that we are dealing with a local sample. If we apply a distance cut d_max = 2.5 kpc onto the data we obtain a sample of 15 RAVE stars and 16 stars from the B00 catalog, and we compute a median estimate of $\mathrm{52}{\mathrm{6}}_{-43}^{\mathrm{+}\mathrm{63}}$ km s^-1. A lower value is expected because the distance criteria mainly removes stars from the inner Galaxy where stars generally have higher velocities. The reason for this is that RAVE is a southern hemisphere survey and therefore observes mostly the inner Galaxy.

5.2. The local escape speed

As described as option (2) in Sect. 2.1, we can estimate for all stars in the catalogs what their radial velocity would be if they were situated at the position of the Sun. We then create two samples using the new velocities. For the first sample we select all stars with rescaled velocities ${\nu }_{\mathrm{\parallel }}^{\mathrm{\prime }}\mathit{>}\mathrm{300}$ km s^-1. S07 shows that such a high velocity threshold predominantly yields halo stars. The resulting sample contains 53 stars (34 RAVE stars), and we refer to it as V300. The second sample has a lower velocity threshold of 200 km s^-1, but stars are preselected, in analogy to the simulation analysis, considering only stars classified as “halo” (Sect. 4.2). This sample we call V200, and it contains 86 stars (69 RAVE stars). Most of the stars are located closer to the Galactic center than the Sun and thus the correction mostly leads to decreased velocity values. In both samples, about 7% of the stars have repeat observations. The maximum difference between two velocity measurements is 2.5 km s^-1.

The resulting likelihood distribution in the (ν_esc,k) parameter plane is shown in the lower panel of Fig. 11. The maximum likelihood pairs for the different samples agree very well, except for the pair constructed from RAVE-only V300 sample, which is located near ν_esc ≃ 410 km s^-1 and k ≃ 0. In all cases a clear degeneracy between k and the escape speed is visible. This was already seen by S07 and reflects that a similarly curved form of the velocity DF over the range of radial velocities available by different parameter pairs.

We go further and compute the posterior probability DF for ν_esc, p(ν_esc) using Eq. (8), which effectively means that we marginalize over the optimized k-interval derived in Sect. 3.1.1. For the medians of these distributions, we obtain higher values than the maximum likelihood value for all samples. This behavior is consistent with our findings in Sect. 2.2 where we showed that the maximum likelihood analysis tends to yield pairs with too low values for k and ν_esc. These median values can be found in Table 3 (“Localized”).

5.3. Binning in Galactocentric distance

For halo stars with original | ν_∥| ≥ 200 km s^-1, we are able to fill several bins in Galactocentric distance r and thereby perform a spatially resolved analysis as described as option (1) in Sect. 2.1. We chose six overlapping bins with a radial width of 2 kpc between 4 and 11 kpc. This bin width is larger than the uncertainties of the projected radius estimates for almost all our sample stars (cf. Fig. 8). The number of stars in the bins are 11, 28, 44, 52, 35, and 8. The resulting median values (again after marginalizing over the optimal k-interval) of the posterior PDF and the 90% confidence intervals are plotted in Fig. 12. The values near the Sun are in very good agreement with the results of the previous section. We find a rather flat escape-speed profile except for the outmost bins, which contain very few stars, though, and thus have large confidence intervals.

6. Discussion

6.1. Influence of the input parameters

The 90% confidence intervals provided by our analysis technique only reflect the statistical uncertainties resulting from the finite number of stars in our samples. In this section we consider systematic uncertainties. In Sect. 3.1.2 we already showed that our adopted interval for the power-law index k introduces a systematic scatter of about 4%.

A further source of uncertainties comes from the motion of the Sun relative to the Galactic center. While the radial and vertical motion of the Sun is known to very high precision, several authors have come to different conclusions about the tangential motion, V_⊙ (e.g., Reid & Brunthaler 2004; Bovy et al. 2012; Schönrich 2012). In this study we used the standard value for V_LSR = 220 km s^-1 and the V_⊙ = 12.24 km s^-1 from Schönrich et al. (2010). We repeated the whole analysis using V_LSR = 240 km s^-1 and compared the resulting escape speeds with the values of our standard analysis (cf. the lower part of Table 3). The magnitudes of the deviations are statistically not significant, but we find systematically lower estimates of the local escape speed for the higher value of V_LSR. The shift is close to 20 km s^-1 and thus comparable to the difference ΔV_LSR. This can be understood if we consider that most stars in the RAVE survey and also in our samples are observed at negative Galactic longitudes and thus against the direction of Galactic rotation (see Fig. 7). In this case correcting the measured heliocentric line-of-sight velocities with a higher solar tangential motion leads to lower ν_∥ which eventually reflects into the escape speed estimate. This systematic dependency is induced by the half-sky nature of the RAVE survey, while this effect might cancel out for an all-sky survey. In contrast, the exact value of R₀ does not influence our results, as long as it is kept within the range of proposed values around 8 kpc.

The quantity with the largest uncertainties used in this study is the heliocentric distance of the stars. In Sect. 4.3 we found a systematic difference between the distances derived for the RAVE stars and for the stars in the B00 catalog. Such systematic shifts can arise for various reasons: e.g., different sets of theoretical isochrones, systematic errors in the stellar parameter estimates, or different extinction laws. Again we repeated our analysis, this time with all distances increased by a factor 1.5, practically moving to the original distance scale of B00. Again we find a systematic shift to lower local escape speeds of the same order as for the alternative value of V_LSR.

We finally also tested the influence of the Galaxy model we use to rescale the stellar velocities according to their spatial position. We changed the disk mass to 6.5 × 10¹⁰   M_⊙ and decreased the disk scale radius to 2.5 kpc, in this way preserving the local surface density of the standard model. The resulting differences in the corrected velocities are below 1%, and no measurable difference in the escape speed estimates was found, illustrating the robustness of our methods to reasonable changes in the Galaxy parameters.

6.2. A critical view of the input assumptions

Our analysis stands and falls with the reliability of our approximation of the velocity DF given in Eq. (1). The conceptual underpinning of this approximation is very weak for four reasons:

• In many analytic equilibrium models of stellar systems at any spatial point, there is a non-zero probability density of finding a star right up to the escape speed ν_esc at that point, and zero probability at higher speeds. For example, the Jaffe (1983) and Hernquist (1990) models have this property but King-Michie models (King 1966) do not: in these models the probability density falls to zero at a speed that is less than the escape speed. There is thus an important counterexample to the proposition that n(ν) first vanishes at ν= ν_esc.

• All theories of galaxy formation, including the standard ΛCDM paradigm, predict that the velocity distribution becomes radially biased at high speed, so in the context of an equilibrium model, there must be significant dependence of the DF on the total angular momentum J in addition to E.

• As Spitzer & Thuan (1972) pointed out, in any stellar system, as E → 0, the periods of orbits diverge. Consequently the marginally-bound part of phase space cannot be expected to be phase mixed. Specifically, stars that are accelerated to speeds just short of ν_esc by fluctuations in Φ in the inner system take arbitrarily long times to travel to apocenter and return to radii where we may hope to study them. As a result different mechanisms populate the outgoing and incoming parts of phase space at speeds ν~ ν_esc: while parts are populated by cosmic accretion (Abadi et al. 2009; Teyssier et al. 2009; Piffl et al. 2011), the outgoing part in addition is populated by slingshot processes (e.g., Hills 1988; Brown et al. 2005) and violent relaxation in the inner galaxy. It follows that we cannot expect the distribution of stars in this portion of phase space to conform to Jeans theorem, even approximately. Equation (1) is founded not just on Jeans theorem but on a very special form of it.

• Counts of stars in the Sloan Digital Sky Survey (SDSS) have most beautifully demonstrated that the spatial distribution of high-energy stars is very non-smooth. The origin of these fluctuations in stellar density is widely acknowledged to be the impact of cosmic accretion, which ensures that at high energies the DF does not satisfy Jeans theorem.

From this discussion it should be clear that to obtain a credible relationship between the density of fast stars and ν_esc, we must engage with the processes that place stars in the marginally bound part of phase space. Fortunately, sophisticated simulations of galaxy formation in a cosmological context do just that. Figure 2 illustrates that Eq. (1) catches the general shape of the velocity DF very well. That we find a relatively small interval for the power-law index k that fits all simulated galaxies with their variety of morphologies argues for the appropriateness of the functional form by Leonard & Tremaine (1990).

The question remains whether the applied simulation technique influences the range of k-values we find, since all eight galaxy models were produced with the same simulation code. In particular, the numerical recipes for so-called sub-grid physics like star formation and stellar energy feedback can have a significant impact on the simulation result as recently demonstrated in the Aquila code comparison project (Scannapieco et al. 2012). However, the main differences were found in the formation of galaxy disks, while in this study we explicitly focus on the stellar halo that was built up from infalling satellite galaxies. Differing implementations of subgrid physics might change the amount of stellar and gas mass being brought in by small galaxies, but it appears unlikely that the phase-space structure of the Galactic halo will change significantly. This view is confirmed by the very similar k-interval found by S07 using simulations with a completely different implementation of subgrid physics.

6.3. Estimating the mass of the Milky Way

We now attempt to derive the total mass of the Galaxy using our escape speed estimates. Doing this we exploit the fact that the escape speed is a measure of the local depth of the potential well $\mathrm{\Phi }\mathrm{\left(}{\mathit{R}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\nu }_{\mathrm{esc}}^{\mathrm{2}}$. A critical point in our methodology is the question whether the velocity distribution reaches up to ν_esc or whether it is truncated at some lower value. S07 used their simulations to show that the level of truncation in the stellar component cannot be more than 10%. However, to test this they first had to define the local escape speed by fixing a limiting radius beyond which a star is considered unbound. The authors state explicitly that the choice of this radius as 3R_vir is rather arbitrary. More stringent would be to state that the velocity distribution in the simulations point to a limiting radius of ~3R_vir beyond which stars do not fall back onto the galaxy or fall back only with significantly altered orbital energies, e.g., as part of an infalling satellite galaxy.

It is not a conceptual problem to define the escape speed as the high end of the velocity distribution in disregard of the potential profile outside the corresponding limiting radius. Then it is important, however, to use the same limiting radius while deriving the total mass of the system using an analytic profile. This means we have to redefine the escape speed to ${\nu }_{\mathrm{esc}}\mathrm{(}\mathit{r}\mathrm{|}{\mathit{R}}_{\max }\mathrm{)}\mathrm{=}\sqrt{\mathrm{2}\mathrm{\left|}\mathrm{\Phi }\mathrm{\right(}\mathit{r}\mathrm{)}\mathrm{-}\mathrm{\Phi }\mathrm{(}{\mathit{R}}_{\max }\mathrm{\left)}\mathrm{\right|}}$(11)where R_max = 3R₃₄₀ seems to be an appropriate value (cf. Sect. 3).

This leads to somewhat higher mass estimates. For example, S07 found an escape speed of 544 km s^-1 and derived a halo mass of 0.85 × 10¹²   M_⊙ for an NFW profile, practically using R_max = ∞. If one consistently applies R_max = 3R_vir, the resulting halo mass is 1.05 × 10¹²   M_⊙, an increase by more than 20%. This is the reason our mass estimates are higher than those by S07 even though we find a similar escape speed. These values represent the masses of the dark matter halo alone, while in the remainder of this study we mean the total mass of the Galaxy when we refer to the virial mass M₃₄₀. Keeping this in mind it is then straightforward to compute the virial mass corresponding to a certain local escape speed. As already mentioned, we use the simple mass model presented in Sect. 2.

In the case of the escape speed profile obtained via the binned data, the procedure becomes slightly more elaborate. We have to compute the escape speeds at the centers of the radial bins R_i and then take the likelihood from the probability distributions PDF_Ri(ν_esc) in each bin. The product of all these likelihoods⁷ is the general likelihood assigned to the mass of the model, i.e., $\mathit{L̂}\mathrm{\left(}{\mathit{M}}_{\mathrm{340}}\mathrm{\right)}\mathrm{=}\underset{\mathit{i}}{\mathrm{􏽙}}{\mathrm{PDF}}_{{\mathit{R}}_{\mathit{i}}}\mathrm{\left(}{\nu }_{\mathrm{esc}}\mathrm{\right(}{\mathit{R}}_{\mathit{i}}\mathrm{|}{\mathit{M}}_{\mathrm{340}}\mathrm{)}\mathrm{)}\mathit{.}$(12)The results of these mass estimates are presented in Table 3. As already seen in Fig. 6 for the simulations the adiabatically contracted halo model always yields higher results than the unaltered halo.

6.4. Fitting the halo concentration parameter

Up to now we have assumed a fixed value for the local standard of rest, V_LSR = 220 km s^-1, to reduce the number of free parameters in our Galaxy model to one. Recently, several authors have found higher values for V_LSR of up to 240 km s^-1 (e.g., Bovy et al. 2012; Schönrich 2012). If we change the parametrization in the model and use the halo concentration c as a free parameter, we can compute the likelihood distribution in the (M₃₄₀,c)-plane in the same way as described in the previous section. Figure 13 plots the resulting likelihood contours for an NFW halo profile (left panel) and the adiabatically contracted NFW profile (right panel).

Navarro et al. (1997) showed that the concentration parameter is strongly related to the mass and the formation time of a dark matter halo (see also Neto et al. 2007; Macciò et al. 2008; Ludlow et al. 2012). With this information we can constrain the range of likely combinations (M₃₄₀,c) further. We use the relation for the mean concentration as a function of halo mass proposed by Macciò et al. (2008). For this we converted their relation for c₂₀₀ to c₃₄₀ to be consistent with our definition of the virial radius. There is significant scatter around this relation reflecting the variety of formation histories of the halos. This scatter is reasonably well fitted by a log-normal distribution with σ_log c = 0.11 (e.g. Macciò et al. 2008; Neto et al. 2007). If we apply this as a prior to our likelihood estimation we obtain the black solid contours plotted in Fig. 13. In the adiabatically contracted case the concentration parameters we are quoting are the initial concentrations before the contraction. Only these are comparable to results obtained from dark matter-only simulations.

The maximum likelihood pair of values (marked by a black “+” in the figure) for the normal NFW halo is M₃₄₀ = 1.37 × 10¹²   M_⊙ and c = 5, which implies a circular speed of 196 km s^-1 at the solar radius. The adiabatically contracted NFW profile yields the same c but a somewhat lower mass of 1.22 × 10¹²   M_⊙. Here the resulting circular speed is only 236 km s^-1.

If we marginalize the likelihood distribution along the c-axis we obtain the one-dimensional posterior PDF for the virial mass. The median and the 90% confidence interval we find to be $\begin{array}{ccc}{\mathit{M}}_{\mathrm{340}}\mathrm{=}\mathrm{1.}{\mathrm{3}}_{-0.3}^{\mathrm{+}\mathrm{0.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}{\mathit{M}}_{\mathrm{\odot }}& & \end{array}$for the unaltered halo profile. For the adiabatically contracted NFW profile we find $\begin{array}{ccc}{\mathit{M}}_{\mathrm{340}}\mathrm{=}\mathrm{1.}{\mathrm{2}}_{-0.3}^{\mathrm{+}\mathrm{0.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}{\mathit{M}}_{\mathrm{\odot }}\mathit{,}& & \end{array}$in both cases almost identical to the maximum likelihood value. It is worth noting that in this approach, the adiabatically contracted halo model yields the lower mass estimate, while the opposite was the case when we fixed the local standard of rest as done in the previous section.

There are several definitions of the virial radius used in the literature. In this study we used the radius which encompasses a mean density of 340 times the critical density for closure in the universe. If one adopts an overdensity of 200, the resulting masses M₂₀₀ increase to $\mathrm{1.}{\mathrm{6}}_{-0.4}^{\mathrm{+}\mathrm{0.5}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}$ M_⊙ and $\mathrm{1.}{\mathrm{4}}_{-0.3}^{\mathrm{+}\mathrm{0.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}$ M_⊙ for the pure and the adiabatically contracted halo profiles, respectively. For an overdensity of 340   Ω₀ ~ 100 (Ω₀ = 0.3 being the cosmic mean matter density), as used, say, by Smith et al. (2007) or Xue et al. (2008), the values even increase to $\mathrm{1.}{\mathrm{9}}_{-0.5}^{\mathrm{+}\mathrm{0.6}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}$ M_⊙ and $\mathrm{1.}{\mathrm{6}}_{-0.3}^{\mathrm{+}\mathrm{0.5}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}$ M_⊙. The corresponding virial radii are $\begin{array}{ccc}{\mathit{R}}_{\mathrm{340}}\mathrm{=}\mathrm{180}\mathrm{\pm }\mathrm{20}\mathrm{kpc}& & \end{array}$for both halo profiles (R₂₀₀ = 225 ± 20 kpc).

6.5. Relation to other mass estimates

As further constraints we can include literature estimates of total masses inside various Galactocentric radii by Xue et al. (2008), Gnedin et al. (2010), and Kafle et al. (2012). Gnedin et al. (2010) obtained an estimate of a mass of 6.9 × 10¹¹ M_⊙ ± 20% within 80 kpc via Jeans modeling using radial velocity measurements of halo stars between 25 and 80 kpc from the Hypervelocity Star Survey. Xue et al. (2008) find a mass within 60 kpc of 4.0 ± 0.7 × 10¹¹ M_⊙ by reconstructing the circular velocity curve using a radial velocities of halo BHB stars from the SDSS combined with cosmological simulations. Kafle et al. (2012) measured a Galactic mass of 2.1 × 10¹¹ M_⊙ within 25 kpc of the Galactic center using a similar data set to Xue et al. (2008), but restricting themselves to stars closer than 25 kpc for which proper motion measurements were available. In this way Kafle et al. (2012) did not have to rely on additional simulation data. We use a 68% confidence interval of [1.8,2.3]  × 10¹² M_⊙ for this last estimate (green shaded area; P. Kafle, priv. comm.). Models fulfilling these constraints are marked in Fig. 14 with colored shaded areas. In the case of the unaltered NFW halo we find excellent agreement with Gnedin et al. (2010) and Kafle et al. (2012), while for the adiabatically contracted model the combination of these estimates favor higher virial masses. The estimate by Xue et al. (2008) is only barely consistent with our results on a 1σ-level for both halo models.

Tests with a different model for the Galactic disk (M_d = 6.5 × 10¹⁰ M_⊙, R_d = 2.5 kpc, similar to the one used by Kafle et al. 2012 and Sofue et al. 2009) resulted in decreased mass estimates (10%), well within the uncertainties. This model changes the values for the circular speed (223 km s^-1 and 264 km s^-1 for the unaltered and the contracted cases, respectively) but not the consistency with the mass estimates by Kafle et al. (2012), Gnedin et al. (2010), or Xue et al. (2008).

Another important constraint for the Galactic halo is the space motion of the satellite galaxy Leo I. Boylan-Kolchin et al. (2013) show that in the ΛCDM paradigm, it is extremely unlikely that a galaxy like the Milky Way has an unbound close-by satellite galaxy. If we take the recent estimates for the Galactocentric distance of 261 ± 13 kpc and the absolute space velocity of $\mathrm{20}{\mathrm{0}}_{-29}^{\mathrm{+}\mathrm{22}}$ km s^-1 (Sohn et al. 2013), we can identify those combinations of M₃₄₀ and c that leave Leo I on a bound orbit. The line separating models in which Leo I is bound from those where it is not bound is also plotted in Fig. 14. All models below this line are consistent with a bound orbit of Leo I. The uncertainties are also shown in the figure in the sense that they mark the ridge lines for the extreme cases that Leo I is slower and closer by 1σ and that it is farther and faster by 1σ. In the case of the unaltered halo profile, our mass estimate is consistent with Leo I being on a bound orbit, while in the contracted case, the mass of the Galaxy would be too low.

Finally, Przybilla et al. (2010) found a star, J1539+0239, with a velocity of $\mathrm{69}{\mathrm{4}}_{-221}^{\mathrm{+}\mathrm{300}}$ km s^-1 at a Galactocentric distance of ~8 kpc moving inwards to the Galaxy. The authors argue that this star should therefore be bound to the Milky Way (see also Irrgang et al. 2013). The star is not in the solar vicinity since its heliocentric distance was measured as 12 ± 2.3 kpc, but its Galactocentric distance is comparable to R₀. We can therefore directly compare our results. Owing to the large uncertainties in the velocity estimate, it is not surprising that our most likely value for ν_esc is consistent with J1539+0239 being on a bound orbit. However, if their median velocity is correct, this star is clearly unbound in our model of the Galaxy and must have obtained its high speed via some other mechanism or be of extragalactic origin.

6.6. Dark matter halo profile

The two halo models, unaltered and adiabatically contracted NFW halos, are rather extreme cases, and the true shape of the Galactic halo is most likely intermediate to these options (Abadi et al. 2010). When we fixed the circular speed at the Sun’s position (as done for the estimates shown in Table 3), the resulting halo masses strongly depended on the shape of the profile. However, when we loosened this constraint using a prior on the halo concentration c (as in Sect. 6.4), our mass estimates became fairly robust to changes of the halo model. In this approach the tension between the constraints coming from the circular speed at the solar radius, and the mass estimates at larger distances are likely to be resolved by an intermediate halo model as proposed by Abadi et al. (2010).

6.7. Future prospects

The ESO cornerstone mission Gaia (Prusti 2012) will soon revolutionize the field of Galactic astronomy. It will deliver the full 6D phase space information for more 100 million stars in the extended solar neighborhood. With these data we will no longer be restricted to the use of radial velocities alone because tangential velocities with similar or even smaller uncertainties will be available. Repeating our analysis with Gaia observations will therefore deliver much more precise results.

On the other hand, we expect that the full complexity of the Galaxy will appear in these data as well. The comparatively sparse RAVE data allowed neglecting many of the details of the Galactic structure, in particular the clumpy nature of the stellar halo. This might no longer be possible with the Gaia data, or in other words, the precision of the estimate might no longer be limited by the data, but by the assumptions in the analysis method itself. It is thus possible that the gain is smaller than one might naively expect if the analysis is repeated in the exact same manner. More robust knowledge about the structure of the inner galaxy obtained, for example, via the analysis of cold tidal streams (Koposov et al. 2010; Sanders & Binney 2013) might make it possible to refine these assumptions.

7. Conclusions

In the present study we analyzed the latest data release of the RAVE survey (fourth data release, Kordopatis et al. 2013), together with additional literature data, to estimate the Galactic escape speed (ν_esc) at various Galactocentric radial bins and, through this, the virial mass of our Galaxy. For it we defined the escape speed as the minimum speed required to reach 3R₃₄₀. To break a degeneracy between our fitting parameters, we had to calibrate our method on a set of cosmological simulations of disk galaxy formation. The 90% confidence interval for our best estimate of the local escape speed is 492 < ν_esc < 587 km s^-1, with a median value of 533 km s^-1.

With our new ν_esc value, we can estimate the virial mass of the Galaxy (baryons and dark matter) by assuming a simple mass model of the baryonic content of the Galaxy and a spherical (adiabatically contracted) NFW halo profile and by using the local standard of rest, V_LSR, as an additional constraint. The resulting values can be found in Table 3.

The value of V_LSR is still under debate. If we loosen our constraint on V_LSR and use a prior on the halo concentration parameter, c, coming from large cosmological simulations we find a most likely value for the virial mass ${\mathit{M}}_{\mathrm{340}}\mathrm{=}\mathrm{1.}{\mathrm{3}}_{-0.3}^{\mathrm{+}\mathrm{0.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}$ M_⊙ for the pure NFW profile and $\mathrm{1.}{\mathrm{2}}_{-0.3}^{\mathrm{+}\mathrm{0.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}$ M_⊙ for an adiabatically contracted halo profile.

In Sect. 6.5 we compared our results to other mass estimates. We found good agreement with estimates based on distant halo stars as well as the space motion of the satellite galaxy Leo I.

---

Acknowledgments

T.P. and M.S. acknowledge support by the German Research Foundation under grant DFG-STE-710/4-3. Funding for RAVE has been provided by the Australian Astronomical Observatory, the Leibniz-Institut für Astrophysik Potsdam (AIP), the Australian National University, the Australian Research Council, the French National Research Agency, the German Research Foundation (SPP 1177 and SFB 881), the European Research Council (ERC-StG 240271 Galactica), the Istituto Nazionale di Astrofisica at Padova, The Johns Hopkins University, the National Science Foundation of the USA (AST-0908326), the W. M. Keck foundation, the Macquarie University, the Netherlands Research School for Astronomy, the Natural Sciences and Engineering Research Council of Canada, the Slovenian Research Agency, the Swiss National Science Foundation, the Science & Technology Facilities Council of the UK, Opticon, Strasbourg Observatory, and the Universities of Groningen, Heidelberg, and Sydney. The RAVE web site is at http://www.rave-survey.org. This research made use of the cross-match service provided by the CDS, Strasbourg. This research made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France.

References

Appendix A: Defining the potential in the simulations

In this section we briefly describe how we consistently define the potential in each of the eight simulations we use in this study. Due to the non-spherical symmetry of the mass distribution in the simulation box the gravitational potential shows a spread at a given galactocentric radius. To obtain a robust estimate of the escape speed we redefine the gravitational potential by assuming that the density profile follows a spherically symmetric NFW profile (Navarro et al. 1997) beyond a radius r_aux: $\begin{array}{ccc}\mathrm{\Phi ̂}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}\mathit{a}\mathrm{\Phi }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{-}\mathrm{median}\mathrm{\left(}\mathrm{\Phi }\mathrm{\right(}{\mathit{r}}_{\mathrm{aux}}\mathrm{\left)}\mathrm{\right)}& \\ & \mathit{a} \mathrm{+}{\mathrm{\Phi }}_{\mathrm{NFW}}\mathrm{\left(}{\mathit{r}}_{\mathrm{aux}}\mathrm{\right)}\mathrm{-}{\mathrm{\Phi }}_{\mathrm{NFW}}\mathrm{\left(}\mathrm{3}{\mathit{R}}_{\mathrm{200}}\mathrm{\right)}& \end{array}$(A.1)where Φ_NFW is the gravitational potential of an NFW sphere with virial mass 10¹⁰   M_⊙ and concentration c′ = 10. The radius r_aux was chosen large enough that the approximation of an NFW sphere is well justified, but small enough that the angular variation of the potential is still small. Furthermore, there must not be any satellite galaxy near this radius in any of the simulations. For our suite of simulations, r_aux = 80 kpc turned out to be a good choice. The value of the concentration parameter is arbitrary since the potential profile is insensitive at large radii for any realistic value of c′. Figure A.1 shows the potential profiles and the approximated profile. The dips in the lower envelopes of the potentials are caused by satellite galaxies orbiting the main halo.

The spread in the potential reflects the fact that these halos live in an anisotropic environment as well as the triaxial shape of the halos. In our simulations we can translate this spread into a maximum deviation from the local escape speed of about 25 km s^-1, which is about 5%. The maximum is 40 km s^-1 in simulation D and minimum is 7 km s^-1 for simulation A.

All Tables

All Figures