The Mass Profile of the Milky Way to the Virial Radius from the Illustris Simulation
Corbin Taylor, Michael Boylan-Kolchin, Paul Torrey, Mark Vogelsberger, Lars Hernquist
MMNRAS , 3483–3493 (2016) Preprint 28 July 2016 Compiled using MNRAS L A TEX style file v3.0
The Mass Profile of the Milky Way to the Virial Radius from theIllustris Simulation
Corbin Taylor (cid:63) , Michael Boylan-Kolchin † , Paul Torrey , , Mark Vogelsberger ,and Lars Hernquist Department of Astronomy, University of Maryland, 1113 Physical Sciences Complex (Building 415), College Park, MD 20742-2421, USA Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712-1205, USA Department of Physics, Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA, 02138, USA
28 July 2016
ABSTRACT
We use particle data from the Illustris simulation, combined with individual kinematic con-straints on the mass of the Milky Way (MW) at specific distances from the Galactic Centre,to infer the radial distribution of the MW’s dark matter halo mass. Our method allows us toconvert any constraint on the mass of the MW within a fixed distance to a full circular velocityprofile to the MW’s virial radius. As primary examples, we take two recent (and discrepant)measurements of the total mass within 50 kpc of the Galaxy and find that they imply verydifferent mass profiles and stellar masses for the Galaxy. The dark-matter-only version ofthe Illustris simulation enables us to compute the effects of galaxy formation on such con-straints on a halo-by-halo basis; on small scales, galaxy formation enhances the density rela-tive to dark-matter-only runs, while the total mass density is approximately 20% lower at largeGalactocentric distances. We are also able to quantify how current and future constraints onthe mass of the MW within specific radii will be reflected in uncertainties on its virial mass:even a measurement of M ( <
50 kpc) with essentially perfect precision still results in a 20%uncertainty on the virial mass of the Galaxy, while a future measurement of M ( <
100 kpc) with 10% errors would result in the same level of uncertainty. We expect that our techniquewill become even more useful as (1) better kinematic constraints become available at largerdistances and (2) cosmological simulations provide even more faithful representations of theobservable Universe.
Key words:
Galaxy: fundamental parameters–Galaxy: halo–Galaxy: structure–dark matter.
While living within the Milky Way (MW) galaxy does have itsvirtues, easily and accurately determining the mass distribution ofthe Galaxy’s dark matter halo is not one of them. This is not forlack of trying, naturally; a variety of techniques have been craftedfor just this purpose, and multiple classes of kinematic tracers areavailable.The difficulty in measuring the MW’s mass distribution is two-fold. First, only line-of-sight information is available for the vastmajority of kinematic measurements. While great strides are be-ing made in measuring the proper motions of both individual stars(Cunningham et al. 2015) and dwarf galaxies (e.g., Piatek et al.2007; Sohn et al. 2013; van der Marel et al. 2014; Pryor et al. 2015)at large Galactocentric distances in the MW’s halo, the number of (cid:63) [email protected] † [email protected] tracers at ∼ − kpc with full 6D phase space informationwill remain small even in the Gaia era (de Bruijne et al. 2014). Per-haps more importantly, the level of precision desired for the MW’smass is simply higher than is the case for other galaxies. Whereasa factor of ± uncertainty in the mass of a typical galaxy’s halowould be considered an excellent measurement, it is often thoughtof more as an embarrassment in the case of the MW.For example, if we take a dark matter halo mass of M (cid:12) as a fiducial estimate for the MW, changes by a factor of 2 in eitherdirection are the difference between: (1) an implied conversion ef-ficiency of baryons into stars of ≈ (at M = 5 × M (cid:12) )and (at × M (cid:12) ); (2) eliminating the too-big-to-fail prob-lem (Boylan-Kolchin et al. 2011b, 2012) and severely exacerbatingit (Wang et al. 2012; Vera-Ciro & Helmi 2013; Jiang & van denBosch 2015); and (3) placing the Large Magellanic Cloud and theLeo I dwarf spheroidal on unbound versus bound orbits (Kallivay-alil et al. 2006, 2013; Besla et al. 2007; Boylan-Kolchin et al. 2013). c (cid:13) a r X i v : . [ a s t r o - ph . C O ] J u l C. Taylor et al.
Our understanding of the MW in cosmological context relies on ourability to know its mass to high precision.While uncertainties are most pronounced in the outer darkmatter halo of the MW, where there are few tracers of the totalmass, they also persist at small Galactocentric distances: there aredisagreements about the mass within the solar circle at the 25%level (e.g., Bovy et al. 2012; Sch¨onrich 2012). At 40–80 kpc, esti-mates differ at the 50% level (see, e.g., Williams & Evans 2015).In this paper, we take an alternate approach to constraining themass distribution of the MW. Cosmological hydrodynamic simu-lations are now producing galaxies that match a variety of obser-vations both for statistical samples of galaxies and for individualgalaxies themselves. In particular, both the Illustris (Vogelsbergeret al. 2014b) and Eagle (Schaye et al. 2015) simulations use ∼ particles within ∼ Mpc boxes, meaning they contain thousandsof haloes with masses comparable to that of the MW, each with ofthe order 1 million particles within the virial radius. The successesof these models, and the underlying successes of the Λ cold darkmatter ( Λ CDM) model, motivate using the results of cosmologicalsimulations to constrain the mass distribution of the MW.There are a number of ways one could use cosmological simu-lations for this purpose. Indeed, several previous works on the massof the MW have used cosmological simulations in some capacity.One possibility is to use dark matter haloes from large cosmolog-ical simulations as point particles and calibrate the timing argu-ment (Kahn & Woltjer 1959) for measuring the total mass of theMW (Li & White 2008; Gonz´alez et al. 2014). Alternately, proper-ties of satellites from cosmological simulations can be compared tothose of MW satellites such as the Magellanic Clouds, yielding es-timates of the virial mass of the MW (Boylan-Kolchin et al. 2011a;Busha et al. 2011; Gonz´alez et al. 2013; Fattahi et al. 2016). Yet an-other possibility is to use individual, high-resolution simulations ofMilky Way-sized haloes in conjunction with kinematic informationabout dwarf satellites of the MW (Boylan-Kolchin et al. 2013; Bar-ber et al. 2014). Cosmological hydrodynamic simulations of indi-vidual MW-mass haloes have also been used to calibrate kinematicanalyses of tracer populations in order to measure the mass of theMW (Xue et al. 2008; Rashkov et al. 2013; Piffl et al. 2014; Wanget al. 2015).Our approach is to use importance sampling in ahomogeneously-resolved, large-volume cosmological simula-tion, weighing each simulated halo by its level of consistencywith the MW; for a clear description of this technique applied tocosmological simulations, see Busha et al. (2011). By taking anyindividual constraint and using it to perform importance samplingfrom simulations, we can find the mass distributions of haloes thatare consistent with the imposed constraint. An advantage of thistechnique is that it allows us to easily map different constraints,with different errors, on to mass distributions for the MW and itsdark matter halo.Variants of importance sampling have been used to measurethe mass of the MW (Li & White 2008; Boylan-Kolchin et al.2011a; Busha et al. 2011; Gonz´alez et al. 2014). However, previ-ous work has generally focused on using dark-matter-only (DMO)simulations to measure the total (virial) mass of the MW. With hy-drodynamic simulations, we are able to make two improvements.First, we are able to measure the mass distribution of the MW insimulations that self-consistently model the effects of galaxy for-mation on the dark matter haloes of galaxies. Secondly, we are ableto compare our constraints directly to those obtained from DMOsimulations, as a DMO version of Illustris is also publicly available.By matching objects between the two simulations, we are able to investigate, in detail, the effects of baryonic physics on inferencesof the mass distribution of the MW from cosmological simulations.We generally use the mass within 50 kpc as our primary con-straint, as this is approximately the largest radius where stellar kine-matic tracers are found in large enough numbers to facilitate a massmeasurement. We also provide estimates for how a measurement ofthe mass within 100 kpc – which future surveys may provide – willimprove our knowledge of the mass distribution at even larger radii.This paper is structured as follows. Section 2 describes ourbasic approach, provides information about the Illustris simulation,and describes our primary analysis of the simulation. Section 3 con-tains our main results regarding the mass distribution of the MW asderived from haloes taken from the Illustris simulation. We alsoquantify how inferences on the enclosed mass at large scales (at250 kpc and various spherical overdensity values) depend on themeasured mass within 50 kpc and quantify the stellar masses ofgalaxies having haloes consistent with the adopted mass constraint.A discussion of our results and prospects for future improvementsis given in Section 4; our primary conclusions are given in Sec-tion 5. Throughout this paper, error bars give 68% confidence in-tervals unless otherwise noted.
Our analysis is based on the Illustris suite of cosmological simula-tions (Vogelsberger et al. 2014b), which consists of paired hydro-dynamic and DMO simulations at three different resolution levels.Each simulation uses a periodic box of length h − Mpc andan initial redshift of z = 127 . The highest resolution simulation,Illustris-1, uses dark matter particles and an equal number ofhydrodynamic cells initially, with a spatial resolution of h − kpcfor the dark matter. The DMO version of this simulation, Illustris-Dark-1, uses identical initial conditions but treats the baryonic com-ponent as collisionless mater. Two lower resolution simulation ofthe same volume, Illustris-2 and Illustris-3, were also performed,with 8 and 64 times fewer particles, respectively. The backgroundcosmology for all of the simulations was chosen to be consistentwith Wilkinson Microwave Anisotropy Probe-9 results (Hinshawet al. 2013): Ω m , = 0 . , Ω Λ , = 0 . , Ω b , = 0 . , σ = 0 . , n s = 0 . , and h = 0 . . Haloes and subhaloes inthe Illustris simulations were identified using a friends-of-friendsalgorithm followed by SUBFIND (Springel et al. 2001). For fur-ther information about the Illustris suite, including details aboutthe implementation of galaxy formation physics, see Vogelsbergeret al. (2013, 2014a).Using Illustris to inform our understanding of the mass distri-bution of the MW requires calculations of the mass profiles of anunbiased sample of dark matter haloes within the simulation. Al-though the halo catalogues provide the centres for each halo (weonly consider central halos, not subhaloes, as possible centres), abrute-force calculation of the mass profile for each halo is pro-hibitively expensive, as it requires repeated searches through the ∼ particles of the simulation. We instead use a K-D tree algo-rithm, taking into account the periodic boundary conditions of thesimulation volume. The algorithm was verified against brute-forcecalculations applied to Illustris-3 and Illustris-2. , 3483–3493 (2016) ass of the Milky Way from Illustris In addition to considering the mass within spherical apertures,we also compute spherical overdensity masses with respect to threecommon overdensity choices: M , c (measured with respect to ρ crit ), M , m (measured with respect to ρ m ≈ ρ crit for the Illustris cosmology at z = 0 ), and M vir (measured with re-spect to ∆ vir ρ crit ; for the Illustris cosmology at z = 0 , ∆ vir ≈ ;Bryan & Norman 1998). Our basic framework is to consider the Illustris simulation a plau-sible model of galaxies in our Universe, then to assign each haloin the simulation a weight based on how closely its enclosed massat some radius (we typically use 50 kpc in what follows) matchesobservational data. The resulting weights for the halo sample thenprovide a constraint on the enclosed mass of the MW at other radii.In more detail, we take an observational measurement of thetotal MW mass within a specific radius and assign a weight to eachhalo in the Illustris galaxy catalog: assuming the observed masshas a value of µ and an associated (Gaussian) error of σ , then theweight W i contributed by an individual halo i with enclosed mass M i at the specified radius is W i = 1 √ πσ exp (cid:18) − ( M i − µ ) σ (cid:19) . (1)We can then construct the full mass or circular velocity profile andcompute the total stellar or halo mass that is consistent with theobserved constraint by using the distribution of weights assignedto the haloes. In this analysis, we assume that observed constraintsall follow Gaussian distributions, consistent with the analyses weincorporate, but this technique can be easily extended to any otheranalytic or numerical probability distribution. In what follows, wequote median values and confidence intervals that are centred onthe median and contain 68% of the probability distribution.The primary observational constraint we use is the total massof the MW within 50 kpc, M ( <
50 kpc) . There are many litera-ture estimates of the MW’s mass at approximately this scale (e.g.,Wilkinson & Evans 1999; Battaglia et al. 2005; Xue et al. 2008;Brown et al. 2010; Gnedin et al. 2010; Kafle et al. 2014; Eadieet al. 2015), in large part because (1) this is approximately thedistance to which large samples of blue horizontal branch (BHB)stars are currently available from surveys such as the Sloan Dig-ital Sky Survey, and (2) the LMC lies at a Galactocentric dis-tance of ≈
50 kpc , meaning estimates of the MW mass based onLMC’s dynamics directly constrain M ( <
50 kpc) . We focus ontwo recent and disparate measurements of M ( <
50 kpc) : Dea-son et al. (2012, hereafter D12), who used BHB stars and found M ( <
50 kpc) = 4 . ± . × M (cid:12) , and Gibbons et al.(2014, hereafter G14), who used the Sagittarius stream to measure M ( <
50 kpc) = 2 . ± . × M (cid:12) (G´omez et al. 2015). Thesemeasurements are clearly incompatible at the σ level and there-fore are useful for showing the effects of varying M ( <
50 kpc) on the inferred mass distribution at larger radii. In Sec. 3.2, weexplicitly show how estimates of M , c vary as a function of M ( <
50 kpc) .In principle, a complete analysis would include every darkmatter halo in Illustris. In practice, however, only a relatively nar-row range of masses contribute any weight to our inferences. We Here and throughout this work, we use ’mass’ to refer to the enclosedmass (as opposed to mass within a spherical shell) therefore restrict our analysis to all haloes with M , c = (0 . − × M (cid:12) , which includes
14 192 haloes for Illustris-1,
14 316 haloes for Illustris-2, and
12 885 haloes for Illustris-3. As we showbelow, this mass range is more than sufficient for including all rel-evant haloes in our analysis and does not bias our results in anyway.
Fig. 1 presents the mass distributions obtained using the constraintson M ( <
50 kpc) from D12 (left-hand panel) and G14 (right-handpanel) from the Illustris-1 sample. The best-fitting Navarro-Frenk-White (1997, hereafter, NFW) profiles for the total mass distribu-tion are given in the figure as well. The fits were performed over theradial range of 40–300 kpc, as we find a lack of convergence amongdifferent resolution versions of Illustris on smaller scales (see be-low; convergence in density profiles should occur at smaller scales,as density is a differential quantity while mass and circular velocityare cumulative quantities). Unsurprisingly, given the significantlyhigher value of M ( <
50 kpc) found in D12 relative to G14, thebest-fitting NFW value of M , c for D12 is much larger than forG14, . × M (cid:12) versus . × M (cid:12) . The best fitting con-centration parameters are similar: c , c = 12 . ± . for D12and c , c = 13 . ± . for G14. Both of these concentrationsare larger than those derived from large DMO simulations, whichtypically find c , c ≈ . for haloes of M , c ≈ M (cid:12) (e.g.,Dutton & Macci`o 2014).The lower panels of the figures show the fractional differencesof Illustris-2 and Illustris-3 with respect to their high-resolutioncounterpart, with error bars representing 68% confidence intervals.There are relatively large differences between the different levels ofresolution at relatively small radii ( r < kpc), while differencesare much less substantial farther away from Galactic Centre. Witha gravitational softening length ∼ ≈ kpc (their fig. 3).We therefore strongly caution against extrapolating the NFW fitspresented in this paper to small radii ( r (cid:46) kpc). If future gen-erations of simulations provide well-converged results at smallerradii, the dark matter fraction within ∼ disc scale lengths willlikely provide important constraints on feedback models (Courteau& Dutton 2015).The circular velocity profiles, V circ ( r ) , corresponding to thecumulative mass profiles of Fig. 1 are shown in the left-hand panelof Fig. 2. This highlights the large difference in the two determina-tions of the MW potential, as well as how this difference persists inpredicted profiles out to 300 kpc. It is only at distances > kpcthat the 68% confidence intervals begin to overlap.The distribution of mass among dark matter, stars, and gaswithin any given radius is interesting to consider: observationally,we can measure the stellar mass with reasonable accuracy and in-fer the dark matter mass, but constraining the distribution of theGalaxy’s gaseous component at large distances is much more diffi-cult (see, e.g., Gupta et al. 2012; Fang et al. 2013). In the right-handpanel of Fig. 2, we plot the circular velocity profile decomposedinto the contributions from each of these components. Dark mat-ter dominates the potential at all radii we study, and while stars MNRAS461
50 kpc) found in D12 relative to G14, thebest-fitting NFW value of M , c for D12 is much larger than forG14, . × M (cid:12) versus . × M (cid:12) . The best fitting con-centration parameters are similar: c , c = 12 . ± . for D12and c , c = 13 . ± . for G14. Both of these concentrationsare larger than those derived from large DMO simulations, whichtypically find c , c ≈ . for haloes of M , c ≈ M (cid:12) (e.g.,Dutton & Macci`o 2014).The lower panels of the figures show the fractional differencesof Illustris-2 and Illustris-3 with respect to their high-resolutioncounterpart, with error bars representing 68% confidence intervals.There are relatively large differences between the different levels ofresolution at relatively small radii ( r < kpc), while differencesare much less substantial farther away from Galactic Centre. Witha gravitational softening length ∼ ≈ kpc (their fig. 3).We therefore strongly caution against extrapolating the NFW fitspresented in this paper to small radii ( r (cid:46) kpc). If future gen-erations of simulations provide well-converged results at smallerradii, the dark matter fraction within ∼ disc scale lengths willlikely provide important constraints on feedback models (Courteau& Dutton 2015).The circular velocity profiles, V circ ( r ) , corresponding to thecumulative mass profiles of Fig. 1 are shown in the left-hand panelof Fig. 2. This highlights the large difference in the two determina-tions of the MW potential, as well as how this difference persists inpredicted profiles out to 300 kpc. It is only at distances > kpcthat the 68% confidence intervals begin to overlap.The distribution of mass among dark matter, stars, and gaswithin any given radius is interesting to consider: observationally,we can measure the stellar mass with reasonable accuracy and in-fer the dark matter mass, but constraining the distribution of theGalaxy’s gaseous component at large distances is much more diffi-cult (see, e.g., Gupta et al. 2012; Fang et al. 2013). In the right-handpanel of Fig. 2, we plot the circular velocity profile decomposedinto the contributions from each of these components. Dark mat-ter dominates the potential at all radii we study, and while stars MNRAS461 , 3483–3493 (2016)
C. Taylor et al. M ( < r ) [ M fl ] NFW Fit:c = 12.2 ± M , c = 1.11 ± × M fl D12 r [kpc] - ( M i / M ) Illustris-2Illustris-3 10 M ( < r ) [ M fl ] NFW Fit:c = 13.2 ± M , c = 0.61 ± × M fl G14 r [kpc] - ( M i / M ) Illustris-2Illustris-3
Figure 1.
Top: the mass distribution M ( < r ) derived from Illustris-1, using the D12 (left) and G14 (right) constraints on M ( <
50 kpc) . Error bars represent68% confidence intervals. The grey lines show best-fitting NFW profiles for the full mass distribution (dark matter and baryonic); the NFW fit parameters aregiven in the figure. The two constraints result in very different estimates of M , c ; the concentration parameters are less disparate. Bottom: residuals betweenthe mass distribution obtained from Illustris-1 and Illustris-2 (black circles) or Illustris-3 (grey triangles). At r < kpc, systematic differences are evident;these likely result from a combination of numerical resolution and differences in the stellar masses at fixed halo mass. Small differences of 2-5% exist at largerradii; however, such deviations are much smaller than the uncertainties we derive in Table 1. r [kpc] V c i r c ( r ) [ k m / s ] D12G14 r [kpc] V c i r c ( r ) [ k m / s ] Dark MatterGasStars
Figure 2.
Circular velocity curves. Left: V circ ( r ) for the mass profiles given in Fig. 1. The overall mass profile using G14’s constraint is lower at everyGalactocentric distance compared to the profile derived using the D12 constraint, and the 68% confidence intervals are disjoint up to
280 kpc . Right: adecomposition of the circular velocity profile derived using the D12 constraint (black points in the left panel of Fig. 1) into separate contributions from darkmatter (black), stars (blue) and gas (red). At all radii probed here, dark matter dominates. The contribution from gas matches that from stars near a halo-centricdistance of 100 kpc. MNRAS , 3483–3493 (2016) ass of the Milky Way from Illustris log [ M (< 100 kpc)/ M fl ] P r ob a b ilit y Figure 3.
The probability distribution of M ( <
100 kpc) derived from theG14 (squares, connected by dashed lines) and D12 (circles, connected bysolid lines) constraints on M ( <
50 kpc) . The colours represent the indi-vidual resolution levels: Illustris-1 (black), Illustris-2 (blue), and Illustris-3(red). The excellent agreement across the three levels of resolution indicatesthat the total mass profiles are well-converged in Illustris. substantially outweigh the gas for r (cid:46) kpc, the two contributeapproximately the same mass by r ≈ kpc. In this subsection, we explore predictions for enclosed masses atspecific radii in more detail. In particular, we are interested in un-derstanding how observational constraints at M ( <
50 kpc) trans-late into inferences on masses at other radii. We consider both in-dividual physical radii (in particular, 100 and 250 kpc) and vari-ous definitions of spherical overdensity masses ( M , c , M vir , and M , m ).Fig. 3 presents the probability distribution for M ( <
100 kpc) ,with black, blue, and red symbols representing Illustris-1, Illustris-2, and Illustris-3 respectively. The results using the D12 constrainton M ( <
50 kpc) are presented as circles connected with solidlines, while those using the G14 constraint are shown as squareswith dashed connecting lines. As expected, and shown previously,the D12 constraint on M ( <
50 kpc) results in a significantly higherpredicted total mass within 100 kpc (approximately 0.2 dex). Thesmaller (relative) error quoted in D12 also results in a narrower dis-tribution for M ( <
100 kpc) .Perhaps the most important aspect of Fig. 3 is the excellent convergence seen across the three Illustris simulations (a factor of64 in mass resolution and 4 in force resolution). Not only is thepeak or median value well converged, the entire distribution is es-sentially identical in each case. This indicates that, while masseson small scales (10–30 kpc) are affected by resolution and bary-onic physics, enclosed masses at larger radii are not subject to sucheffects. The consistency of the mass distributions at large radii, sub-ject to a constraint at 50 kpc, points to robustness of our techniquefor constraining the mass distribution of the MW.Inferred values of aperture masses within 100 and 250 kpc andthree different spherical overdensity masses, along with 68% and
Table 1.
Median values, along with 68% and 90% confidence intervals,for mass measures explored in this paper; all masses are expressed in unitsof M (cid:12) . In each case, we calculate values using constraints from bothD12 (column 2) and G14 (column 3) on each of the three Illustris resolutionlevels. Good convergence across the three levels of resolution is evident. D12 G14Illustris-1 M , c . +0 .
370 (0 . − .
240 (0 . . +0 .
196 (0 . − .
148 (0 . M vir . +0 .
511 (1 . − .
304 (0 . . +0 .
251 (0 . − .
179 (0 . M , m . +0 .
642 (1 . − .
361 (0 . . +0 .
306 (0 . − .
213 (0 . M ( <
100 kpc) 0 . +0 .
091 (0 . − .
090 (0 . . +0 .
087 (0 . − .
076 (0 . M ( <
250 kpc) 1 . +0 .
334 (0 . − .
236 (0 . . +0 .
209 (0 . − .
164 (0 . Illustris-2 M , c . +0 .
296 (0 . − .
196 (0 . . +0 .
166 (0 . − .
134 (0 . M vir . +0 .
419 (0 . − .
243 (0 . . +0 .
213 (0 . − .
163 (0 . M , m . +0 .
490 (1 . − .
306 (0 . . +0 .
253 (0 . − .
186 (0 . M ( <
100 kpc) 0 . +0 .
090 (0 . − .
078 (0 . . +0 .
077 (0 . − .
077 (0 . M ( <
250 kpc) 1 . +0 .
283 (0 . − .
194 (0 . . +0 .
184 (0 . − .
148 (0 . Illustris-3 M , c . +0 .
332 (0 . − .
193 (0 . . +0 .
170 (0 . − .
135 (0 . M vir . +0 .
441 (0 . − .
249 (0 . . +0 .
220 (0 . − .
160 (0 . M , m . +0 .
562 (1 . − .
293 (0 . . +0 .
256 (0 . − .
190 (0 . M ( <
100 kpc) 0 . +0 .
092 (0 . − .
080 (0 . . +0 .
078 (0 . − .
073 (0 . M ( <
250 kpc) 1 . +0 .
300 (0 . − .
192 (0 . . +0 .
178 (0 . − .
148 (0 .
90% confidence intervals, are given in Table 1. The estimated virialmass, M vir , using D12 is . × M (cid:12) , with a 90% confidenceinterval of . − . × M (cid:12) . This is similar to the result ofBoylan-Kolchin et al. (2013), who found a 90% confidence intervalof . − . × M (cid:12) for M vir based on the dynamics of theLeo I satellite galaxy. Using the G14 estimate of M ( <
50 kpc) ,we find a median value of M vir = 0 . × M (cid:12) with a 90%confidence interval of . − . × M (cid:12) , both of which aresubstantially lower than our inference based on the results of D12.These results highlight the importance of accurate determinationsof M ( <
50 kpc) for understanding the large-scale properties ofthe MW. We note that the 99.95% confidence interval for haloesconsistent with the D12 constraint is . × < M , c < . × M (cid:12) (the range for the G14 constraint is . × 50 kpc) , it is also importantto understand how inferences of spherical overdensity masses de-pend on M ( < 50 kpc) . To do this, we assume that M ( < 50 kpc) can be measured with an accuracy of 10% (i.e., X ± . X ) andcompute the resulting median value and 68% confidence intervalsfor M , c . The resulting dependence of M , c on M ( < 50 kpc) is shown in Fig. 4, where the error bars show 68% confidence in- MNRAS461 MNRAS461 , 3483–3493 (2016) C. Taylor et al. log [ M (< 50 kpc) / M fl ] l og [ M , c / M fl ] G D Figure 4. The dependence of the inferred value of M , c on the input(measured) value of M ( < 50 kpc) . The data points with error bars showvalues of M , c based on our weighting procedure, assuming a 10% errorin M ( < 50 kpc) . The G14 and D12 determinations of M ( < 50 kpc) arehighlighted with yellow and magenta vertical bands, respectively, with thewidths of the bands showing the 68% confidence intervals. The best-fittinglog-quadratic relation between M ( < 50 kpc) and M , c (given in equa-tion 2) is plotted as a solid grey line, while the dashed grey line shows thefit to the unweighted data; see the text for details. This relation can be usedto map any constraint on M ( < 50 kpc) to an inferred value of M , c . tervals. It is clear that there is a strong correlation between M ( < 50 kpc) and M , c . We fit this with a quadratic function in logspace: log (cid:18) M , c M (cid:12) (cid:19) = A + B µ + C µ , (2) µ = log (cid:18) M ( < 50 kpc)4 × M (cid:12) (cid:19) . Fitting to the weighted results plotted in Fig. 4, we find A =12 . , B = 1 . , C = 0 . with an rms scatter of 0.069, whereasfitting the unweighted data, we find A = 12 . , B = 1 . , C =0 . with an rms scatter of 0.067. The latter is offset slightlyhigher at fixed M ( < 50 kpc) , as the weighted results naturallyinvolve averaging over the dark halo mass function within eachbin, which is a steeply declining function of mass, whereas the un-weighted results do not.Equation 2 can be used to convert any constraint on M ( < 50 kpc) into a constraint on M , c . It is also straightforward toconvert this fit to a constraint on M vir or M , m , as M vir ≈ . M , c and M , m ≈ . M , c for the typical mass pro-files in Illustris. If Equation 2 or a similar relation holds broadlyfor other hydrodynamic simulations with different galaxy forma-tion physics implementations, then it will be of tremendous valuefor MW mass inference studies. We plan to examine this issue inmore detail in future work (and see further discussion below). Our primary analysis, presented over the previous subsections,makes use of the highest resolution Illustris simulation. This, and all other hydrodynamic simulations of the evolution of a represen-tative galaxy population over cosmic time, require a number of as-sumptions in order to produce a realistic set of galaxies. One ofthe primary calibrations for Illustris, for example, was to match the z = 0 galaxy stellar mass function. As shown in fig. 7 of Vogels-berger et al. (2014a), the galaxy formation prescriptions in Illustrisresult in notable changes in the total masses of dark matter haloesover a wide range in halo mass. Moreover, these changes dependon specific choices made in the galaxy formation modelling, as thegalaxy formation modelling within the Eagle simulation results insubstantially different effects on halo masses (see fig. 1 of Schalleret al. 2015).It is not a priori obvious whether using the DMO run shouldresult in similar or different predictions from the fully hydrody-namic simulation, and if the results are different, it is not clearwhether they will be higher or lower. Certainly, we expect that theformation of a galaxy will lead to a more centrally concentratedmass distribution relative to the DMO run, to some extent. Adi-abatic contraction of the dark matter in response to gas coolingwill also tend to increase the amount of dark matter in the cen-tral regions of the halo. On the other hand, it is well establishedthat galaxy formation must be inefficient in Λ CDM (e.g., Fukugita& Peebles 2004), meaning that only a relatively small fraction ofthe baryonic allotment of a dark matter halo ( ∼ for MW-mass haloes) will be converted into stars by z = 0 . Strong feed-back from galaxy formation can change the structure of dark mat-ter haloes, reducing their mass within a given radius compared towhat would be obtained in a DMO version (e.g., Vogelsberger et al.2014b; Schaller et al. 2015). It is therefore of great interest to studyprecisely how inferences about the MW’s mass profile change fromusing DMO simulations – which, for given cosmological parame-ters, are uniquely predicted – to using cosmological hydrodynamicsimulations.The first test we perform to gauge the effects of includinggalaxy formation physics on mass inferences is to rerun our anal-ysis on the DMO versions of Illustris. Fig. 5 shows the resultsof applying the D12 constraint to Illustris-Dark-1. It can be di-rectly compared to Fig. 1, in which the D12 constraint was ap-plied to the hydrodynamic version of Illustris-1. Relative to thefull Illustris simulation, inferences based on the DMO version re-sult in a significantly higher estimate of M , c ( . × ver-sus . × M (cid:12) ) and a significantly lower version of the NFWconcentration ( c = 7 . versus 12.3). Table 2 provides an alternateversion of Table 1 in which all constraints are obtained using theDMO version of Illustris-1. In all cases, the net effect of using theDMO run rather than the hydrodynamic version is to infer higher values for a given aperture mass.We can use the Illustris suite to perform an additional test ofthe effects of galaxy formation on the mass distribution within darkmatter halos (and for accompanying inferences on the mass distri-bution of the MW): since Illustris and Illustris-Dark share the sameinitial conditions, individual dark matter halos can be matched be-tween the two simulations (for details, see section 3.2 of Vogels-berger et al. 2014a). In this way, we can study the effects of galaxyformation on a halo-by halo basis by identifying the DMO analogueof each halo in the full Illustris run and comparing the resultingmass distributions.Fig. 6 shows the results of this comparison, for which weuse haloes in Illustris-1 falling within the 68% confidence inter-val of M , c computed using the D12 constraint (see Table 1) –assigning equal weight to all such haloes – and their counterpartsin the DMO run. The left-hand panel shows how the density pro- MNRAS , 3483–3493 (2016) ass of the Milky Way from Illustris r [kpc] M ( < r ) [ M fl ] Illustris-Dark-1 NFW Fit:c = 7.32 ± M , c = 1.55 ± × M fl Figure 5. The cumulative mass distribution from Illustris-Dark-1 (blackpoints with error bars), along with the best-fitting NFW profile (grey line),derived assuming the D12 constraint on M ( < 50 kpc) . This figure can bedirectly compared to the left-hand panel of Fig. 1, which shows the samequantities from the full hydrodynamic run. While both versions of Illus-tris are well-fitted by NFW profiles, the fit parameters differ substantiallybetween the two: the DMO run is fitted by a higher-mass (37% higher),lower-concentration (40% lower) halo. If DMO runs are used for modellingthe MW mass distribution based on M ( < 50 kpc) , or a similar constraint,this effect must be taken into account. files are affected at each radius. On small scales ( r (cid:46) 30 kpc ),the hydrodynamic run has higher densities on a halo-by-halo basis.This is caused by the formation of the central galaxy, both throughits mass and through any adiabatic contraction. On larger scales( r (cid:38) 40 kpc ), a given halo in the hydrodynamic run is less densethan its equivalent in the DMO run by approximately 20%. Thisreduction in density is likely caused by outflows and the loss ofgas mass (or the prevention of gas accretion). The effect on thecumulative mass distribution is shown in the right-hand panel ofFig. 6. On a halo-by-halo basis, the hydrodynamic run results inlarger masses out to ≈ kpc; on larger scales, the masses inthe DMO run are larger, with the difference reaching an asymptoticvalue of ≈ at 250-300 kpc. As discussed in Section 4, the de-tails of the reduction in mass may depend on the adopted modelsof galaxy formation modelling. We can also use the technique explored in the previous sections tocompute the galaxy stellar masses from Illustris that are consistentwith the adopted mass constraints at 50 kpc. Table 3 gives themedian values as well as 68% and 90% confidence intervals basedon the D12 and G14 constraints in each of the three Illustris reso-lution levels. Unlike the total enclosed mass at large radii, whichis well-converged across the three different Illustris resolutions,the stellar masses in these haloes increase by a factor of ∼ fromIllustris-3 to Illustris-1. This difference is not large enough to bereflected in stellar mass functions (which are reasonably similar forthe different resolution levels studied here; see, e.g., Vogelsbergeret al. 2013 and Torrey et al. 2014). It is larger than the uncertainty Table 2. Median values, along with 68% and 90% confidence intervals,for a variety of mass measures explored in this paper (similar to Table 1);all values are in units of M (cid:12) . In contrast to Table 1, however, theseresults use the Illustris-Dark simulations. D12 G14Illustris-Dark-1 M , c . +0 . 460 (1 . − . 343 (0 . . +0 . 296 (0 . − . 220 (0 . M vir . +0 . 676 (1 . − . 445 (0 . . +0 . 382 (0 . − . 275 (0 . M , m . +0 . 841 (1 . − . 545 (0 . . +0 . 466 (0 . − . 318 (0 . M ( < 100 kpc) 0 . +0 . 116 (0 . − . 095 (0 . . +0 . 105 (0 . − . 095 (0 . M ( < 250 kpc) 1 . +0 . 357 (0 . − . 299 (0 . . +0 . 283 (0 . − . 231 (0 . Illustris-Dark-2 M , c . +0 . 515 (1 . − . 344 (0 . . +0 . 296 (0 . − . 219 (0 . M vir . +0 . 712 (1 . − . 452 (0 . . +0 . 399 (0 . − . 268 (0 . M , m . +0 . 920 (2 . − . 565 (0 . . +0 . 476 (1 . − . 312 (0 . M ( < 100 kpc) 0 . +0 . 116 (0 . − . 095 (0 . . +0 . 103 (0 . − . 094 (0 . M ( < 250 kpc) 1 . +0 . 396 (0 . − . 290 (0 . . +0 . 290 (0 . − . 221 (0 . Illustris-Dark-3 M , c . +0 . 567 (1 . − . 349 (0 . . +0 . 302 (0 . − . 231 (0 . M vir . +0 . 760 (1 . − . 488 (0 . . +0 . 389 (0 . − . 271 (0 . M , m . +1 . 00 (2 . − . 601 (0 . . +0 . 459 (1 . − . 316 (0 . M ( < 100 kpc) 0 . +0 . 115 (0 . − . 113 (0 . . +0 . 107 (0 . − . 090 (0 . M ( < 250 kpc) 1 . +0 . 417 (0 . − . 311 (0 . . +0 . 292 (0 . − . 225 (0 . Table 3. Inferred values of M (cid:63) , in units of M (cid:12) , using the D12 (col-umn 2) and G14 (column 3) constraints on M ( < 50 kpc) . The quoted er-rors are the 68% and 90% confidence intervals. The Illustris feedback pre-scriptions were calibrated for the highest-resolution simulation (Illustris-1),so perfect convergence in M (cid:63) across the three simulations is not expected.D12 G14Illustris-1 . +1 . 47 (2 . − . 32 (2 . . +0 . 98 (1 . − . 72 (1 . Illustris-2 . +1 . 19 (2 . − . 15 (1 . . +0 . 76 (1 . − . 59 (0 . Illustris-3 . +0 . 79 (1 . − . 71 (1 . . +0 . 48 (0 . − . 35 (0 . on the measured M (cid:63) of the MW, however: most recent estimatesfor the Galaxy fall in the range M (cid:63) = 5 − . × M (cid:12) (e.g.,McMillan 2011; Bovy & Rix 2013; Licquia & Newman 2015).Differences in the simulated stellar masses at the factor of ∼ level are unsurprising, as the galaxy formation models usedin the Illustris suite were calibrated at the resolution of Illustris-1;we would not expect the same models to work identically at signif-icantly lower resolution. Specifically, the minimum resolution re-quired for the feedback implementation in Illustris to produce a re-alistic galaxy population is not achieved in Illustris-3 (Vogelsbergeret al. 2013). We therefore consider the results from Illustris-1 to bethe most reasonable comparison to make with observations.We adopt the measurement of Licquia & Newman (2015, MNRAS461 35 (0 . on the measured M (cid:63) of the MW, however: most recent estimatesfor the Galaxy fall in the range M (cid:63) = 5 − . × M (cid:12) (e.g.,McMillan 2011; Bovy & Rix 2013; Licquia & Newman 2015).Differences in the simulated stellar masses at the factor of ∼ level are unsurprising, as the galaxy formation models usedin the Illustris suite were calibrated at the resolution of Illustris-1;we would not expect the same models to work identically at signif-icantly lower resolution. Specifically, the minimum resolution re-quired for the feedback implementation in Illustris to produce a re-alistic galaxy population is not achieved in Illustris-3 (Vogelsbergeret al. 2013). We therefore consider the results from Illustris-1 to bethe most reasonable comparison to make with observations.We adopt the measurement of Licquia & Newman (2015, MNRAS461 , 3483–3493 (2016) C. Taylor et al. r [kpc] - ρ D M O ( r ) / ρ ( r ) r [kpc] - M D M O ( < r ) / M ( < r ) Figure 6. Fractional differences in the density (left) and enclosed mass (right) profiles between Illustris-1 and Illustris-Dark-1, where haloes are individuallymatched across the two simulations (see the text for details). Data points represent the median differences between Illustris-1 and Illustris-Dark-1, while errorbars show the central 68% range of the data. On small scales, the inclusion of baryonic physics results in more mass at a given radius owing to the formationof the central galaxy. On large scales, however, feedback causes an overall reduction in mass on a halo-by-halo basis for the full hydrodynamic simulationrelative to the DMO run. The effect in the density profile is ∼ at large radii, while the effect in the cumulative mass profile is ∼ at large distances. hereafter, LN15), in which the authors used results derived in Bovy& Rix (2013) to obtain M (cid:63) = 6 . ± . × M (cid:12) , as a rep-resentative value of the stellar mass of the MW and use it as a ref-erence point in what follows. Comparing this number to the resultsfor Illustris-1 in Table 3, we see that D12 agrees well with the ob-served value, while G14 is substantially lower. This is not surpris-ing, given the results of Table 1. The very low value of M , c ob-tained based on G14 is much lower than the typical value found forhaloes with the stellar mass of the MW via either abundance match-ing (Guo et al. 2010; Behroozi et al. 2013; Moster et al. 2013),galaxy-galaxy lensing (Mandelbaum et al. 2016), or satellite kine-matics (e.g., Watkins et al. 2010; Boylan-Kolchin et al. 2013). Evenaccounting for possible differences in halo masses of red and bluegalaxies at fixed stellar mass (Mandelbaum et al. 2016), the MWwould be a strong outlier if its mass is as low as the median value in-dicated by our analysis using the G14 constraint on M ( < 50 kpc) .As noted in Section 2, our methodology for constraining themass profile of the MW is quite general. While we have focused onconstraining the total mass at large radii based on measurements ofthe total mass within 50 kpc, we can instead use other quantities –for example, M (cid:63) – for our inference. Following the same procedureoutlined in Section 2, we estimate M ( < 50 kpc) , M ( < 100 kpc) ,and the three spherical overdensity masses used above based onLN15’s determination of M (cid:63) ; the results are presented in the sec-ond column of Table 4. The results are very similar to those ob-tained using the D12 determination of M ( < 50 kpc) , with LN15-based estimates being 5-7% higher (the results are approximately afactor of 1.6–2 larger than G14-based estimates).Since M (cid:63) and M ( < 50 kpc) can be considered independentvariables, we can also study the joint probability of obtaining var-ious mass measures conditioned on M (cid:63) and M ( < 50 kpc) . Thesejoint constraints, using D12’s estimate of M ( < 50 kpc) , are givenin the third column of Table 4. The joint constraints are simi-lar to both the estimates using M (cid:63) alone and the estimate using M ( < 50 kpc) (from D12) alone, which is a result of the good M (< 50 kpc) [ M fl ] M [ M fl ] G D LN 15 l og [ M , c / M fl ] Figure 7. The correlation between M (cid:63) and M ( < 50 kpc) for all haloes inthe Illustris-1 sample; the haloes are coloured by the value of M , c . Ver-tical shaded bands show G14 (yellow) and D12 (magenta) determinationsof M ( < 50 kpc) while the horizontal band shows the LN15 determinationof M (cid:63) for the MW. Very few haloes agree with both the G14 measurementof the total mass at 50 kpc and the MW’s stellar mass; many more of thesimulated galaxies match the D12 value for M ( < 50 kpc) and the LN14 M (cid:63) value simultaneously. agreement of each of these estimates individually. Had we used theG14 value of M ( < 50 kpc) , the constraints would have shifted sub-stantially. This is highlighted in Fig. 7, which shows the Illustris-1data in M (cid:63) − M ( < 50 kpc) space; each halo assigned a colour MNRAS , 3483–3493 (2016) ass of the Milky Way from Illustris Table 4. The 68% and 90% confidence intervals of various mass measuresof the MW (all in units of M (cid:12) ) inferred using the LN15 measurementof the MW’s M (cid:63) alone (column 2) and jointly with the D12 constraint on M ( < 50 kpc) (column 3). P ( M | M (cid:63) ) P ( M | M (cid:63) , M D12 ) M , c . +0 . 62 (1 . − . 35 (0 . . +0 . 36 (0 . − . 22 (0 . M vir . +0 . 81 (1 . − . 43 (0 . . +0 . 49 (1 . − . 29 (0 . M , m . +1 . 00 (2 . − . 53 (0 . . +0 . 65 (1 . − . 35 (0 . M ( < 100 kpc) 0 . +0 . 179 (0 . − . 145 (0 . . +0 . 081 (0 . − . 087 (0 . M ( < 250 kpc) 1 . +0 . 51 (1 . − . 34 (0 . . +0 . 34 (0 . − . 21 (0 . Fractional Error on M (< r ) F r ac ti on a l E rr o r on M , c r = 50 kpc r = 80 kpc r = 100 kpc Figure 8. The precision attained in measurements of M , c as a functionof the precision in the input constraint. We consider input constraints of thetotal mass within 50, 80, and 100 kpc (black circles, red squares, and bluetriangles, respectively) to show how more precise determinations of masseswithin larger radii can affect the inferred value of M , c . The figure showsthe trade-off between precision and distance: for example, an error of 15%in M ( < 100 kpc) results in the same precision in the estimate of M , c as an error of 9% in M ( < 50 kpc) . according to its value of M , c . The intersection of the D12 andLN15 constraints falls along the main locus of the points while theG14 constraint intersects the LN 15 constraint in a part of parame-ter space with very few haloes. The D12 and LN15 measurementsare therefore in good agreement based on the Illustris haloes, whilethe G14 and LN15 measurements are not. As larger samples of halo stars at greater distances become avail-able, it may become possible to constrain the mass of the MW en-closed within 80 or even 100 kpc (see, e.g., Gnedin et al. 2010;Cohen et al. 2015 for initial work in this direction). Such measure-ments would have the benefit of providing stronger constraints onthe virial mass of the MW. Fig. 8 shows the fractional uncertaintyin M , c as a function of the error in the mass contained within 50 (black circles), 80 (red squares), and 100 kpc (blue triangles). Ata fixed uncertainty in M ( < r ) , the implied uncertainty in M , c does indeed become smaller as one moves to greater Galactocentricdistance.The figure quantifies how improving uncertainties at a givendistance will be reflected in uncertainties on M , c : for example,reducing the error on M ( < 50 kpc) from 10 to 5% would reducethe error on M , c from 28 to 23%. On the other hand, an mea-surement of the mass within 80 kpc that is accurate to 10% resultsin an error of 22% in M , c , while the same accuracy on a mea-surement of the mass within 100 kpc of the Galaxy would yielderrors of 19% in M , c . The figure also shows the fundamentallimitations in extrapolating to M , c based on measured aperturemasses within smaller radii. Some level of irreducible uncertaintyis unavoidable in standard cosmological models, as extrapolationfrom mass at a given radius to the virial radius depends on the haloconcentration (e.g. Navarro et al. 1997; Bullock et al. 2001). Forexample, consider the recent study of Williams & Evans (2015),who found that M ( < 50 kpc) = 4 . +0 . − . × M (cid:12) , or an er-ror of approximately 3% on M ( < 50 kpc) . Using this constraint,we obtain M , c = 1 . +0 . − . × M (cid:12) ; the uncertainty on thederived value of M , c remains large in spite of the high precisionof the input measurement. Fig. 8 makes it clear that measurementsof the mass within 50 (80, 100) kpc will result in an uncertainty on M , c of no better than 23% (17%, 14%).A central assumption of the techniques we employ here is thatIllustris-1 provides a faithful representation of galaxies and the ef-fects of galaxy formation on dark matter halo structure. Since cos-mological hydrodynamic simulations are still at the point of relyingon subgrid models of physics, and will be for the foreseeable future,a logical extension of our work would be to investigate predictionsin future generations of simulations to test the robustness of our re-sults. It would also be interesting to compare the results we haveobtained with Illustris to the Eagle simulations, as the galaxy for-mation modelling employed there is somewhat different. Given thedifferences seen in the ratio of masses in hydrodynamic to DMOsimulations in Illustris versus Eagle (compare fig. 7 of Vogelsbergeret al. 2014a and fig. 1 of Schaller et al. 2015), such a comparisonwould be timely.One effect that appears to be particularly important for set-ting the amount of mass reduction for a given halo in the hydro-dynamic run relative to its counterpart in the DMO version is theunderlying model of AGN feedback. Vogelsberger et al. (2014b)adopted an AGN model that drives very strong outflows, perhapsunrealistically so (Genel et al. 2014). Forthcoming updates to theIllustris suite will use modified versions of AGN feedback that areless powerful and may result in different modifications of the large-scale halo properties of galaxies, which may in turn affect how M ( < 50 kpc) maps on to M , c .To explore the potential impact of this effect on our results,we use the current generation of Illustris and compare the effects ofBH mass for galaxies of a fixed halo mass (we use the haloes thatare closest to the median value of M , c found in Illustris-Dark-1using the D12 constraint). We rank this sample according to blackhole mass and then compute the difference in mass in the hydrody-namic simulation relative to the DMO run. There is indeed a differ-ence: the galaxies with the highest-mass black holes show a 20%reduction in their overall mass, on average, while the galaxies withthe lowest mass black holes see a 10% reduction in mass comparedto their DMO counterparts. The total halo mass therefore appearsto depend somewhat on the choice of black hole feedback model,although this does not appear to be a large source of uncertainty MNRAS461 50 kpc) maps on to M , c .To explore the potential impact of this effect on our results,we use the current generation of Illustris and compare the effects ofBH mass for galaxies of a fixed halo mass (we use the haloes thatare closest to the median value of M , c found in Illustris-Dark-1using the D12 constraint). We rank this sample according to blackhole mass and then compute the difference in mass in the hydrody-namic simulation relative to the DMO run. There is indeed a differ-ence: the galaxies with the highest-mass black holes show a 20%reduction in their overall mass, on average, while the galaxies withthe lowest mass black holes see a 10% reduction in mass comparedto their DMO counterparts. The total halo mass therefore appearsto depend somewhat on the choice of black hole feedback model,although this does not appear to be a large source of uncertainty MNRAS461 , 3483–3493 (2016) C. Taylor et al. in our predictions. Future generations of Illustris-like simulationswith modified black hole feedback models will allow us to directlytest the effects on inferences regarding the MW mass.It is not entirely obvious how the effects of vigorous feedbackpropagate through our analysis, as this will depend on the changein enclosed mass within 50 kpc relative to the change in enclosedmass within larger radii. However, given that the black hole feed-back in the current version of Illustris may be too effective and thatthe larger-mass black holes correlate with larger reductions in halomass as compared to lower-mass black holes, it is likely that anymodified prescriptions will result in slightly higher inferences onthe total halo mass compared to our current results, should there bea difference.Future work would also benefit significantly from cosmolog-ical simulations with larger volumes and higher mass resolution.Importance sampling relies on having a well-sampled parameterspace, which can be an issue if not many haloes match the desiredconstraint(s) (see Busha et al. 2011 and Gonz´alez et al. 2014 formore details). Our current analysis has many haloes contributingsignificant weights: 870 and 2196 haloes contribute weights that areat least 10% of the maximum possible weight ( W max = 1 / √ πσ from equation. 1) for the D12 and G14 constraints, respectively.However, if we wish to add additional restrictions – based on mor-phology, disc size, star formation history, or specific star formationrate, for example – the sample would likely become significantlysmaller, which would be the limiting factor in the conclusions wecould draw. With larger sample sizes, such concerns would beeliminated. From Fig. 7, joint constraints on M ( < 50 kpc) and M (cid:63) are unlikely to be strongly affected by sample size unless amuch larger volume produced many haloes with much larger stel-lar masses at fixed halo mass [in which case, the G14 measurementof M ( < 50 kpc) would be more consistent with the simulation re-sults than it is at present]. In this paper, we have explored how the Illustris suite can be used toinform our understanding of the mass distribution around the MW.Our main conclusions are as follows. • The mass profiles of haloes consistent with a given constrainton M ( < 50 kpc) differ substantially between DMO and hydrody-namic versions of Illustris. Using DMO simulations to extrapolatefrom 50 kpc to larger radii results in an overestimate of the halomass and an underestimate of the halo concentration. • The effects of baryonic physics on the mass distribution ofMW-like systems in Illustris are substantial: by matching haloesbetween the DMO and hydrodynamic simulations, we find that thelatter have more mass on small scales and less mass on large scales.The asymptotic difference in the total mass density at large radii isapproximately 20%. • Since different feedback models result in very different effectson the mass distribution of dark matter even at large distances fromhalo centres (e.g., fig. 7 of Vogelsberger et al. 2014a compared tofig. 1 of Schaller et al. 2015), it is imperative to test how inferenceson the mass of the MW depend on galaxy formation modelling. • The mass distribution in the inner ∼ kpc is not convergedin the Illustris suite [see Schaller et al. 2016 for similar results inthe Eagle simulations]; this is a much larger distance than the for-mal convergence radius for the dark matter simulations. Resultsregarding the density distribution for r (cid:46) kpc must thereforebe interpreted with caution, and our best-fitting NFW profiles for the hydrodynamic simulations, which were obtained over the radialrange of 40-300 kpc, should not be extrapolated to smaller radii. • The relationship between M ( < 50 kpc) and M , c inIllustris-1 is well-described by a log-quadratic relationship (equa-tion. 2). This relationship enables the translation of any existing orfuture constraint on M ( < 50 kpc) into a measurement M , c . • The constraints on M ( < 50 kpc) derived by D12 ( . ± . × M (cid:12) ) and G14 ( . ± . × M (cid:12) ) predict very differentvalues for the virial mass of the Galaxy’s halo when using Illustris:for D12, we find M , c = 1 . +0 . − . × M (cid:12) (68% confi-dence), while for G14, we find M , c = 0 . +0 . − . × M (cid:12) (68% confidence). The values for M vir and M , m are 17% and32% larger, respectively. • Illustris haloes that have galaxies with stellar masses consis-tent with measurements of the MW’s M (cid:63) have significantly moremass within 50 kpc than the result of G14; the measurements ofD12 and Williams & Evans (2015) are in much better agreementwith Illustris haloes that match the observed value of M (cid:63) . In partic-ular, almost no haloes in Illustris jointly satisfy the G14 constraintand the LN15 measurement of M (cid:63) for the MW. • From our analysis of the Illustris simulation, even an infinitelyprecise measurement of M ( < 50 kpc) would result in an uncer-tainty of > 20% in M , c . The same uncertainty can be achievedfor 10% errors on M ( < 80 kpc) or 12% errors on M ( < 100 kpc) .A measurement of M ( < 100 kpc) that is accurate to 5% will trans-late into 15% uncertainties on M , c .As ever larger and ever more realistic hydrodynamic simula-tions become available, so too will better statistical constraints onthe mass profile of our Galaxy. ACKNOWLEDGEMENTS We thank Joss Bland-Hawthorn, Nitya Kallivayalil, Julio Navarro,and Annalisa Pillepich for helpful conversations. The analysis ofthe Illustris data sets for this paper was done using the Odysseycluster, which is supported by the FAS Division of Science, Re-search Computing Group at Harvard University. MB-K acknowl-edges support provided by NASA through a Hubble Space Tele-scope theory grant (programme AR-12836) from the Space Tele-scope Science Institute (STScI), which is operated by the Associa-tion of Universities for Research in Astronomy (AURA), Inc., un-der NASA contract NAS5-26555. PT acknowledges support fromNASA ATP Grant NNX14AH35G. LH acknowledges support fromNASA grant NNX12AC67G and NSF grant AST-1312095. REFERENCES Barber C., Starkenburg E., Navarro J. F., McConnachie A. W., Fattahi A.,2014, MNRAS, 437, 959Battaglia G., et al., 2005, MNRAS, 364, 433Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57Besla G., Kallivayalil N., Hernquist L., Robertson B., Cox T. J., van derMarel R. P., Alcock C., 2007, ApJ, 668, 949Bovy J., Rix H.-W., 2013, ApJ, 779, 115Bovy J., et al., 2012, ApJ, 759, 131Boylan-Kolchin M., Besla G., Hernquist L., 2011a, MNRAS, 414, 1560Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011b, MNRAS, 415,L40Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2012, MNRAS, 422,1203 MNRAS , 3483–3493 (2016) ass of the Milky Way from Illustris Boylan-Kolchin M., Bullock J. S., Sohn S. T., Besla G., van der Marel R. P.,2013, ApJ, 768, 140Brown W. R., Geller M. J., Kenyon S. J., Diaferio A., 2010, AJ, 139, 59Bryan G. L., Norman M. L., 1998, ApJ, 495, 80Bullock J. S., Kolatt T. S., Sigad Y., Somerville R. S., Kravtsov A. V.,Klypin A. A., Primack J. R., Dekel A., 2001, MNRAS, 321, 559Busha M. T., Wechsler R. H., Behroozi P. S., Gerke B. F., Klypin A. A.,Primack J. R., 2011, ApJ, 743, 117Cohen J. G., Sesar B., Banholzer S., PTF Consortium t., 2015,arXiv:1509.05997,Courteau S., Dutton A. A., 2015, ApJ, 801, L20Cunningham E. C., Deason A., Guhathakurta P., Rockosi C., Kirby E., vander marel r. p., Sohn S. T., 2015, IAU General Assembly, 22, 55864Deason A. J., Belokurov V., Evans N. W., An J., 2012, MNRAS, 424, L44Dutton A. A., Macci`o A. V., 2014, MNRAS, 441, 3359Eadie G. M., Harris W. E., Widrow L. M., 2015, ApJ, 806, 54Fang T., Bullock J., Boylan-Kolchin M., 2013, ApJ, 762, 20Fattahi A., et al., 2016, MNRAS, 457, 844Fukugita M., Peebles P. J. E., 2004, ApJ, 616, 643Genel S., et al., 2014, MNRAS, 445, 175Gibbons S. L. J., Belokurov V., Evans N. W., 2014, MNRAS, 445, 3788Gnedin O. Y., Brown W. R., Geller M. J., Kenyon S. J., 2010, ApJ, 720,L108G´omez F. A., Besla G., Carpintero D. D., Villalobos ´A., O’Shea B. W., BellE. F., 2015, ApJ, 802, 128Gonz´alez R. E., Kravtsov A. V., Gnedin N. Y., 2013, ApJ, 770, 96Gonz´alez R. E., Kravtsov A. V., Gnedin N. Y., 2014, ApJ, 793, 91Guo Q., White S., Li C., Boylan-Kolchin M., 2010, MNRAS, 404, 1111Gupta A., Mathur S., Krongold Y., Nicastro F., Galeazzi M., 2012, ApJ,756, L8Hinshaw G., et al., 2013, ApJS, 208, 19Jiang F., van den Bosch F. C., 2015, MNRAS, 453, 3575Kafle P. R., Sharma S., Lewis G. F., Bland-Hawthorn J., 2014, ApJ, 794, 59Kahn F. D., Woltjer L., 1959, ApJ, 130, 705Kallivayalil N., van der Marel R. P., Alcock C., Axelrod T., Cook K. H.,Drake A. J., Geha M., 2006, ApJ, 638, 772Kallivayalil N., van der Marel R. P., Besla G., Anderson J., Alcock C., 2013,ApJ, 764, 161Li Y.-S., White S. D. M., 2008, MNRAS, 384, 1459Licquia T. C., Newman J. A., 2015, ApJ, 806, 96Mandelbaum R., Wang W., Zu Y., White S., Henriques B., More S., 2016,MNRAS, 457, 3200McMillan P. J., 2011, MNRAS, 414, 2446Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Nelson D., et al., 2015, arXiv:1504.00362,Piatek S., Pryor C., Bristow P., Olszewski E. W., Harris H. C., Mateo M.,Minniti D., Tinney C. G., 2007, AJ, 133, 818Piffl T., et al., 2014, A&A, 562, A91Pryor C., Piatek S., Olszewski E. W., 2015, AJ, 149, 42Rashkov V., Pillepich A., Deason A. J., Madau P., Rockosi C. M., GuedesJ., Mayer L., 2013, ApJ, 773, L32Schaller M., et al., 2015, MNRAS, 451, 1247Schaller M., et al., 2016, MNRAS, 455, 4442Schaye J., et al., 2015, MNRAS, 446, 521Sch¨onrich R., 2012, MNRAS, 427, 274Sohn S. T., Besla G., van der Marel R. P., Boylan-Kolchin M., MajewskiS. R., Bullock J. S., 2013, ApJ, 768, 139Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS,328, 726Torrey P., Vogelsberger M., Genel S., Sijacki D., Springel V., Hernquist L.,2014, MNRAS, 438, 1985Vera-Ciro C., Helmi A., 2013, ApJ, 773, L4Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L.,2013, MNRAS, 436, 3031Vogelsberger M., et al., 2014a, MNRAS, 444, 1518Vogelsberger M., et al., 2014b, Nature, 509, 177 Wang J., Frenk C. S., Navarro J. F., Gao L., Sawala T., 2012, MNRAS, 424,2715Wang W., Han J., Cooper A. P., Cole S., Frenk C., Lowing B., 2015, MN-RAS, 453, 377Watkins L. L., Evans N. W., An J. H., 2010, MNRAS, 406, 264Wilkinson M. I., Evans N. W., 1999, MNRAS, 310, 645Williams A. A., Evans N. W., 2015, MNRAS, 454, 698Xue X. X., et al., 2008, ApJ, 684, 1143de Bruijne J. H. J., Rygl K. L. J., Antoja T., 2014, in EAS PublicationsSeries. pp 23–29 ( arXiv:1502.00791 ), doi:10.1051/eas/1567004van der Marel R. P., et al., 2014, in Seigar M. S., Treuthardt P., eds, As-tronomical Society of the Pacific Conference Series Vol. 480, Structureand Dynamics of Disk Galaxies. p. 43 ( arXiv:1309.2014 )MNRAS461