A phenomenological model for dark matter phase space distribution

Understanding the nature of dark matter is among the top priorities of modern physics. However, due to its inertness, detecting and studying it directly in terrestrial experiments is extremely challenging. Numerical N-body simulations currently represent the best approach for studying the particle properties and phase space distribution, assuming the collisionless nature of dark matter. These simulations also address the lack of a satisfactory theory for predicting the universal properties of dark matter halos, including the density profile and velocity distribution. In this work, we propose a new phenomenological model for the dark matter phase space distribution. This model aims to provide an NFW-like density profile, velocity magnitude distribution, and velocity component distributions that align closely with simulation data. Our model is relevant both for theoretical modeling of dark matter distributions, as well as for underground detector experiments that rely on the dark matter velocity distribution for experimental analysis.


INTRODUCTION
Dark matter is an invisible and mysterious substance that makes up a great portion of the universe.While its existence is solely inferred from its gravitational effects, its true nature and composition remain one of the most significant questions in modern astrophysics and particle physics.Detecting dark matter has proven to be incredibly challenging because it does not emit or interact with light.Based on these limitations, cosmological N-body simulations may, at the moment, be the only possibly way to study dark matter structures and distributions.These simulations are performed assuming the collisionless property of dark matter, see the review of simulations Vogelsberger et al. (2020); Angulo & Hahn (2022).
Another important aspect is the velocity distributions of dark matter particles.The direct detection experiments of dark matter, which aim to detect the low energy nuclear recoil from rare scattering events between target nuclei and dark matter particles, relies on the knowledge of velocity distributions near the location of solar system.Therefore, it is essential to know the velocity distributions of dark matter halos.The inverse-Eddington method is very efficient to perform this task Eddington (1916), however, it is only valid under certain limiting circumstances Lacroix et al. (2018).There are also degeneracies in the velocity distributions for a given density profile when the system is not ergodic, because single variable function () is not sufficient to determine the velocity distribution which involves multiple variables.An ideal velocity distribution function should also provide the distribution of velocity components.As an example, the standard halo model (SHM) is known to not provide the right description on the velocity distributions Vogelsberger. et al. (2009); Kuhlen. et al. (2010); Ling et al. (2010); Fairbairn & Schwetz (2009b), especially the behavior of high velocity tails.Therefore, many velocity distribution models have been proposed, either from the first principle or theories Hansen et al. (2005) (2015), guessing from the empirical fitting King (1966); Lisanti et al. (2011); Mao et al. (2013); Callingham et al. (2020); Bozorgnia et al. (2016), inferring from the observational data Fairbairn & Schwetz (2009a); Bhattacharjee et al. (2013); Mandal et al. (2019), or parameterizing velocity distributions Kavanagh & Green (2013); Lee (2014); Kavanagh (2015).However, none of these are fully satisfactory and widely accepted.Furthermore, these velocity distribution models also have degeneracy in the velocity components since they all use the velocity magnitude rather than the velocity components as principle variables.
Instead of working on the density profile and velocity distribution separately, we suggest applying an approach where the density profile (as well as other universal properties) and velocity distributions are consistent with each other.One promising way is to propose a phase space distribution function from which the density profile and velocity distributions could easily be derived and thus automatically be consistent with each other.The density could be obtained through the velocity integral of the phase space distribution function, and the velocity distribution could also be derived by fixing the radius and potential in the phase space distribution function Binney & Tremaine (2008).In previous works, there are action-based phase space distribution that were derived using the action-angel method Posti et al. (2015); Piffl et al. (2015); Binney & Vasiliev (2023).However, the reasoning employed to create the distribution function does not apply to the NFW case.Consequently, there is no action distribution function for NFW density profile, instead they proposed an empirical model for NFW Posti et al. (2015).In addition, the action-based models generally do not provide analytical formula for the derived quantities, making their application more challenging.In this work, we propose a new analytical phase space distribution, which could result in a NFW-like density profile, the radial varied anisotropy as well as power law phase-space density (although certain disagreements with simulation data exist).What's more, it can also provide us with the velocity magnitude distribution as well as its components, with natural cut off at the local escape velocity.We will also compare our analytical velocity distributions with simulation data, which has radial and tangential velocity data separately allowing to break the degeneracy of the velocity magnitude distributions.Our model fit quite well to the radial velocity data and also the low velocity regime of tangential velocity data, giving consistent results on the distributions as well as the fitting parameters within an acceptable range.This provides support that our analytical model may be relevant for understanding the gravitational dynamics of cosmological dark matter structures.This paper is organized as follows: In Section.2, we will present our dark matter phase space distribution function and study its predictions on the density profile, anisotropy parameter, as well as the phase space density.Additionally, we will provide the formulas for the velocity magnitude , the radial and tangential velocity components distributions.Next, in Section.3, we will compare our radial and tangential velocity distributions with simulation data.The velocity data were extracted from different radial bins (shell) of equilibrium dark matter halos, with both radial and tangential velocity data available in each bin.For comparison, we use data from two different simulation schemes, namely the cold collapse and explicit energy exchange.The fitting results and estimated parameters are also presented in this section.Finally, we will conclude and discuss our findings in Section.4.

THE DISTRIBUTION FUNCTION AND ITS PROPERTIES
Based on the observations of the simulation data introduced in Section.3, a reasonable distribution function should be spherically symmetric due to the simulated halos have no preferred spatial directions.Therefore, it is convenient to write out the distribution function of dark matter halos in spherical coordinates (, , ).The distribution function should also be anisotropic because the simulation data of radial and tangential velocity distributions are different.That means we need take into account the angular momentum of dark matter particles.It would be a good choice to put the angular momentum in the exponent because when we integrate the distribution function to obtain the density profile, it will result in a term of 1 +  2 in the denominator.Then if we wish to obtain an NFW-like profile, we could just multiply it by the radial distance .Meanwhile, we would like the distribution function to have a natural cutoff, and it could be realized by multiplying the binding energy.After considering the above conditions, our distribution function is given by where  0 is the parameter to normalize radial distance, and  is the parameter to normalize the binding energy  and exponent  which are given by where  is the positive potential,   and   = √︃  2  +  2  are respectively the radial and tangent velocity, and  is the angular momentum of dark matter particles.One could notice that this distribution function does not satisfy the Jeans theorem since there is a factor  0 /, which means that the full form is not only a function of the integrals of motion.This could be seen as the limitation of our model.However, this limitation does not undermine the primary objective of our work.We have proposed a phenomenological or effective model rather than presenting a comprehensive theory.Therefore, it is acceptable for breaching the Jeans theorem as long as it is consistent with the simulation data.
For simplicity, in the following discussion of this section, we will denote  = / 0 and set  = 1.Then the potential, binding energy and  are all in unit of  2 , the velocity is in unit of .The density of the above distribution is given by the integral with respect to velocities where   is the normalization constant, which is usually the inverse of the halo mass,   () and   () are called error function and imaginary error function.To obtain (), we first need to know the dependence of  on , which requires solving the spherically symmetric Poisson equation, however, it has to be solved numerically.Starting from  = 0 for some chosen value of (0) and    = 0, we thus can solve for ().In practice, we need to apply a softening to the centre point  = 0 due to the divergent behavior of 1/ of the Poisson equation.If we assume that the central potential is nearly flat within a small radial distance, let's say  = 0.01.Then, we could solve the Poisson equation by choosing that the initial conditions at  = 0.01 is the same as the centre point.Thus, the numerical results are not sensitive to the choice of  = 0.01 as long as the approximation of flat central potential is valid.In the following discussions, we mainly choose to start the numerical solution from  = 0.01 with different initial potentials, i.e., (0.01)= 5, 7, 10 and    | =0.01 = 0.As a comparison, we also show the numerical results of potential when we choose to start with  = 0.001.The gravitational constant  in the numerical solutions is set to be 4.301 × 10 −6   −1 ⊙ (/) 2 , and the normalization factor   is set to   = 1.

potential, density, anisotropy and pseudo phase pace density
The numerical results of the potential are shown in Fig. 1 and Fig. 2 for given initial potentials.We can see that all the potentials with different initial conditions are quite smooth when close to the center and at a significant distance from the center, with a drastic transition between them.A larger initial potential, a greater difference between the central and distant regions, and a faster drop in the transition region.By comparing the potentials in Fig. 1 and Fig. 2, we can see that results are not sensitive to the initial value  in the numerical solutions.The potentials are almost the same regardless of the different initial values of .Therefore, we could safely use one of the numerical results of potential (the case with  = 0.01) to continue our discussion.
With the results of potential, we can obtain the density (6) corresponding to the distribution function (1).The results are shown in Fig. 3.For comparison, we also plot the NFW density profile by setting the characteristic radius of NFW as   = √ 2 0 .All the density profiles are normalized by demanding that they have the same total mass   = 1 within their  200 radius where the average density within this radius are 200 times than the mean density of the universe which is set to be 5×10 −5 for simplicity.We can see that they are very similar to the well-known NFW density profile, which has √ 2  0 .All profiles are normalized by ensuring that they share the same total mass within their  200 radius as   = 1.The concentrations for the curves with  (0.01) = 5, 7, 10 and NFW are 222.7,645.5, 14866.7 and 231.9 respectively.logarithmic density slope  = / approximately −1 near the center and −3 in the outer region of the halo, referring to Fig. 4 for the variation of  with radial distances.However, there are also differences between NFW and our model, especially in the transition region where  is from −1 to −3.Our model seems more sharper than the NFW density profile around the transition region.The initial potentials will affect the logarithmic slope near the center region.The concentrations  ≡  200 / −2 ( −2 is the radius at which the logarithmic slope is −2) for the curves in Fig. 3 with (0.01)= 5, 7, 10 and NFW are respectively 222.7, 645.5, 14866.7 and 231.9.It is evident that larger initial potential values result in higher concentrations.To closely resemble the NFW profile, the initial potential should be sufficiently small.The smaller the initial potential, the more similar the profile to the NFW one.It is remarkable that we can obtain the NFW-like density profile from a phase space distribution function.
Since our model considers the anisotropic case, it is also important to investigate the anisotropy parameter  which is defined as where   and   are respectively the tangential and radial velocity dispersion.They could be compute through We show  as a function of  in Fig. 5.Despite the different initial potentials, we observe that the values of  are nearly identical.The initial conditions do not affect the anisotropy. starts at  = 0, i.e.,  2  /2 =  2  =  2  =  2  , indicating isotropy around the center.As the radius grows,  increase rapidly around the scale radius  ≈ √ 2, reaching  = 1 at larger radial distances.This implies  2  ≪  2  .This type of  behaves similarly to the commonly-used Osipkov-Merritt model Osipkov (1979); Merritt (1985), which is expressed as with   as the scale parameter in Osipkov-Merritt model.We also illustrate the behavior of () when  2  = 2 √ 2  2 0 in Fig. 5.We can see that they are almost the same to our model across a large range of radial distances.This could be attributed to the exponent of our distribution function  which is similar to the Osipkov-Merritt model.() is only related to the scale parameter   and independent of the other parameters of the Osipkov-Merritt model.A similar observation holds for our model, possibly explaining why all the () values in our model are almost identical.Eq.( 12) could serve as an approximation for () in our model, achieved by setting  2  = 2 √ 2  2 0 .We also investigate the pseudo phase space density which is another important empirical laws.In Fig. 6, we plot the pseudo phase space density as a function of .We can see that it has a pretty much the same logarithmic density slope as the density, see Fig.  function of radius Dehnen & McLaughlin (2005) in Fig. 7.The slopes in the outer region are the same as the pseudo phase space density, while there is an overall shift to higher values in the center region.

velocity distributions
The velocity magnitude distribution  () could be obtained by transformation and integrating over the distribution function (1).We get where   is the normalization constant, DawsonF() is the so-called Dawson integral.We have plotted the velocity magnitude distribution for different radii  = 1, 5, 10 in Fig. 8, by choosing the initial potential (0.01)= 5.As comparison, we also show the SHM velocity distributions which corresponds to the Maxwell-Boltzmann phase space distribution.We can see that our velocity distribution is suppressed more than the SHM model in the high velocity tails.The velocity distributions are more concentrated at lower values as the radial distance becomes larger.We also compare the velocity distribution function ( 14) with simulation data from the Aquarius Project Springel et al. (2008);Vogelsberger. et al. (2009).We fitted the data with our model as well as the SHM model in Fig. 9.The velocity distribution data are taken from many 2   boxes located between 7 and 9   from the center of halo Aq-A-1 in Aquarius Project.In each 2   box, there are about 10 4 to 10 5 particles.In Fig. 9, the black line is the median value of velocity distribution measured over all 2   boxes, and the green band encloses 95% of the measured velocity distributions for each velocity bin, and each velocity bin has a width of 5 /.We can see that our model almost overlaps the SHM model, with the exception that our model exhibits a relatively larger suppression in the high velocity tails, resulting in a better fit to the simulation data.
The above velocity magnitude has some extent degeneracy over the velocity components.Therefore, it is crucial to understand the radial and tangential velocity distributions and validate them through simulations.to the velocity distribution data extracted from 2   boxes between 7 and 9   from the center of halo Aq-A-1 in Aquarius Project.The black line is the median value while the green band enclose 95% of the velocity distributions over all 2   boxes.
By integrating over the tangent velocity component in the distribution function (1), this yields the radial velocity distribution.Similarly, we can integrate over the radial velocity component in the distribution (1) to obtain the tangential velocity distribution.They are respectively given by where   and   represent the normalization constants.To fit the real simulation data in the next section, we have kept  as a free parameter rather than  = 1 in the above formulas of  (  ) and  (  ).It can be verified that the radial and tangential velocity distribution functions approach zero at the local escape velocity   =   =    = √︁ 2, indicating a cutoff at the local escape velocity.

the simulation data
We use the velocity data extracted from the simulation in Sparre & Hansen (2012).These numerical simulations were described in detail previously Hansen & Sparre (2012); Eilersen et al. (2017).To test our model, we compare the velocity distribution function with two different scheme of simulation performed in Sparre & Hansen (2012).
The first is called the cold collapse simulation.This aims to simulate the violent relaxing process that occurs in the early universe when structure collapsed.They create a main halo with a number of compact and condensed substructures.All particles started with zero velocities and evolve under gravity.Once it has attained equilibrium, we divide a halo structure into radial bins (thin spherical shell) and thus we can extract the radial and tangential velocity components of the particles from selected bins.We sample the same number of particles in three radial bins whose radial and the tangential velocity data are shown in Fig. 10.To facilitate the readability of the bin data, some of the bin data have been shifted vertically.The three radial bins were chosen near the slope of  0 = / = −1.6,−2.0, −2.4.
The second set of numerically simulated data is the so-called explicit energy exchange simulation.Various types of energy exchange take place between collisionless particles, particularly through violent relaxation and dynamical friction.Therefore, the numerical set-up take into account a perturbation in which the spherical symmetry is preserved, but they permit the particles to exchange energy.Each radial bin is designed to preserve energy with instantaneous energy exchange, ensuring that the perturbation itself has no impact on the density or dispersion profiles.Afterward, the system evolves with normal collisionless dynamics.Similarly, we divide the equilibrium structure into radial bins and then extract the radial and the tangential velocity data.We plot the velocity data in Fig. 11.The three radial bins were chosen near the slope of  0 = −1.7,−2.4 − 3.0.

fitting and parameter estimation for the velocity data
When fitting our radial and tangential velocity distribution model to the real simulation data, we adopt error-bars on these data.We assume that the horizontal (velocity) error for each data point corresponds to the width of velocity bin Δ, while the vertical (frequency) errors are taken as the square root of the counts in each velocity bin.This assumption is based on considering the particles in each velocity bin follow a Poisson distribution.In total, we have four parameters to fit for radial or tangential velocity distribution data: the normalization factor   (  ), the radial distance , the potential (), the parameter .
From an ideal perspective, if our model is the truth, when fitting the radial and tangential data independently in each bin, there should be two constraints for the fitting parameters of the data: first, since the three radial bins data are taken from one halo structure, the value of  for the three radial bins should be the same; second, for each radial bin, the fitting parameters (, ) for the radial and tangential data should also be the same.Instead of fitting the radial and tangential data independently, we choose to jointly fit the radial and tangential data for each bin by using the same values of  and () in the fitting algorithm, i.e. we define the total logarithmic likelihood function as the sum of the radial and tangential parts L (rad and tan data|  ,   , , (), ) = L 1 (rad data|  , , (), ) + L 2 (tan data|  , , (), ) where ran and tan data represent the radial and tangential velocity data in each radial bin.With this joint fitting, the second constraint could be easily implemented.Then, we employ Bayesian inference with flat priors and MCMC sampling method to explore the posterior distributions of parameters given the velocity data from each radial bin.We use 300 random walkers, and each walker takes 5000 steps.
After excluding 10% of the samples in the burning phase, we finally obtain the samples that can be used to analyze our model.The fitting parameters for each bin in the collapse simulation data are presented in Table .1.We show the median values as well as the upper and lower bound that enclose 68% of the MCMC samples..1 that the fitting potential is approximately equal to the local escape energy.This supports the self-consistency of the fitting () in different radial bins.It is noteworthy that our model appears to violate the first constraint mentioned above, as the fitting  values are not uniform.However, in practice, due to the size of radial bins and the limitation of particle numbers in the simulation, some flexibility in adhering to these constraints is expected.
To be more explicit, we present the fitting results to the three radial bins of collapse simulation in Fig. 12. From top to bottom, we present the data and fits for the bins with  0 = −1.6,−2.0, −2.4.In each plot, we show both the radial and tangential velocity data along with the corresponding error-bars, as mentioned above.The solid blue curve represents the radial velocity distribution, while the solid red curve represents the tangential velocity distribution, both predicted by the median values of the MCMC samples.Additionally, we include the 1-sigma band of the velocity distribution predicted by the MCMC samples.The band is colored in shallow blue for the radial fit and shallow red for the tangential fit, though it may be relatively small for some bins.The relative residuals, defined as the difference between the data and the solid curves, are calculated as data − solid curve data (18) are also shown in the bottom of each subplot in Fig. 12.The blue and red dots respectively represent the radial and tangential relative residuals.We can see from Fig. 12, our model gives a good fit to all three bins especially for the radial data, despite some residuals remaining.The radial fits of our model are much better than the tangential fits.Our model only performs well in the low velocity regime of tangential data and underestimate their high velocity tails.The relative residuals increase as  0 and velocity increase.Nevertheless, as expected, the radial and tangential velocity converge to each other at higher velocities.The smaller the value of  0 , the higher the similarity of radial and tangential velocity distributions.
We have also shown the fitting parameters of the three radial bins for the explicit energy exchange simulation data in Table .2.Similar to the situation with the collapse simulation data, we also show the upper and lower bounds that encompass 68% of the MCMC the radial and tangential velocity distributions in analytical forms (15) (16).To assess the utility of our velocity distributions, we have respectively compared with radial and tangential velocity data which were extracted from different radial bins of equilibrium dark matter halos.We verify our results in two different simulation schemes, the cold collapse and explicit energy exchange which effectively cover a wide range of possible equilibrium processes.The fits show good performance and the results are promising, especially for the radial data (across the whole velocity range) and the lower velocity regime of tangential data.There is agreement in estimating the relevant parameters of the radial and tangential velocity data.Despite the advantages of providing good fits to the velocity distributions, there are also some potentially important limitations to our model that we wish to address.Our model faces the first challenge in predicting the pseudo phase space density, which does not follow a single, unbroken power law.The second limitation arises when comparing with the data, as significant deviations are observed, especially in the high-velocity tails of tangential data.Additionally, the predicted parameter  differs among three bins of the same halo.These observations could indicate that our model is missing some additional features of dark matter halos.However, it is essential to verify this issue using another independent simulation program since the error could originate from the data rather than the model.If these limitations have been confirmed by other simulations, we should explore possible ways to improve our model.Our model only uses the linear order of binding energy  and exponent , which represents the simplest case.One could consider using   or   in the exponent of (1) with  ≥ 0. Introducing non-linearity could potentially compensate for these limitations.Furthermore, the  0 / factor in our phase space distribution (1) could also be a simplification of some complicated function involving integrated quantities such as  and .If  0 / could be replaced by functions of integrated quantities, the model would be much plausible, indicating its adherence to the Jeans theorem and transcending a mere phenomenological model.

Figure 4 .
Figure 4.The logarithmic slope of the dark matter density  as a function of radial distance , with three different initial conditions  (0.01) = 5, 7, 10, and the NFW density slope by setting   = √ 2  0 .

Figure 8 .
Figure8.The normalized velocity magnitude distribution  2  () as a function of velocity , at three different radial locations  = 1, 5, 10, with initial potential  (0.01) = 5.We also show the Maxwell-Boltzmann distribution which is also known as the standard halo mode (denote as SHM) at  = 1, as comparison.

Figure 9 .
Figure9.The best fit of our model(red dashed) and SHM model(blue dashed) to the velocity distribution data extracted from 2   boxes between 7 and 9   from the center of halo Aq-A-1 in Aquarius Project.The black line is the median value while the green band enclose 95% of the velocity distributions over all 2   boxes.

Figure 10 .
Figure10.The three radial bins data for radial and tangential velocities in the cold collapse simulation.The three radial bins were chosen near the slope of  0 = −1.6,−2.0, −2.4 (From top to bottom), and the  0 = −1.6,−2.0 data were shifted vertically for easy reading.Each data point represents the number of particles within a velocity bin Δ = 0.1.

Figure 11 .
Figure11.The three radial bins data for radial and tangential velocities in the explicit energy exchange simulation.The three radial bins were chosen near the slope of  0 = −1.7,−2.4 − 3.0 (From top to bottom), and the  0 = −1.7,−2.4 data were shifted vertically for easy reading.Each data point represents the number of particles within a velocity bin Δ = 0.07.

Figure 12 .
Figure12.The radial and tangential velocity fitting results to the three radial bins in the cold collapse simulation.From top to bottom are respectively  0 = −1.6,−2.0, −2.4.The radial data and the median radial distributions are colored blue, while the tangential results are colored red.The 1-sigma band of radial distributions are colored shallow blue, and tangential bands are colored shallow red.The relative residuals are also given for the radial (colored blue) and tangential data (colored red).

Table 1 .
The best estimated parameters for the three radial bins in the collapse simulation.We show the median values with upper and lower bounds enclosing 68% of the MCMC samples.Additionally, we display the local escape energy  2  /2 for each bin as reference to compare with the potential.The normalization factors increase with  0 , showing a positive correlation.The median values of radial distance  and potential () are reasonable, as they exhibit growth and decrease respectively with increasing  0 . Aditionally, for comparison with the fitting potential, we also show the local escape kinetic energy  2  /2 for each bin, where   represents the maximum velocity (≈ local escape velocity) of each bin.Ideally, the fitting potential should be equal to local escape kinetic energy, represented by  ≈  2  /2.We could see from Table