The effects of dynamical substructure on Milky Way mass estimates from the high velocity tail of the local stellar halo

TL;DR

This study examines how dynamical substructures like streams affect estimates of the Milky Way's mass derived from local high velocity stars, revealing biases and proposing corrected mass values.

Contribution

It demonstrates the significant impact of phase space substructure on mass estimates and provides a revised Milky Way mass accounting for these biases.

Findings

Substructure causes variability in local velocity distributions.

Mass estimates are biased low by about 20% without correction.

Corrected Milky Way mass estimate is approximately 1.29 x 10^{12} solar masses.

Abstract

We investigate the impact of dynamical streams and substructure on estimates of the local escape speed and total mass of Milky Way-mass galaxies from modelling the high velocity tail of local halo stars. We use a suite of high-resolution, magneto-hydrodynamical cosmological zoom-in simulations, which resolve phase space substructure in local volumes around solar-like positions. We show that phase space structure varies significantly between positions in individual galaxies and across the suite. Substructure populates the high velocity tail unevenly and leads to discrepancies in the mass estimates. We show that a combination of streams, sample noise and truncation of the high velocity tail below the escape speed leads to a distribution of mass estimates with a median that falls below the true value by $\sim 20 \%$, and a spread of a factor of 2 across the suite. Correcting for these biases, we derive a revised value for the Milky Way mass presented in Deason et al. of $1.29 ^{+0.37}_{-0.47} \times 10^{12}$ $\rm M_{\odot}$.

Figures

Full text

††pubyear: 2019††pagerange: The effects of dynamical substructure on Milky Way mass estimates from the high velocity tail of the local stellar halo–References

1 Introduction

The mass of the Milky Way (MW) is one of the most fundamental, basic parameters of our Galaxy, yet remains uncertain by about a factor of 2. Historically, most mass estimates are inferred from stellar tracers of the Galactic potential, which include the kinematics of halo stars (e.g. Leonard & Tremaine, 1990; Piffl et al., 2014; Monari et al., 2018), the kinematics of satellite galaxies (e.g. Callingham et al., 2019) and globular clusters (e.g. Watkins et al., 2019; Posti & Helmi, 2019), modelling of stellar streams (e.g. Malhan et al., 2018) and a timing argument applied to the orbit of the distant satellite Leo I (Zaritsky et al., 1989; Li & White, 2008; Peñarrubia et al., 2016). Weakly bound halo stars have long been thought to be valuable probes of the Galactic potential at large radii, and have, since the 1980s, been used to infer the local Galactic escape speed (e.g. Caldwell & Ostriker, 1981). As more of these stars became observable, statistical techniques were developed to constrain the local escape speed from parametrized modelling of the velocity distribution (Leonard & Tremaine, 1990). The most recent applications of this approach utilize data sets with precise 3-dimensional velocities and distances, and yield mass determinations of $\sim 1.00^{+0.31}_{-0.24}\times 10^{12}$ $\rm M_{\odot}$ (Deason et al., 2019, D19 hereafter). These results compare well with estimates based on other tracers provided by the Gaia second data release, which range from $\sim 1-1.5\times 10^{12}$ $\rm M_{\odot}$.

The popular method of Leonard & Tremaine (1990) relies on the assumption that the stellar distribution function is smooth, well mixed and has a high velocity tail which extends to the local escape speed with the form of a power-law, $f\propto(v_{e}-v)^{k}$. However, the stellar halo contains a lot of highly structured features, including the well known Sagittarius stream, that wraps more than one full ring around the sky (e.g. Sesar et al., 2017; Fardal et al., 2019). In addition to this dominant stream, the local halo contains a wealth of smaller scale, clumpy substructure (Helmi et al., 2017; Koppelman et al., 2018), some of which may arise from recently accreted progenitors (e.g. $\omega\rm Cen$, see Myeong et al., 2018b; Ibata et al., 2019). This range of structure in the stellar halo is supported by recent cosmological simulations (Helmi et al., 2003; Gómez et al., 2013; Cooper et al., 2015; Monachesi et al., 2019; Simpson et al., 2019), which show that progenitors of the stellar halo can produce streams in the local phase space distribution, particularly if they undergo several pericentric passages as dynamical friction drags them to the galaxy centre. It is reasonable to expect that such a distribution would: i) be populated unevenly; ii) vary as a function of position within the Galaxy; and iii) lead to potentially significant biases in the escape speed estimate and therefore the total mass estimate of the MW.

To understand the effects of substructure on MW mass estimates, we use a suite of high-resolution cosmological simulations of the formation of MW-mass galaxies in the $\Lambda$CDM paradigm (Grand et al., 2017). The high resolution allows for exquisite sampling of the gravitational potential such that a range of stellar substructure in local volumes is captured. We apply the analysis of D19 to local volumes centred on four equidistant solar positions in the disc of each simulated galaxy, and derive escape speeds and MW mass estimates. We find that the estimates for the escape speed and total mass vary significantly between positions and individual haloes, and in some cases deviate strongly from the true values. We link these trends with dynamical substructure and the high velocity tail, quantify the biases and discuss the implications for the MW.

2 Simulations

We analyse local samples of accreted halo stars in a suite of high resolution, cosmological magneto-hydrodynamical simulations for the formation of the MW (from the auriga project, Grand et al., 2017; Grand et al., 2018). Our suite comprises simulations from two different resolution levels: level 3 ($\sim 6\times 10^{3}$ $\rm M_{\odot}$ per baryonic element; softening length of $184$ pc after $z=1$) and level 4 ($\sim 5\times 10^{4}$ $\rm M_{\odot}$ per baryonic element; softening length of $369$ pc after $z=1$), to maximize statistics. The haloes range between $0.5$-$2\times 10^{12}\rm M_{\odot}$ in total mass ($M_{200}$), which we define as the mass contained inside the radius at which the mean enclosed mass volume density equals 200 times the critical density for closure. Each halo was initially selected from the $z=0$ snapshot of a parent dark matter only cosmological simulation of comoving periodic box size 100 Mpc, with the standard $\Lambda$CDM cosmology. The adopted cosmological parameters are $\Omega_{m}=0.307$, $\Omega_{b}=0.048$, $\Omega_{\Lambda}=0.693$ and a Hubble constant of $H_{0}=100h$ km s*-1* Mpc*-1*, where $h=0.6777$, taken from Planck Collaboration (2014). At $z=127$, the resolution of the dark matter particles of this halo is increased and gas is added to create the initial conditions of the zoom, which is evolved to present day with the magneto-hydrodynamics code arepo (Springel, 2010).

The simulations include a comprehensive galaxy formation model (Grand et al., 2017), which treats many relevant processes, including gravity, gas cooling, star formation, mass and metal return from stellar evolutionary processes, energetic stellar and AGN feedback, and magnetic fields. Large cosmological box simulations (e.g. Pillepich et al., 2018; Nelson et al., 2018) have demonstrated that a similar galaxy formation model reproduces many of the global properties of the observed galaxy population, such as the galaxy stellar mass function, galaxy sizes, the cosmic star formation rate density and galaxy morphological mix. The combination of a globally successful model with the ability to resolve internal galactic structures, such as bars and spiral arms, makes the simulation suite ideal to study Galactic dynamics within the full cosmological setting.

3 Estimating the escape speed and mass of the Galaxy

The mass estimation procedure we test is based on that of D19 and consists of two parts: i) the estimation of the local escape speed from local, accreted star particles, which is sensitive to the outer halo; and ii) the estimation of the total mass using the circular velocity at the Solar radius as an additional constraint. We briefly summarise the method below, and refer the reader to D19 for more details.

For a given solar position, we select all accreted star particles within a 3 kpc sphere. We assume that the velocity distribution follows a power law: $f\propto(v-v_{e})^{k}$, where $v$ is the total velocity of the particle in the Galactic rest frame, $v_{e}$ is the escape velocity at the position of the star and $k$ is the power law slope. We further assume that the escape velocity varies as a function of galactocentric radius as (e.g. Williams et al., 2017): $v_{e}=v_{e,0}(r/r_{0})^{\gamma/2}$, where $r_{0}=8$ kpc is the solar radius, and $v_{e,0}$ is the escape speed at the position of the Sun. We use a log-likelihood analysis to evaluate probability density functions for $v_{e,0}$, $\gamma$ and $k$, on uniform grids in the range $v_{e,0}\in[300,900]$, $\gamma\in[0,1]$ and $k\in[k_{\rm min},k_{\rm max}]$. Here, $k_{\rm min}$ and $k_{\rm max}$ are defined to be equal to $k_{\rm av}\pm 1$, where $k_{\rm av}$ is measured from all accreted star particles in the radial range $4<r/\rm kpc<12$ for each simulation, which mimics the construction of the prior on $k$ adopted by D19 for the MW. For further consistency, we adopt the prior $P(v_{e,0})\propto 1/v_{e,0}$.

To constrain the dark matter mass, $M_{\rm DM,200}$, we first calculate the escape speed at the solar radius: $v_{e}(r_{0})=\sqrt{2(\Phi(r_{0})-\Phi(2R_{200}))}$, where $\Phi$ is the total gravitational potential from all baryonic and dark material. We assume an NFW profile (Navarro et al., 1997), and calculate its potential at $r_{0}$ and $2R_{200}$ for a range of values for $M_{\rm DM,200}$ and concentration, $c$, that span uniform grids of $\log_{10}M_{\rm DM,200}\in[11.5,12.5]$ and $c\in[1,30]$. The baryonic contribution to the potential is calculated from the simulation data. We then evaluate the probability distribution of the dark halo parameters using the escape speed derived from the high velocity tail, after marginalising over $k$ and $\gamma$, the latter of which is poorly constrained given the narrow radial range of stars under consideration. To yield a final mass estimate, we combine the probability distribution based on the escape speed with constraints on $M_{\rm DM,200}$ and $c$ from the circular velocity, $v_{c}(r_{0})=\sqrt{GM(<r_{0})/r_{0}}$.

Our method differs from that of D19 in two main ways: firstly, they model the baryonic component with fixed analytic expressions, whereas we calculate the circular velocity and escape speed from the simulated particle data. Secondly, D19 consider counter-rotating stars from the Gaia DR2, whereas we consider local, accreted star particles at four solar positions - each located at a radius of 8 kpc in the disc midplane, and equidistant in galacto-centric azimuth - per simulation. Our selection of accreted star particles is appropriate given that counter-rotating stars are mainly ex-situ, with the additional benefit of maximising the number of solar vicinities that contain at least 240 accreted star particles within a 3 kpc sphere (or 5 kpc, if a 3 kpc sphere contains fewer particles), which was the number of stars selected by D19 in the same volume.

4 Results

4.1 Illustrative cases

We begin by calculating the true escape speed in each simulation, for which we have the total gravitational potential of each particle available. We calculate the true escape speed from the gravitational potential at $r_{0}$ and the mean graviational potential of particles in a thin spherical shell at $2R_{\rm 200}$. To assess the accuracy with which the method described above is able to recover the true escape speed from a smooth, well-mixed distribution, we apply the method to all dark matter particles contained within each volume, which are expected to be virialized. The upper panels of Fig. 1 show the velocity distributions of these dark matter particles within each volume for simulation Au 23. As expected, the distributions are smooth and are populated all the way up to the true escape speed, which is recovered within 1%. This verifies that the method works well for smooth distributions that reach the escape speed.

Contrary to the dark matter particles, the upper panels of Fig. 1 show that velocity distributions of accreted star particles exhibit multiple bumpy features at both low ($\sim 300$ $\rm km\,s^{-1}$) and high ($>400$ $\rm km\,s^{-1}$) velocities. The shapes of the distributions, as well as the estimated escape speeds, vary with solar position: the escape speed shows a variation of $\sim 70$ $\rm km\,s^{-1}$, with the lowest and highest estimates bracketing the true escape speed. This variation in escape speed estimates between the different solar positions can be understood from inspection of the distribution of these particles in phase space. We define the total specific energy of a particle as $E_{\rm orb}=0.5v^{2}+\Phi$, where $v$ is the velocity of the particle, $\Phi$ is the gravitational potential. We define the magnitude of the vertical component of the specific angular momentum of particles, in a coordinate frame whose $z$-axis is aligned with the spin-axis of the disc, as $L_{z}=(\mathbf{x}\times\mathbf{v})\cdot\mathbf{\hat{z}}$, where $\mathbf{x}$ and $\mathbf{v}$ are the three dimensional position and velocity of the particle. The lower panels of Fig. 1 show the distribution of selected accreted star particles in the $E_{\rm orb}$-$L_{z}$ plane at each of the solar positions. These diagrams reveal a highly structured phase space distribution, characterised by a variety of distinct stream-like features such as clumps at high binding energies, and narrow structures at low binding energies. The overestimate in the $\theta_{\rm ref}=90$ volume is largely driven by the under-populated velocities in the range $\sim 320$-$400$ $\rm km\,s^{-1}$ relative to the peaks at lower and higher velocities, which skews the likelihood function to estimate too high an escape speed. We find that removing the most weakly bound stream stars in this volume (yellow box in Fig. 1) alleviates the skew, but also removes the highest velocity stars from the distribution, resulting in an underestimate of the escape speed. This illustrates the extreme sensitivity of the escape speed estimate to the details of dynamical substructure.

We emphasise that almost all of the substructure in each of the lower panels of Fig. 1 originates from the same progenitor system that merged with the primary galaxy approximately 3 Gyr ago. The peak stellar mass of this object (i.e. the highest stellar mass it attained over its history) is $\sim 5.4\times 10^{9}$ $\rm M_{\odot}$, and it experienced two pericentric passages before it merged. The abundance and diversity of substructure in this example can therefore be explained by repeated stellar stripping events of this satellite that occur during infall; stars stripped at different times and places, as well as multiple wraps of orbits from a single stripping event, have different binding energies and populate different regions of phase space as streams. The upper panels of Fig. 1 show that this debris is responsible for the bumps in the velocity distribution at each solar position.

In contrast to highly structured phase space distributions, we show in Fig. 2 an example of a smooth local halo star distribution. All four solar volumes appear similarly smooth, and as a result, the high velocity tail is relatively smooth. As such, the estimated escape speed is consistent across all volumes, although it falls $\sim 20$ $\rm km\,s^{-1}$ short of the true value. This halo has experienced a quiet merger history since $z\sim 1$ relative to the halo shown in Fig. 1.

4.2 Simulation sample

We now investigate the estimated escape speeds and total halo masses with respect to their true values for a suite of simulations. We include haloes that contain at least 240 accreted star particles in solar-like vicinities (Section. 3), which results in a total of 30 haloes: 23 at level 4 resolution; and 7 at level 3 resolution. We have verified that the results do not change if the level 3 simulations are replaced with their level 4 counterparts. The number of accreted star particles in each volume varies between several hundred to several thousand across the suite. To match the Gaia DR2 sample size of D19, we draw random samples of 240 accreted star particles from each full sample volume until they are used up (e.g., a volume with 1000 star particles is sampled 4 times). This approach increases the effective size of our sample of local accreted star volumes to a total of 892, which yields better statistics and incorporates the same degree of noise as in the D19 sample.

In Fig. 3, we show the range of estimated escape speeds (left panel) and total masses (right panel) for the random samples (crosses) for each simulation. The distributions of the escape speed and total mass estimates are similarly smooth and non-Gaussian, with median values of $-0.03$ dex and $-0.11$ dex, respectively. The 16th-84th percentile range is $\sim 0.1$ dex and $\sim 0.3$ dex for the escape speed and total mass, respectively, indicating the sensitivity of the mass estimate to that of the escape speed. The scatter between positions in each halo is comparable to that between simulations, which indicates that noise (from 240 star particles) and variation of substructure between positions is at least as important as the merger history for the estimate accuracy. For many haloes, the range of estimates of both quantities from all accreted particles at each solar position (plusses) is smaller (but significant) than that of the random samples, which illustrates the joint effect of noise and substructure on the scatter. Fig. 3 shows also the estimated mass distribution if the true escape speed is used to constrain $M_{\rm DM,200}$. This shows that the assumption of a spherical NFW profile leads to a negligible bias and a scatter of $\sim 0.1$ dex.

For the full accreted particle data at each position, we discern a range of accuracies for both escape speed and mass that reflect the significant variation in dynamical structure (from smooth to highly structured) both between positions within the same halo and between different haloes. For example, haloes with the smoothest structure (as in Fig. 2) at all positions (e.g., haoles 6, 9 and 26) typically either return good estimates or underestimates with a low degree of scatter between positions. Evidently, the high velocity tail does not extend to the escape speed for the latter. In contrast, haloes with the largest variations (e.g., haloes 16, 21) tend to have prominent dynamical substructure (as in Fig. 1) created by systems of peak stellar mass greater than $10^{9}$ $\rm M_{\odot}$ that merged within the last several gigayears of evolution after several pericentric passages.

5 Discussion

The local stellar halo in the MW is thought to be dominated by a single component, characterised by a radially anisotropic velocity distribution (“the Sausage”), sculpted long ago from a single progenitor (Belokurov et al., 2018; Helmi et al., 2018). Fattahi et al. (2019) found that the auriga haloes with analogous properties formed most of their stellar halo at early times and experienced their last major merger between 6 and 10 Gyr ago (haloes 5, 9, 10, 17, 18, 22, 26). For our full local halo samples (crosses in Fig. 3), the mass estimates are consistently within $\sim 10\%$ of the true mass at each position for some haloes (6, 9), whereas others produce mass estimates up to $\sim 40\%$ different from (usually lower than) the true mass. The latter are consistent with the high velocity tail truncating below the true escape speed. This may be understood from the orbits of early merger ($z>1$) progenitors, which need to be sufficiently tight to be captured by the main galaxy, such that its stars are not expected to go out as far as today’s $2R_{200}$ after disruption. This is consistent with the mass estimates in this work falling $20\%$ below, on average, the true mass. This reasoning applies also to a similar underestimate reported in Sec. 5.2 of Wilkinson & Evans (1999) related to streams that fall short of the escape speed.

With respect to dynamical substructure, the “Sausage” haloes with the largest spreads in Fig. 3 (5, 10) do contain moderate substructure at some positions, despite their early halo assembly. More importantly, the continual discovery of dynamical substructure of halo stars in the local vicinity (e.g. Helmi et al., 1999, 2017; Myeong et al., 2018a, b; Koppelman et al., 2018; Ibata et al., 2019) indicates that the scatter (and therefore, uncertainty) in the mass measurements linked to streams found in this paper are an important consideration for the MW. This is pertinent for future larger volume samples, especially if they include any member stars of the richly structured Sagittarius stream (Sesar et al., 2017). It is therefore clear that, in addition to larger and more accurate future datasets, a complete mapping of the streams and assembly history of the stellar halo is required to fully understand the extent of these biases. For the present situation, we may use the biases found in Fig. 3 (median and spread) to revise the MW mass derived in D19 to $1.29^{+0.37}_{-0.47}\times 10^{12}$ $\rm M_{\odot}$, which agrees well with estimates based on satellites (Callingham et al., 2019) and globular clusters (Posti & Helmi, 2019).

6 Conclusions

We have used a total of 30 magneto-hydrodynamical cosmological zoom-in simulations for the formation of MW-mass galaxies to investigate the connection between the phase-space substructure and high velocity tail of accreted stars in the vicinity of solar-like positions. We investigate how this substructure and velocity distribution varies between different solar-like positions and across the simulation set. We have applied a common method for estimating the escape speed and total mass of each halo from the high velocity tail D19. Our conclusions are summarised as follows:

• •

The properties of dynamical substructure depend on the merger history and vary dramatically (from smooth to highly structured) between positions in individual haloes and between different haloes. The most significant phase space substructure originates from massive progenitors, that had: peak stellar masses $\gtrsim 10^{9}$ $\rm M_{\odot}$; and several pericentric passages before merging at late times, whereas the quietest merger histories have the smoothest phase space structure.

• •

The method to estimate the escape velocity and halo mass from the high-velocity tail is extremely sensitive to the dynamical structure of halo stars. Dynamical substructure (streams) leads to un-evenly populated velocity distributions that act to increase the scatter of mass estimates between positions within individual haloes. Smooth phase space distributions are conducive to better estimates, although they have a tendency to underestimate the mass owing to the truncation of the high velocity tail below the true escape speed. However, the effects of noise and streams lead to significant scatter in the estimated quantities, even in cases for which the velocity distribution is relatively smooth.

• •

The distribution of escape speed and mass estimates from all simulations and positions have medians lower than the true values by $7\%$ and $20\%$, and spreads of factors of $\sim 1.3$ and $\sim 2$, respectively. Based on these biases, the MW mass estimate derived in D19 is revised to $1.29^{+0.37}_{-0.47}\times 10^{12}$ $\rm M_{\odot}$.

Acknowledgements

RG thanks Guinevere Kauffmann, Wilma Trick and Vasily Belokurov for useful comments. AD is supported by a Royal Society University Research Fellow-ship. AD also acknowledges the support from the STFC grant ST/P000541/1. FAG acknowledges financial support from CONICYT through the project FONDECYT Regular Nr. 1181264, and funding from the Max Planck Society through a Partner Group grant. FM is supported through the Program “Rita Levi Montalcini” of the Italian MIUR.