§ INTRODUCTION

Dark matter accounts for more than eighty percent of the matter in the universe <cit.>, but has managed to evade efforts to determine its precise makeup.
Among the most popular ideas is that dark matter is made up of as-of-yet undiscovered elementary particles of nature, and there are significant efforts underway to directly detect new particles with the characteristics necessary for them to be dark matter.

While so far direct detection attempts have been unsuccessful, the lack of signal provides important constraints on the allowed microphysical properties of theoretical dark matter particles.  In particular, direct detection experiments place joint limits on the allowed mass and nucleon-interaction cross section for particle dark matter (e.g. <cit.>). Two astrophysical assumptions are critical for deriving these constraints (see ref. <cit.> for a nice discussion): 1) the local mass density of dark matter in the Earth's vicinity, $\rho$, and 2) the local speed distribution of dark matter particles, $f(v)$. Moreover, the interpretation of any future detection signal from such an experiment will rely heavily on both $\rho$ and $f(v)$.  These facts motivate considerable effort to determine local dark matter properties as accurately as possible. Furthermore, in indirect detection, the dark matter density is fundamental in determining whether a WIMP dark matter model for the Galactic Center extended gamma-ray excess is compatible with the lack of annihilation signal from dwarf galaxies (e.g. <cit.>).

Methods for determining the local dark matter density fall into two broad categories: Local and Global.  Local measures rely on tracers in the Solar vicinity, while Global estimates rely on fitting a density model to measurements taken throughout the Milky Way. See ref. <cit.> for a review of the various methods. We summarize eleven estimates of $\rho$ in Table <ref> and see that even among measurements that adopt similar approaches, estimates for the local density vary more than can be explained by quoted statistical error bars. This suggests that there are systematic errors at play. There is also a trend for Local methods to give slightly higher dark matter densities than Global methods, with averages of $\rho = 0.47 \pm 0.11\ \denunit$ for Local and $\rho = 0.39 \pm 0.10\ \denunit$ for Global[The error bars on these averages are based on the standard deviation of the $N$ measurements in each group, weighted by $\sqrt{1 + 1/N}$.]. However, this is not an exhaustive literature search; more work could be done to determine whether there is a true bias with a physical basis. Ref. <cit.> found that assuming the Milky Way's disk is axisymmetric and in a steady state, as most Local measures do, can impart ${\sim}30\%$ uncertainty. Ref. <cit.>, which takes into account the Galaxy's time-varying structure and phase-space spiral, produces the lowest density estimate of the Local measurements we highlight. It will be interesting to see whether further Local measures that allow for time-varying structure and axi-asymmetry produce similar results. We further note that ref. <cit.> adopted an approach qualitatively similar to the method we present below, in that they rely on only the local equation of centrifugal equilibrium, not the full rotation curve.  They find a preferred density of $\rho \simeq 0.43 \, \denunit$, which is similar to average of all the studies summarized in Table <ref>:   $\rho = 0.41 \pm 0.10 \, \denunit$.

Eleven of the most recent studies measuring the local dark matter density $\rho$, grouped by Global and Local methods. The margins of error for the averages are calculated as $s\sqrt{1+1/N}$ where $s$ is the standard deviation of the measurements, and $N$ is the number of measurements. In deSales+19 <cit.>, B1 and B2 refer to two different baryonic mass models.

4cGlobal methods

Model basis                                                 Note        $\rho\,/\,[\denunit]$

Eilers+19 <cit.>               Fitting the circular speed curve                                     $0.3\pm0.03$

deSales+19 <cit.>       Fitting the circular speed curve                         B1          $0.3\pm0.03$

B2          $0.38\pm0.04$

Nitschai+20 <cit.>                  Anisotropic jeans modeling of disk stars                                $0.437\pm0.076$

Petac20 <cit.>                   Galactic mass model with circular speed as one input     NFW         $0.357^{+0.02}_{-0.021}$

Burkert     $0.381^{+0.02}_{-0.022}$

Benito+21 <cit.>              Fitting the circular speed curve                                     $0.6\pm0.1$

Hattori+21 <cit.>             Galactic mass model from halo stars                                     $0.342\pm0.007$

Ou+24 <cit.>                           Fitting the circular speed curve                                     $0.447\pm0.004$

Average                                                                                $0.39\pm0.10$

4cLocal methods

Model basis                                                 Note        $\rho\,/\,[\denunit]$

Sivertsson+18 <cit.>          Vertical kinematics                                                     $0.46^{+0.07}_{-0.09}$

Buch+19 <cit.>                          Very local kinematics                                       A stars     $0.61\pm0.38$

G stars     $0.42^{+0.38}_{-0.34}$

Guo+20 <cit.>                        Vertical kinematics                                                     $0.50^{+0.09}_{-0.08}$

Salomon+20 <cit.>               Vertical kinematics                                                     $0.37\pm0.09$

Widmark+21 <cit.>                 Time-varying structure of the local phase-space spiral                  $0.32\pm0.15$

Average                                                                                $0.47\pm0.11$

Average of all above                                                                        $0.41\pm0.10\,\denunit$

This paper                                                                       $\boldsymbol{\data{rho_GeV}\pm\data{drho_GeV}\,\denunit}$

For the speed distribution, the standard has been to assume a Maxwellian (Maxwell-Boltzmann) distribution that peaks at the Milky Way's circular speed and truncates at the escape speed (e.g. <cit.>).
This assumption is often called the “standard halo model" or SHM, usually with a circular speed of $220\,\mathrm{km\,s^{-1}}$. There has been extensive research into whether the SHM appropriately represents the Milky Way or even whether a Maxwellian is the correct functional form to adopt. Investigators typically carry out this evaluation on specific simulated galaxies, chosen from a larger pool of candidates based on how closely they match certain properties of the Milky Way <cit.>.

While a few papers using high-resolution galaxy formation simulations have found speed distributions consistent with Maxwellians <cit.>,
a Maxwellian usually fails to reproduce the shape of dark matter speed distributions as seen in modern simulations, usually over-predicting the number of high-velocity particles <cit.>.

These studies tend to achieve better fits with more complex models, such as generalized Maxwellian distributions, Tsallis distributions, Eddington formulas, or various other modifications to the Maxwellian distribution such as those developed by ref. <cit.> or <cit.>. The prevailing methodology has been to then analyze data from dark matter direct detection experiments assuming that the speed distribution for a specific simulated Milky Way analogue represents the Milky Way itself.

In this paper we analyze twelve zoom-in simulations of Milky Way size disk galaxies from the FIRE-2 project and attempt to develop formulaic methods that allow us to infer local dark matter properties of the Milky Way itself. We obviate the need to choose one particular simulated galaxy as a Milky Way analogue. Instead, we correlate the density, velocity dispersion, and speed distribution using a sliding scale of one single, simple observable: the local circular speed.

In the next section we provide a short overview of dark matter detection and introduce some nomenclature.  In 3 we discuss our simulations, and in 4 we describe our analysis.  All of our results are presented in 5, and 6 is reserved for discussion and conclusions.

§ BACKGROUND AND NOMENCLATURE

The differential event rate for interactions between dark matter particles and a detector is

\begin{equation}
\label{eq:drde}
\dv{\mathcal{R}}{E_\mathrm{R}}=\frac{\rho}{m_\chi m_\mathrm{N}}\int_{v_\mathrm{min}}^\infty vf(v)\dv{\sigma_\mathrm{WN}(v,E_\mathrm{R})}{E_\mathrm{R}}\mathrm{d}v ,
\end{equation}

where $\rho$ is the local dark matter density, $m_\chi$ is the mass of the dark matter particle, $m_\mathrm{N}$ is the mass of the nucleus used in the detector, $\sigma_\mathrm{WN}$ is the interaction cross section between the nucleus and the dark matter particle, and $v_\mathrm{min}$ is the minimum dark matter particle speed required for the detector nucleus to recoil with energy $E_\mathrm{R}$. The value of $v_\mathrm{min}$  is given by

\begin{equation}
\left(\frac{v_\mathrm{min}}{c}\right)^2={\frac{m_\mathrm{N}E_\mathrm{R}}{2\mu^2}}\left(1+\frac{\mu}{m_\mathrm{N}E_\mathrm{R}}\delta\right)^2 .
\end{equation}

Here, $\mu$ is the reduced mass between the dark matter particle and the detector nucleon, and $\delta$ encodes the mass splitting that may occur during an interaction. For standard and low-mass dark matter, $\delta=0$; for inelastic DM, $\delta=\mathcal{O}(100\,\mathrm{keV})$ <cit.>.
Note that for the canonical leading order spin-independent and spin-dependent DM-nucleus interactions, the differential cross section is inversely proportional to the square of the WIMP's speed $\left(\mathrm{d}\sigma_\mathrm{WN}/\mathrm{d}E_\mathrm{R}\propto v^{-2}\right)$ <cit.>. Therefore, the differential event rate's dependence on velocity is entirely captured by the halo integral

\begin{equation}
\label{eqn:halo_integral}
g(v_\mathrm{min}) = \int_{v_\mathrm{min}}^\infty\frac{f(v)}{v}\mathrm{d}^3v.
\end{equation}

A traditional assumption for the speed distribution of dark matter particles is a Maxwellian:

\begin{equation}
\label{eqn:simple_max}
\end{equation}

where $v = |\vec{v}|$,  $v_0$ is the most-probable speed or “peak speed", and
$ N(v_0) = \pi^{3/2}v_0^3$ is a normalization factor that gives $\int_0^\infty f(v)4\pi v^2dv=1$. It is common to set the peak speed equal to the circular speed of the Milky Way at the Solar radius: $v_0=v_\mathrm{c}$. Alternatively, $v_0$ can be expressed in terms of the 3D velocity dispersion:

\begin{equation}
v_0 = \sqrt{2/3}\,\sigma_\mathrm{3D}.
\label{eqn:23sigma}
\end{equation}

A slight modification of the Maxwellian assumption is to introduce a sharp truncation in Equation <ref> above some escape speed, $v_\mathrm{esc}$.  In this case, the normalization factor must be modified:

\begin{equation}
N(v_0) = \pi v_0^2
\left[
\sqrt{\pi}v_0\erf
\left(\frac{v_\mathrm{esc}}{v_0}\right)
\end{equation}

An alternative to Equation <ref> motivated by galaxy formation simulations is a parameterization presented by Mao et al. (2013) <cit.>:

\begin{equation}
f(v) =
\begin{cases}
N\exp(-\frac{v}{v_0})(v^2_\mathrm{esc}-v^2)^P, & 0 \leq v\leq v_\mathrm{esc}\\
0, & \text{otherwise}.
\end{cases}
\label{eqn:mao}
\end{equation}

In Section <ref>, we introduce a new model for $f(v)$ and tie its parameters to the circular speed. We then use the Milky Way's observed circular speed to determine its most likely parameters. In the Appendix Section <ref>, we provide a comparison between our model and the functional form in Equation <ref>.

§ SIMULATIONS

We analyze cosmological zoom-in simulations run with the gizmo code <cit.> from the Feedback in Realistic Environments (FIRE) project[<https://fire.northwestern.edu>]. Specifically, we use the FIRE-2 model <cit.> to implement stellar feedback and star formation. An effective feedback implementation is important for regulating star formation, producing galaxies that are the correct mass for their dark matter halo size, and also in enabling the formation of galaxies with disk morphologies <cit.>. Additionally, feedback can affect galaxies' density profiles via “feedback-induced core formation" (e.g. <cit.> and references therein).  Galaxy formation is also known to increase the central dark matter velocity dispersion compared to dark-matter-only simulations <cit.>, thus achieving realistic galaxy masses is crucial to producing dark matter properties that are accurate. The stellar feedback in FIRE includes stellar winds, radiation pressure, photoelectric heating, and photoheating from ionizing radiation. Each star particle assumes a Kroupa IMF <cit.>, has an age determined from its formation time, and inherits its metallicity from its parent gas particle. Within the gas particles, the simulations track 11 elemental abundances (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) with sub-grid diffusion via turbulence <cit.>.  The conditions for star formation in gas are local self gravitation, sufficient density (${>}1000\rm\, cm^{-3}$), Jeans instability, and molecularity (following ref. <cit.>).

We focus on twelve galaxies with prominent disk components and stellar masses similar to that of the Milky Way. The initialization of these galaxies follows ref. <cit.> using the music package <cit.>. Six of the twelve come from the Latte suite <cit.> and six from the ELVIS on FIRE project <cit.>. The former is a collection of zooms on isolated galaxies with an m12 naming convention. The latter consists of pairs with separations and relative velocities similar to the MW-m31 pair. Each ELVIS pair comprises a single zoom, and they follow the naming convention of famous duos. We exclude galaxies such as m12z and m12w whose fraction of stellar mass in the thin disk is only 5% and 23% respectively. m12r exhibits a similarly low thin-disk fraction of 15%. However, $\sqrt{GM/R_0}$ approximates the observable circular speed of cold gas in m12r as well as those galaxies with higher thin-disk fractions. This indicates some level of diskiness. Therefore, we include m12r in our analysis. The thin-disk fractions of the eleven other disks range from 29% for Thelma to 43% for Romeo.
Dark matter particle mass is $3.5\times10^4\ \rm M_\odot$ in Latte and $\simeq2\times 10^4\ \rm M_\odot$ in ELVIS. Gas and star particles in Latte have initial masses of $7100\ \rm M_\odot$ while ELVIS simulations have approximately twice the resolution at 3500 to 4000 $\rm M_\odot$. Gas softening lengths are fully adaptive down to $\simeq0.5$–$1\,\rm pc$. The star-particle softening length is ${\simeq}4\,\rm pc$ physical. The dark matter force softening is ${\simeq}40\,\rm pc$ physical.

§ METHODS

§.§ Defining Sun-like Regions in the Simulations

We are interested in exploring dark matter characteristics around mock Sun-like regions in the disks of our simulated galaxies.  In each simulation, the plane of the disk is defined as that perpendicular to to the aggregate angular momentum of all gas, stars, and dark matter within 10 kpc of the galaxy's center. We find that the orientation of disk planes remains virtually identical if we include only stars when determining the angular momentum direction.

For analysis we use particles in ring-like regions centered on mock solar radii ($R_0 = 8.3$ kpc) in each simulated disk. One could add a level of complexity to this study by scaling $R_0$ for each galaxy. However, the results below show that a fixed radius yields tight constraints on our conclusions, which we find satisfactory. We define each galaxy's solar ring by first taking a spherical shell 1.5 kpc thick with a midpoint radius 8.3 kpc from the center of the given galaxy. From there, we exclude the parts of the shell 0.5 kpc above and below the plane of the disk. These choices for thickness and height allow the ring to be as small as possible while still allowing appreciable particle counts.

Azimuthal variation in dark matter density (left) and velocity dispersion (right) across 30 different Sun-like regions in each of our simulated galaxies.  In both plots, measurements were made in 30 evenly-divided regions in solar rings for each galaxy and then normalized by the average density ($\overline{\rho}$) or velocity dispersion ($\overline{\sigma}_{\rm3D}$) over the full ring.    Across all galaxies, the standard deviation of $\rho(\phi)/\overline{\rho}$ is stdev_linear_dendiff, and the standard deviation of $\sigma_{\rm3D}(\phi)/\overline{\sigma}_{\rm3D}$ is  stdev_linear_dispdiff. We adopt the aggregate solar ring values as our simulated measurements, so $\rho=\overline{\rho}$, and $\sigma_{\rm3D}=\overline{\sigma}_{\rm3D}$ throughout the text and in Figures <ref> and <ref>.  We use these quoted standard deviations to inform our uncertainties in predicting local dark matter properties near the Sun.

§.§ Circular Speed

All the models in this study make predictions based on a sliding scale of one simple observable: the galaxy's observed circular rotation speed at the solar location ($R_0=8.3$ kpc). In an idealized, spherically-symmetric system, $v_\mathrm{c}=\sqrt{GM/R_0}$. However, in reality, the Milky Way is a flattened disk, and material in such a disk will not rotate at this idealized speed <cit.>.
We determine the circular speeds of our galaxies as the average azimuthal velocity component $\vec{v}\cdot\hat{\phi}$ of gas colder than $10^4$ K in the solar ring, where $\hat{\phi}$ points in the direction of the galaxy's rotation. This temperature ensures that the observed gas is in circular motion, because the implied ${\sim}16\rm\,km s^{-1}$ 3D velocity dispersion is an order of magnitude smaller than $v_{\rm c}$. Note that we find that these “observed" circular speeds end up being very close to what we would have obtained with the idealized spherical approximation. For example, Figure <ref> in the Appendix plots our mock “observed" circular speeds versus the idealized spherical approximation and shows that they obey close to a one-to-one relationship.
Additionally, ref. <cit.> found that, for Milky Way size galaxies in FIRE, observational methods exhibit excellent agreement with the true rotation curve. Even the three largest m11's have observable rotation curves that match the true curve to within 10%.
This suggests that our results will be largely insensitive to how we define $v_\mathrm{c}$ in the simulations.

Ultimately, our goal is to use our inferred relationships between  simulated dark matter properties and $v_{\rm c}$ on the Milky Way itself. In doing so, we adopt $v_\mathrm{c,MW}=229\pm7\ \mathrm{km\,s^{-1}}$ for the Milky Way's rotation speed near the Sun <cit.>.  This error bar includes systematic uncertainties provided by ref. <cit.>. Even so, at 3%, it contributes little to the overall uncertainty in our predictions, which is dominated by galaxy-to-galaxy variations.

§.§ Density and Dispersion Variation within Disks

As discussed above, we use ring-like regions centered on mock solar locations to determine the local density $\rho$, velocity dispersion $\sigma_\mathrm{3D}$, and speed distributions of dark matter particles in each simulation. It is interesting to ask if there is any significant variation in density or velocity dispersion as we move around the plane of the disks at fixed radii $R_0$.

In order to explore this, we divided each ring into 30 $\phi$ slices and calculated $\rho(\phi)$ and $\sigma_\mathrm{3D}(\phi)$ in each piece.
Figure <ref> shows histograms of those measurements at fixed $\phi$ divided by the corresponding galaxy's average for the whole ring.  Density variation is shown on the left, and velocity dispersion variation on the right.  The standard deviation of a $\phi$ slice from the aggregate solar ring is stdev_linear_dendiff for $\rho$ and stdev_linear_dispdiff for $\sigma_\mathrm{3D}$. We conservatively add these standard deviations in quadrature to the other uncertainties in predicting the Milky Way's local dark matter density and velocity dispersion in what follows. Note that this level of variation is quite small and adds an almost negligible amount to the error bar around the prediction.
Moreover, the number of dark matter particles in our twelve FIRE-2 solar rings ranges from 15,881 for m12r to 51,785 for Romulus.
With 30 slices, this implies
shot noise at the ${\sim}2$–$4\%$ level—not far from the standard deviation we measure directly among density slices.  This suggests that our reported variation of stdev_linear_dendiff in density is a conservative upper limit.

We compared shot noise to the standard deviation of $\rho/\overline{\rho}$ to determine the optimal number of slices. With more than ${\sim}30$ slices, shot noise exceeds the standard deviation between slices. In other words, counting statistics could be driving apparent azimuthal variations. We decreased the number of slices in multiples of 5 until standard deviation exceeded shot noise for all galaxies, which occurs at 30.
Also note that our results do not appreciably change if we use an extreme number of slices. Using 200 for example, the resulting 8.9% and 3.0% standard deviations in density and dispersion remain nearly negligible contributors to overall uncertainty in Milky Way predictions. The main contributor to said uncertainty is galaxy-to-galaxy variation.

§.§ Defining Escape Speed

As discussed above, it is common to truncate an assumed dark matter particle speed distribution above some “escape speed" $v_{\rm esc}$, but knowing the precise value of the escape speed to use is non-trivial.  Even if the the full potential is known, the formal escape speed depends on the distance[Qualitatively, in order for a dark matter particle to “escape" the Solar region, it does not need to escape to infinity, but just how far is not entirely clear.] a particle needs to travel in order to “escape" the region of interest (e.g. <cit.>).

In what follows, we take a practical approach and define the simulated escape speed to be the value of $v$ where the local dark matter speed distribution falls to zero. To determine the sensitivity of this measure to resolution and Poisson noise, we tested the variability of the results by sampling only 10% of each galaxy's particles and repeating the analysis ten times. The standard deviation between each galaxy's ten trials ranges from $\data{vesc_min_std}{\sim}\data{vesc_max_std}$. Therefore, the study's measure of escape speed is acceptably robust.
We also compare our “measured" escape speed to the formal escape speed to infinity: $v_\mathrm{esc}(\Phi)=\sqrt{2|\Phi|}$. For a given galaxy, we calculate this at eight points in the solar ring, equidistant in azimuthal angle $\phi$ and take the average of those eight values.
As shown below, this categorically overshoots the measured $v_{\rm esc}$ in our simulations, which is consistent with previous findings (e.g. <cit.>).

§ RESULTS

Relationship between the local dark matter density and circular speed for our simulated galaxies.  The name of each simulated galaxy (as defined in the original FIRE-2 papers) is provided for each respective point.  The color bar maps to the total stellar mass of each galaxy for reference. The blue line shows the best fit from Equation <ref>, which yields $\rho \sim v_{\rm c}^{\data{den_slope}}$ and a log-space coefficient of determination $r^2$ of $0.65$. The black error bars on each data point come from our measured azimuthal deviations in local density, discussed in Section <ref>. The grey band represents the 1$\sigma$ uncertainty around the prediction line. The red point and error bar represents our Milky Way prediction: $\rho(v_\mathrm{c,MW})=\data{rho_1e7msun}^{+\data{drho_1e7msun_plus}}_{-\data{drho_1e7msun_minus}}\,10^7\mathrm{M_\odot\,kpc^{-3}}=\data{rho_GeV}\pm\data{drho_GeV}\,\denunit$.

§.§ Local Dark Matter Density

Figure <ref> shows the relationship between the measured “local” dark matter densities at mock solar locations in our simulations versus the local circular speed. We find that the relationship between density and $v_\mathrm{c}$ is well described by a power law:

\begin{equation}
\rho=\rho_0\left(\frac{v_\mathrm{c}}{100\,\mathrm{km\,s^{-1}}}\right)^{\alpha},
\label{eqn:rho}
\end{equation}

where $\rho_0=\data{rho0_1e7msun}^{+\data{drho0_1e7msun_plus}}_{-\data{drho0_1e7msun_minus}}10^7\mathrm{M_\odot\,kpc^{-3}}=\data{rho0_GeV}^{+\data{drho0_GeV_plus}}_{-\data{drho0_GeV_minus}}\denunit$ and $\alpha=\data{den_slope}\pm\data{dden_slope}$.
Section <ref> in the Appendix discusses our treatment of uncertainties in the model.
The solid blue line in Figure <ref> visualizes this density model, which exhibits a log-space coefficient of determination $r^2$ of 0.65. The shaded band shows the one-sigma region, and the color bar maps to stellar mass of the galaxy.

If we apply our model fit to the Milky Way, assuming a circular speed of $v_\mathrm{c,MW}=229\pm7\ \mathrm{km\,s^{-1}}$, from ref. <cit.>, this yields

\begin{align}
\rho(v_\mathrm{c,MW})&=\data{rho_1e7msun}^{+\data{drho_1e7msun_plus}}_{-\data{drho_1e7msun_minus}}10^7\mathrm{M_\odot\,kpc^{-3}}\label{eqn:rho_mw_astro}\\
\end{align}

The quoted error includes the uncertainty from the fit as well as allowance for variance within the disk. The latter is equal to the azimuthal variation found in our simulations in Section <ref>.
Interestingly, our determination of the local density is quite close to the aggregate values in Table <ref>. Not only does this cross-validate this study's density result with others', it also supports using these select FIRE-2 galaxies for making other Milky Way predictions. Note that one could easily adopt a different value for $v_\mathrm{c,MW}$ through Equation <ref>.

This is the same as Figure <ref> but now showing measured velocity dispersion versus circular speed for each simulated galaxy. The best-fit model (Equation <ref>, blue line) has $\sigma \sim v_{\rm c}^{\data{disp_slope}}$ and  a log-space coefficient of determination $r^2$ of 0.90. The relationship between local dark matter velocity dispersion and local circular speed is even tighter than that between the latter and density shown in Figure <ref>. The red point and error bar shows our Milky Way prediction: $\sigma_{\rm 3D}(v_\mathrm{c,MW})=\data{disp}^{+\data{ddisp_plus}}_{-\data{ddisp_minus}}\,\mathrm{km\,s^{-1}}$.

§.§ Local Dark Matter Velocity Dispersion

We find that the dark matter velocity dispersion is even more tightly correlated with $v_\mathrm{c}$ than the density, with

\begin{equation}
\sigma_\mathrm{3D}=\sigma_0\left(\frac{v_\mathrm{c}}{100\,\mathrm{km\,s^{-1}}}\right)^{\gamma},
\label{eqn:sigma_vs_vc}
\end{equation}

where $\sigma_0=\data{disp_amp}^{+\data{ddisp_amp_plus}}_{-\data{ddisp_amp_minus}}\,\velunit$, and $\gamma=\data{disp_slope}\pm\data{ddisp_slope}$, which exhibits a 0.90 log-space $r^2$ coefficient of determination. For the Milky Way, this yields

\begin{equation}
\sigma_\mathrm{3D}(v_\mathrm{c,MW})=\data{disp}^{+\data{ddisp_plus}}_{-\data{ddisp_minus}}\,\mathrm{km\,s^{-1}}.
\label{eqn:sigma_val}
\end{equation}

Figure <ref> visualizes the model and its Milky Way prediction versus the data. Similarly to our density model, one could adopt a different MW circular speed by simply plugging their preferred value into Equation <ref>. Again, when making predictions around $v_\mathrm{c}\approx229\,\mathrm{km\,s^{-1}}$, one can assume the corresponding uncertainties are the same as those we cite.

The escape speed, $v_{\rm esc}$, versus circular speed for each simulated galaxy. Here, the escape speed is set to be the speed of the fastest dark matter particle in the given simulation's solar ring. The blue line shows the power-law model from Equation <ref>, $v_{\rm esc} \sim v_c^{\data{veschat_slope}}$, which has a 0.82 $r^2$ coefficient of determination in log-space. The grey band shows the 1$\sigma$ uncertainty around the prediction. The red point with error bars shows our prediction for the local escape velocity near the Sun in the Milky Way using this model: $v_{\rm esc}(v_\mathrm{c,MW})=\data{vesc_mw(vc)}\pm\data{dvesc_mw(vc)}\,\velunit$.

This figure shows the measured speed distribution data for three example systems in black histograms. The blue lines show the Maxwellian distributions defined by Equation <ref>. These distributions have only one parameter, the peak speed $v_0$. In this figure, we adopt the common assumption that peak speed is equal to the circular speed $v_\mathrm{c}$ measured in each galaxy. This assumption under-predicts the peak speed.
For illustrative purposes, the grey vertical dashed lines show the escape speed, $v_{\rm esc}$, the speed of the fastest dark matter particle in a given simulation's solar ring.
The plots here also show the escape speed based on potential, $v_{\rm esc}(\Phi)$, where we take escape distance to be infinity. Section <ref> provides further detail on the calculation method for $v_{\rm esc}(\Phi)$. Take notice of the difference between these two possible cutoff speeds. If one were to cutoff the distributions at the escape speeds determined by potential, one would assume the presence of too many high-speed particles.
However, even if we were to cut the distributions at the $v_{\rm esc}$ we measure, the Maxwellian would still give too many high-speed particles.
The green lines show best fits for each individual galaxy using the two-parameter function given by Equation <ref>. The green fits to the black data are strikingly tight and demonstrate that this functional form can, in principle, reproduce the speed distributions accurately.

§.§ Escape Speed

As discussed in Section <ref>, we define the escape speed, $v_\mathrm{esc}$, as the smallest speed above which we find no dark matter particles in the simulated local region. Figure <ref> shows the relationship between $v_\mathrm{esc}$ measured this way and the local circular speed. The two parameters correlate as

\begin{equation}
v_{\rm esc}(v_{\rm c})=v_{\rm e,0}\left(\frac{v_\mathrm{c}}{100\,\velunit}\right)^\varepsilon,
\label{eqn:vesc(vc)}
\end{equation}

where $v_{\rm e,0}=\data{veschat_amp}^{+\data{dveschat_amp_plus}}_{-\data{dveschat_amp_minus}}\velunit$ and $\varepsilon=\data{veschat_slope}\pm\data{dveschat_slope}$. Its predictions exhibit a 0.82 coefficient of determination $r^2$ versus the $v_{\rm esc}$ data in log space. Figure <ref> visualizes this fit and shows our $v_\mathrm{esc}(v_\mathrm{c,MW})=\data{vesc_mw(vc)}\pm\data{dvesc_mw(vc)}\,\velunit$ estimate for the Milky Way. As we continue to emphasize, this and the other models in this paper are calibrated to non-observable data but, once fit, require only one single observable feature—the local circular speed—to make predictions for the Milky Way.

§.§ Speed Distribution

A common assumption for the dark matter speed distribution is a simple Maxwellian (Equation <ref>).
Figure <ref> compares the true speed distribution data for three of our galaxies (black) to the Maxwellian assumption (blue), where we have fixed the peak speed, $v_0$, in  Equation <ref> to the local circular speed, $v_c$, as is the standard choice.  We see that the Maxwellian curve peaks at too low of a speed and has an excess high-speed tail. Even if we truncated the Maxwellian distribution at the formal escape speed set by the potential (dashed black line) or the escape speed calculated by our approach (Section <ref>, grey dashed line), the Maxwellian curve would still over-produce counts at high speeds. We find that similar shape mismatches exist for all of our simulated galaxies.

The red lines represent the best-fit speed distributions predicted by this study's final model using only the galaxy's $v_{\rm c}$, specifically using Equations <ref>, <ref>, and <ref>. The yellow bands represent the $1\sigma$ uncertainty around those predictions, determined by drawing from 5000 Gaussian distributions of $D$ and $H$ with means and standard deviations of $D=\data{d}\pm\data{dd}\,\mathrm{km\,s^{-1}}$ and $H=\data{h}\pm\data{dh}\ \mathrm{km\,s^{-1}}$.  By construction, these quoted standard deviations capture 68% of the aggregate black histogram data. The green lines show the speed distributions outside of the 68% band. Section <ref> details the method for simultaneously determining the yellow uncertainty band and the uncertainties on $D$ and $H$.

While the Maxwellian curves shown in Figure <ref> are not best-fit Maxwellian shapes to each galaxy,  we have found that even allowing for a variable $v_0$ for each galaxy, the simple Maxwellian fit fails. As discussed in the introduction, others have found similar results and tried to provide better descriptions of the speed distribution by employing models such as generalized Tsallis distributions <cit.> or other modifications to the Maxwellian distribution <cit.>.

We find a good fit by introducing a speed parameter $v_\mathrm{damp}$, above which each distribution is suppressed:

\begin{equation}
\mathcal{S}(v_\mathrm{damp}-v),
\label{eqn:final_model}
\end{equation}

where $\mathcal{S}$ is the sigmoid function

\begin{equation}
\mathcal{S}(v_\mathrm{damp}-v)=\frac{1}{1+\exp(-k(v_\mathrm{damp}-v))},
\label{eqn:theta}
\end{equation}

and $k$ represents how strongly the sigmoid suppresses the distribution around $v_\mathrm{damp}$. We adopt $k=\data{k}$ for reasons we describe below.
The green lines in Figure <ref> plot this model using the best fit parameter set $(v_0, v_\mathrm{damp})$ for each individual galaxy; it matches the data remarkably well. We find similar success for all of the galaxies in our simulated sample. Note that these fits use the non-observable dark matter speed distributions to directly determine $(v_0, v_\mathrm{damp})$. Such a fit does not yet gauge our ability to define the best pair of parameters to use for the Milky Way itself; however, it shows that, in principle, the model's shape has enough flexibility to reproduce the range of distributions we predict.

Our goal now is to determine the best pair of parameters $(v_0, v_\mathrm{damp})$ to apply to the Milky Way. Our approach is to search for a way to predict each simulated galaxy's values using power law models on their local circular speed:

\begin{equation}
\label{eqn:final_v0}
\end{equation}

\begin{equation}
\label{eqn:final_vdamp}
\end{equation}

We find $D=\data{d}\pm\data{dd}\,\mathrm{km\,s^{-1}}$, $e=\data{e}$, $H=\data{h}\pm\data{dh}\ \mathrm{km\,s^{-1}}$, and $j=\data{j}$ provides a good description for our set of twelve model galaxies. We fit all these parameters, including $k$, by minimizing the aggregate sum of squared errors (SSE) between the model and the data for all simulated galaxies.  In this process we find the best single value of $k$ that provides the minimum SSE for the full analysis, and this results in the adoption of $k=\data{k}$ mentioned earlier.  The individual galaxy fits shown in green in Figure <ref> use the same $k$ as the one that optimizes this universal model.

Figure <ref> shows the speed distribution data for all twelve Milky Way sized FIRE-2 disks (black histograms) along with 5000 realizations each of our model predictions based only on each galaxy's $v_c$ (red, yellow, and green lines). The red lines show best-fit predictions. The yellow bands show $1\sigma$ range speed distributions predicted by the quoted model uncertainties in $H$ and $D$, and the green lines are realizations that fall outside of the $1\sigma$. We generate the green lines by taking Gaussian samples of $D$ and $H$, given their uncertainties quoted above. Section <ref> of the Appendix details the process for determining those uncertainties. Remarkably, this simple model does a good job of predicting the dark matter speed distributions in each simulated galaxy using only $v_c$.

Figure <ref> provides our prediction for the Milky Way's local dark matter speed distribution using $v_c = v_\mathrm{c,MW}=229\pm7\ \mathrm{km\,s^{-1}}$. The implied parameters in Equation
<ref> are $v_0(v_{\rm c, MW})=\data{v0_mw}\pm\data{dv0_mw}\,\velunit$ and $v_{\rm damp}(v_{\rm c, MW})=\data{vdamp_mw}\pm\data{dvdamp_mw}\,\velunit$, and these generate the red line in Figure <ref>. The dashed blue line shows the standard Maxwellian assumption. Notice that the peak speed of our model is higher than in the standard assumption, while the damping behavior of our model pulls the high-speed tail below the Maxwellian starting around $\data{vcrit_mw}\,\velunit$.

Figure <ref> shows the implied halo integral (Equation <ref>) for our predicted speed distribution in the Milky Way along with the same for the standard Maxwellian assumption. Comparing the two in log space (left) shows that the Maxwellian over-predicts the halo integral at the highest velocities.  The dashed blue line cuts the Maxwellian at our best-fit escape speed for the Milky Way and does a better job. The difference at low velocity is more visible in linear space (right), where the Maxwellian over-predicts at low velocity, owing to its lower peak speed.

The red line represents this work's final prediction for the Milky Way's speed distribution. The pink band shows the $1\sigma$ uncertainty range, calculated in the same way as the uncertainty bands in Figure <ref>. The blue dashed line shows a Maxwellian with peak speed equal to the local circular speed and a sharp cutoff at the escape speed. We use Equation <ref> to determine $v_{\rm esc}(v_{\rm c, MW})=\data{vesc_mw(vc)}\,\velunit$. The difference between the blue and red curve points to the benefits of using this work's prediction for the Milky Way's speed distribution over the standard assumption.

Left: Log-space comparison of the halo integral for the Maxwellian model in blue versus this study's model in red. We do not force the solid lines' speed distributions to zero at the escape speed. In this regime, there is significant difference between our model and the Maxwellian. Given that standard practice is to cut the distribution at the escape speed, we also cut both models at $v_{\rm esc}(v_{\rm c, MW})=\data{vesc_mw(vc)}\,\velunit$ and show the results in dashed lines. In that case, the two models are actually similar at high speeds. Right: Linear-space comparison of the halo integrals. At low speeds, the two diverge. The divergence suggests that there is a benefit to adopting this work's more accurate speed distribution model.

§ DISCUSSION AND CONCLUSIONS

We have used FIRE-2 simulations of Milky Way size disk galaxies to inform predictions of the local dark matter density, velocity dispersion, and speed distribution near the Sun using one observable: the circular speed near the Sun, $v_c$.  We then use our best-fit models to inform dark matter expectations in the Milky Way using $v_{\rm c,MW}=229\pm7\,\velunit$, as determined by ref. <cit.>.   Our main results are summarized in Table <ref>. All these parameters are crucial for informing and interpreting  dark matter direct detection experiments.

Local dark matter density: In our simulations, the local dark matter density near solar locations is well characterized by a power law $\rho \propto v_{\rm c}^{\alpha}$, with $\alpha \simeq 0.8$ (Equation <ref> and Figure <ref>).  Assuming that the same model holds for the Solar location in the Milky Way, we find find $\rho(v_{\rm c, MW})=\data{rho_GeV}\pm\data{drho_GeV}\,\denunit$.  This value is similar to those obtained by more complicated kinematic models, and is remarkably close to the average of past estimates since 2018 summarized in Table <ref>. Another relevant comparison is to ref. <cit.>, who find $\data{rho_canon_GeV}\pm\data{drho_canon_GeV}\,\denunit$ using the average of
15 independent $\rho$ determinations from literature published between 2010 and 2020.

Local dark matter velocity dispersion:  We find that the 3D velocity dispersion of dark matter near mock solar locations in our simulations follows a power law given by Equation <ref> (see Figure <ref>).  Note that the best-fit power-law is slightly flatter than linear; $\sigma \propto v_{\rm c}^{\gamma}$, with $\gamma \simeq 0.8$.  This is different[It is not generally expected that the velocity dispersion should be proportional to the circular speed.  Even with an isotropic velocity dispersion ($\beta = 0$), the spherical Jeans equation tells us that $v_c^2 \propto \gamma_n \, \sigma^2$, where $\gamma_n = -\mathrm{d}\ln n/\mathrm{d}\ln r$ is the log-slope of the dark matter profile.   It would not be surprising if $\gamma_n$ varied systematically with $v_c$ such that $\sigma$ and $v_c$ were not linearly proportional. ] than what is expected in an isothermal sphere model, where the 3D velocity dispersion is linear with circular speed as $\sigma_{\rm 3D} = \sqrt{3/2} \, v_{\rm c}$.
Assuming Equation <ref> holds for the Milky Way we find  $\sigma_{\rm 3D}(v_{\rm c,MW})=\data{disp}^{+\data{ddisp_plus}}_{-\data{ddisp_minus}}\,\velunit$. This happens to be very close to $\sqrt{3/2} \,  v_{\rm c,MW}$ for our adopted Milky Way circular speed, though for much larger or smaller circular speeds, our model would give different results.

Local escape speed: Conventional escape-speed estimates depend on an assumption for the escape distance and have recently ranged anywhere from ${\sim}470$ to ${\sim}580\,\rm km\,s^{-1}$ <cit.>. Rather than choosing a distance, we measure the escape speed in solar regions in our simulations as the speed above which no particles are found. With this choice, the escape speed in our simulations is well characterized by a power-law, $v_{\rm esc} \propto v_{\rm c}^\varepsilon$, with $\varepsilon \simeq 0.5$ (Equation <ref> and Figure <ref>).
Assuming this holds for the Solar location in the Milky Way we find
$v_{\rm esc}(v_{\rm c, MW})=\data{vesc_mw(vc)}\pm\data{dvesc_mw(vc)}\,\velunit$.

Speed distribution: We find that distribution of dark matter particle speeds is well-described by a modified Maxwellian (Equation <ref>) with two shape parameters, both of which correlate with the observed $v_{\rm c}$. We use that modified Maxwellian to predict the speed distribution of dark matter near the Sun and find that it peaks at a most probable speed of $\data{v0_mw}\,\mathrm{km\,s^{-1}}$ and begins to truncate sharply above $\data{vcrit_mw}\,\mathrm{km\,s^{-1}}$. (See Figure <ref>.) This peak speed is somewhat higher than expected from the standard halo model, and the truncation occurs well below the formal escape speed to infinity, with fewer very-high-speed particles than is often assumed. The best fit parameters are given in Table <ref>.

Final results for the Milky Way

Parameter     Best estimate     Units

$\rho(v_\mathrm{c,MW})$     $\data{rho_GeV}\pm\data{drho_GeV}$     $\denunit$

$\data{rho_1e7msun}^{+\data{drho_1e7msun_plus}}_{-\data{drho_1e7msun_minus}}$     $10^7\mathrm{M_\odot\,kpc^{-3}}$

$\sigma_\mathrm{3D}(v_{\mathrm{c,MW}})$     $\data{disp}^{+\data{ddisp_plus}}_{-\data{ddisp_minus}}$     $\mathrm{km\,s^{-1}}$

$v_{\rm esc}(v_{\rm c, MW})$
$\data{vesc_mw(vc)}\pm\data{dvesc_mw(vc)}$     $\velunit$

$v_0(v_\mathrm{c,MW})$*     $\data{v0_mw}\pm\data{dv0_mw}$     $\mathrm{km\,s^{-1}}$

$v_\mathrm{damp}(v_\mathrm{c,MW})$*     $\data{vdamp_mw}\pm\data{dvdamp_mw}$     $\mathrm{km\,s^{-1}}$

*Shape parameters for our model of the speed distribution defined in Equation <ref>

This study's major contribution to the local dark matter speed distribution discussion is our methodology for predicting the Milky Way's parameter values, which slide higher or lower based on the assumed local circular speed of the Galaxy.  While we have presented a functional form that works remarkably well, others could apply our methodology to other functions. In fact, we do just that for the Mao parameterization in the Appendix Section <ref>.

One potential subject for further exploration is the effect of Magellanic Clouds and a Sagittarius stream on local dark matter speed distributions. For example, simulations have shown that both the Sagittarius dwarf <cit.> and LMC <cit.> could boost the Milky Way's high-speed tail. There is opportunity to extend an analysis like the one given here to included simulated systems with LMC-size objects and extended streams to more accurately predict the high-speed tail of dark matter particles.

§ ACKNOWLEDGMENTS

PGS and JSB were supported by NSF grant AST-1910965 and NASA grant 80NSSC22K0827. MBK acknowledges support from NSF CAREER award AST-1752913, NSF grants AST-1910346 and AST-2108962, NASA grant 80NSSC22K0827, and HST-AR-15809, HST-GO-15658, HST-GO-15901, HST-GO-15902, HST-AR-16159, HST-GO-16226, HST-GO-16686, HST-AR-17028, and HST-AR-17043 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555.

§ FITTING PROCEDURE AND DEFINING UNCERTAINTIES

This study develops models for predicting local dark matter density, velocity dispersion, and escape speed based on circular-speed power laws. It achieves this by performing linear regressions on the logarithms of both the target variable and circular speed. Each model has a band of uncertainty around its prediction curve. The band has two components: the standard error of the regression and the standard error of the mean. The former is

\begin{equation}
s_{\rm reg}=\sqrt{\frac{1}{N-2}\sum_i^N e_i^2}
\end{equation}

where $N$ is the number of data points, and $e_i$ is the difference between data point $i$ and the model's prediction for that point. The latter captures the increasing uncertainty in the prediction with increasing distance from the average circular speed. It is

\begin{equation}
s_{\rm mean}(\log v_{\rm c})=\frac{s_{\rm reg}}{\sqrt{N}}\sqrt{1+\frac{(\log v_{\rm c}-\overline{\log v}_{\rm{c}})^2}{{\rm var}(
\log\boldsymbol{v_{\rm c}})}},
\end{equation}

\begin{equation}
{\rm var}(\log\boldsymbol{v_{\rm c}})=\frac{\sum_i^N(\log v_{{\rm c},i}-\overline{\log v}_{\rm c})^2}{N}.
\end{equation}

The standard error of a given prediction is then $s_{\rm pred}(\log v_{\rm c})=\sqrt{s_{\rm reg}^2+s_{\rm mean}^2(\log v_{\rm c})}$. From here, the final uncertainty band around a given prediction is

\begin{align}
\delta_{\rm log\,band}(\log v_{\rm c})=&t_{\rm crit}s_{\rm pred}(\log v_{\rm c})\\
\approx&0.42\,s_{\rm pred}(\log v_{\rm c})
\end{align}

where $t_{\rm crit}$ is the critical two-tailed t-value assuming 68% confidence and $N-2$ degrees of freedom. Keep in mind that $\delta_{\rm log\,band}(\log v_{\rm c})$ applies to the prediction line in log space. The transformation from the logarithmic error band to linear is

\begin{equation}
\delta_{\rm band}(v_{\rm c})
= \pm 10^{\log y \pm \delta_{\rm log\,band}(v_{\rm c})} \mp y.
\end{equation}

For our Milky Way prediction, it is also necessary to propagate the $\delta v_{\rm c,MW}=7\,\velunit$ systematic uncertainty in circular speed from ref. <cit.> and the $\delta
\rho_{\rm data}/\rho_{\rm data}=\data{stdev_linear_dendiff}$, $\delta\sigma_{\rm data}/\sigma_{\rm data}=\data{stdev_linear_dispdiff}$ uncertainties in this study's density and dispersion data points, which stem from possible variations across the Solar ring discussed in Section <ref>. Therefore, the final uncertainty is

\begin{equation}
\delta y_{\rm MW} = \sqrt{\delta_{\rm band}^2(v_{\rm c}) + \delta v_{\rm c,MW}^2\left(\diff{y}{v_{\rm c}}\right)^2 + y^2\left(\frac{\delta y_{\rm data}}{y_{\rm data}}\right)^2}
\end{equation}

where $y$ is one of $\rho$, $\sigma_{\rm 3D}$, or $v_{\rm esc}$. Note that we do not evaluate intra-ring variations in escape speed, so we effectively assume $\delta v_{\rm esc}/v_{\rm esc}=0$.

For the speed distribution model, this study performs a least-squares fit to the aggregate of all twelve galaxies. The feature matrix consists of columns 1) dark matter particle speed $v$ and 2) local circular speed $v_\mathrm{c}=\langle\vec{v}\cdot\hat{\phi}\rangle_{T\leq 10^4\mathrm{K}}$ discussed in Section <ref>. Each galaxy contributes 50 rows to the feature matrix and probability-density target vector. These correspond to each galaxy's 50-bin histogram of dark-matter-particle speeds. For a given galaxy's entries into the circular-speed column of the feature matrix, the same $v_\mathrm{c}$ value is repeated 50 times. The 50 data rows from each galaxy yield an aggregate of 600 data points.

We quantify this model's uncertainty using the $D$ and $H$ parameters in Equations <ref> and <ref> with the following process.  Our goal is to determine values of $\delta D$ and $\delta H$ that, when sampled, deliver speed distributions that capture 68% of the speed histogram bins (50 for each galaxy). We sample 5000 pairs of $(D, H)$ values from normal distributions with means equal to the parameters' best-fit values and standard deviations equal to their uncertainties $\delta D$ and $\delta H$. Keeping all other parameters constant at their best-fit values, those $(D, H)$ pairs yield 5000 speed distributions. For each speed $v$, the bottom and top of the yellow $1\sigma$ uncertainty band in Figure <ref> correspond to the 16th and 84th percentiles of all predicted histogram values. The key step in the uncertainty analysis is identifying the $(\delta D, \delta H)$ set that causes the yellow band to capture 68% of the black histogram data. We accomplish this by generating a $26{\times}26$ meshgrid of $\delta D/D$ and $\delta H/H$ values each ranging from 0.00 to 0.16, evaluating the percentage of data captured by the $1\sigma$ band for each point on the grid, and identifying the $(\delta D/D,\, \delta H/H)$ point with the smallest percentage-captured value over 68%. This process yields $(\delta D/D,\, \delta H/H)=(0.096, 0.058)$, corresponding to $(\delta D,\, \delta H)=(\data{dd},\data{dh})\rm\ km\,s^{-1}$.

§ VELOCITY DISPERSION ANISOTROPY

Figure <ref> summarizes the anisotropy of the local dark matter velocity dispersions in cylindrical coordinates $\phi$, $z$, and $r$.
On the left we show the anisotropy parameters, defined as

\begin{align}
&\beta = 1 - \dfrac{\sigma_\phi^2+\sigma_z^2}{2\sigma_r^2} \label{eqn:beta}\\
&\beta_\phi = 1 - \dfrac{\sigma_\phi^2}{\sigma_r^2} \\
&\beta_z = 1 - \dfrac{\sigma_z^2}{\sigma_r^2}. \label{eqn:betaz}
\end{align}

Positive values correspond to particles moving faster radially than vertically and azimuthally, while negative values convey the opposite. Overall anisotropy is low for FIRE disks; the majority have $\beta<0.2$, and only Louise is higher than 0.3. Looking at specific components, azimuthal and radial dispersions tend to match more closely than vertical and radial, as represented by the lower cyan and higher orange bars. The right panel directly shows the ratio of the velocity dispersion components. For most galaxies, the radial dispersion is slightly larger than the other two components, and the $z$ dispersion is the smallest.  One interesting outlier is Romulus, which has a slightly larger $\phi$ velocity dispersion than radial velocity dispersion. Though some level of velocity dispersion anisotropy does exist, the ratios are all close to unity at ${\sim}0.8$–1.

Characterization of the anisotropy of dark matter velocity dispersion over our suite of simulations. The left panel shows the anisotropy parameters in cylindrical coordinates, defined in Equations <ref> – <ref>. On the right we show the ratio of velocity dispersion components. In most cases the radial dispersion is slightly larger than the other two components, and $\sigma_{\phi}$ is larger than $\sigma_{z}$. In only one case (Thelma) is the $z$-component of dispersion larger than the $\phi$ component. Overall, anisotropy is low, with $\beta<0.2$ for the majority of galaxies in the left panel. Additionally, dispersion ratios in the right panel are all approximately between 0.8 and 1.

§ USING THE MAO SPEED DISTRIBUTION

Note that, in the past, studies have evaluated the best-fit Mao <cit.> shape parameters (Equation <ref>) for individual simulated galaxies (e.g. <cit.>), but have not gone beyond implying that the best parameters for the Milky Way lie somewhere in the range of parameter values found for those simulations. In this work have aimed to provide models that can be applied to the Milky Way using the observed value of $v_c$.  In this section we do so for the Mao parameterizaton.

The red curves show the final model that this study recommends, given by Equations <ref> – <ref>. The grey curve shows a straightforward implementation of the Mao parameterization with the peak speed equal to the “observed" circular speed of each galaxy. Specifically, it shows Equation <ref> with parameters $(v_0, v_{\rm esc}, P)=\left(v_{\rm c},\, v_{\rm esc}(v_{\rm c}),\, \data{p_mao_naive_agg}\right)$, where Equation <ref> provides $v_{\rm esc}(v_{\rm c})$ and $P=\data{p_mao_naive_agg}$ minimizes the SSE between the model and the data for the aggregate of all twelve galaxies. The cyan curve is an implementation of the Mao parameterization that allows both $v_{\rm esc}$ and $v_0$ to vary using best-fit power laws on $v_{\rm c}$. Using Equation <ref> to determine $v_{\rm esc}$ is the same as in the standard implementation; however, we also predict $v_0$ with Equation <ref> as opposed to the standard $v_0=v_{\rm c}$ assumption. To be precise, our full Mao implementation uses $(v_0, v_{\rm esc}, P)=\left(v'_0(v_{\rm c}),\, v_{\rm esc}(v_{\rm c}),\, \data{p_mao_lim_fit}\right)$, where $P=\data{p_mao_lim_fit}$ minimizes the aggregate SSE. Each panel shows the root mean square errors of the three models versus the black histogram data, in units of $10^{-3}\,\rm km^{-1}\,s$. The $v_0 = v_c$ version of the Mao model performs surprisingly well, but our model offers an improvement; quantitatively, it is the most accurate of the three models shown. Nevertheless, if one strongly prefers Mao's functional form, one can achieve comparable accuracy by using our Equations <ref>, <ref> and <ref> to determine the parameter values for the Mao function. The root mean squared error for the aggregate of all twelve galaxies is $\rm{RMS}_{\ref{eqn:final_model}}=0.25\times10^{-3}$ for Equation <ref>, $\rm{RMS}_{\rm Mao}=0.27\times10^{-3}$ for Mao with our method, and $\rm{RMS}_{\rm{M}v_{\rm c}}=0.33\times10^{-3}\ \rm km^{-1}\,s$ for Mao with $v_0=v_{\rm c}$.

Figure <ref> shows the dark matter speed distributions for each of our galaxies (black histograms) compared to best-fit models using the Mao parameterization defined by Equation <ref>, along with our preferred model. We first evaluate the predictive strength of the Mao parameterization by using $v_0=v_\mathrm{c}$ and escape speed given by Equation <ref>.  For the third parameter, P, we make the simplifying assumption that one value should work well for all galaxies. Using $P=\data{p_mao_naive_agg}$ minimizes the sum of squared errors of the aggregate of all twelve disks' predicted distributions. This version of Mao performs well. Even so, our model given by Equations <ref> – <ref> does a little better in most cases. It provides a 25% improvement from Mao with $v_0=v_\mathrm{c}$ and $P=\data{p_mao_naive_agg}$ judged by the root mean squared error of the aggregate of the predicted distributions, pushing it from $\rm{RMS_{Mao}}=0.33\times10^{-3}\,\rm km^{-1}\,s$ down to $\rm{RMS_{5.7}}=0.25\times10^{-3}\,\rm km^{-1}\,s$.

Furthermore, this study can improve the Mao implementation by predicting the best value for $v_0$ with a methodology similar to our model. In this, which we call our implementation of the Mao parameterization,

\begin{equation}
\label{eqn:mao_v0}
\end{equation}

$D'=\data{d_mao_lim_fit}\pm\data{dd_mao_lim_fit}\,\velunit$, and $e'=\data{e_mao_lim_fit}$. Additionally, the best fit value for $P$ is

\begin{equation}
\label{eqn:p_ours}
\end{equation}

This framework performs very well, as shown by the cyan line in Figure <ref>, although our preferred model tends to perform slightly better; the root mean square error for our implementation of Mao across the aggregate of all twelve galaxies is $\mathrm{RMS}_{\mathrm{M}v_\mathrm {c}}=0.27\times10^{-3}\,\mathrm{km^{-1}\,s}$. Although we do prefer the model described by Equations <ref> – <ref> given its higher performance, the results here show that if one's preference is to use Mao, one could achieve almost as accurate results by implementing our methods of determining $v_{\rm esc}(v_{\rm c})$ with Equation <ref> and predicting $v_0$ for Mao using Equation <ref>.

§ CIRCULAR SPEED METRIC

For completeness we provide Figure <ref>, which plots this study's observable proxy for circular speed, $v_\mathrm{c}=\langle\vec{v}\cdot\hat{\phi}\rangle_{T\leq 10^4\mathrm{K}}$, versus the spherical ideal $\sqrt{GM/R_0}$.
The two exhibit a near one-to-one relationship, with only a 4% scatter. Therefore, this study's results are likely insensitive to the choice of whether to use the ideal or observable circular speed.

Here we compare measured values of $v_\mathrm{c}$ (the average $\phi$ velocity of cold gas $v_\mathrm{c}=\langle\vec{v}\cdot\hat{\phi}\rangle_{T\leq 10^4\mathrm{K}}$) in our simulations to a spherical idealization of the circular speed $\sqrt{GM/R_0}$. The one-to-one line on this scatter plot has a 0.96 $r^2$ coefficient of determination. This suggests that the detailed choice of how we measure or characterize $v_\mathrm{c}$ in the simulations is probably not crucial to the success of our power-law characterizations in the paper.

§ HALO INTEGRALS

The red lines represent the halo integrals for this work's speed-distribution model given by Equations <ref> through <ref>. The solid blue lines represent that given by the simple Maxwellian without any truncation. Standard practice is to cut the distribution at the escape speed, so we also cut both models at the $v_{\rm esc}(v_{\rm c})$ given by Equation <ref> and show the results in dashed lines.

Figure <ref> compares halo integral performance for our model in red verses a Maxwellian in blue. For both models, the solid lines show standard non-truncated versions while the dashed lines truncate the speed distributions sharply at $v_{\rm esc}(v_{\rm c})$. Even when truncating the Maxwellian, its halo integral tends to exhibit an excessive bump at high speeds while our model's integral provides better agreement with the data.

sofue2020rotation, necib2019under