  We for the first time provide self-consistent axisymmetric phase-space
distribution models for the Milky Way's dark matter (DM) halo which are
carefully matched against the latest kinematic measurements through Bayesian
analysis. By using broad priors on the individual galactic components, we
derive conservative estimates for the astrophysical factors entering the
interpretation of direct and indirect DM searches. While the resulting DM
density profiles are in good agreement with previous studies, implying
$\rho_\odot \approx 10^{-2} \, M_\odot / \mathrm{pc}^3$, the presence of
baryonic disc leads to significant differences in the local DM velocity
distribution in comparison with the standard halo model. For direct detection,
this implies roughly 30% stronger cross-section limits at DM masses near
detectors maximum sensitivity and up to an order of magnitude weaker limits at
the lower end of the mass range. Furthermore, by performing Monte-Carlo
simulations for the upcoming DARWIN and DarkSide-20k experiments, we
demonstrate that upon successful detection of heavy DM with coupling just below
the current limits, the carefully constructed axisymmetric models can eliminate
bias and reduce uncertainties by more then 50% in the reconstructed DM coupling
and mass, but also help in a more reliable determination of the scattering
operator. Furthermore, the velocity anisotropies induced by the baryonic disc
can lead to significantly larger annual modulation amplitude and sizable
differences in the directional distribution of the expected DM-induced events.
For indirect searches, we provide the differential $J$-factors and compute
several moments of the relative velocity distribution that are needed for
predicting the rate of velocity-dependent annihilations. However, we find that
accurate predictions are still hindered by large uncertainties regarding the DM
distribution near the galactic center.
