We use a distribution function analysis to estimate the mass of the Milky Way out to 100 kpc using a large sample of halo stars. These stars are compiled from the literature, and the vast majority (~98%) have 6D phase-space information. We pay particular attention to systematic effects, such as the dynamical influence of the Large Magellanic Cloud (LMC), and the effect of unrelaxed substructure. The LMC biases the (pre-LMC infall) halo mass estimates towards higher values, while realistic stellar halos from cosmological simulations tend to underestimate the true halo mass. After applying our method to the Milky Way data we find a mass within 100 kpc of M(< 100 kpc) = 6.07 +/- 0.29 (stat.) +/- 1.21 (sys.) x 10^11 M_Sun. For this estimate, we have approximately corrected for the reflex motion induced by the LMC using the Erkal et al. model, which assumes a rigid potential for the LMC and MW. Furthermore, stars that likely belong to the Sagittarius stream are removed, and we include a 5% systematic bias, and a 20% systematic uncertainty based on our tests with cosmological simulations. Assuming the mass-concentration relation for Navarro-Frenk-White haloes, our mass estimate favours a total (pre-LMC infall) Milky Way mass of M_200c = 1.01 +/- 0.24 x 10^12 M_Sun, or (post-LMC infall) mass of M_200c = 1.16 +/- 0.24 x 10^12 M_Sun when a 1.5 x 10^11 M_Sun mass of a rigid LMC is included.
Download