Constraining the Milky Way Mass Profile with Phase-Space Distribution of
  Satellite Galaxies
  We estimate the Milky Way (MW) halo properties using satellite kinematic data
including the latest measurements from Gaia DR2. With a simulation-based 6D
phase-space distribution function (DF) of satellite kinematics, we can infer
halo properties efficiently and without bias, and handle the selection function
and measurement errors rigorously in the Bayesian framework. Applying our DF
from the EAGLE simulation to 28 satellites, we obtain an MW halo mass of
$M=1.23_{-0.18}^{+0.21}\times 10^{12} M_\odot$ and a concentration of $c=9.4_{
-2.1}^{ +2.8}$ with the prior based on the $M$-$c$ relation. The inferred mass
profile is consistent with previous measurements but with better precision and
reliability due to the improved methodology and data. Potential improvement is
illustrated by combining satellite data and stellar rotation curves. Using our
EAGLE DF and best-fit MW potential, we provide much more precise estimates of
kinematics for those satellites with uncertain measurements. Compared to the
EAGLE DF, which matches the observed satellite kinematics very well, the DF
from the semi-analytical model based on the dark-matter-only simulation
Millennium II (SAM-MII) over-represents satellites with small radii and
velocities. We attribute this difference to less disruption of satellites with
small pericenter distances in the SAM-MII simulation. By varying the disruption
rate of such satellites in this simulation, we estimate a $\sim 5\%$ scatter in
the inferred MW halo mass among hydrodynamics-based simulations.
