Dark Matter in 3D

We discuss the relevance of directional detection experiments in the post-discovery era and propose a method to extract the local dark matter phase space distribution from directional data. The first feature of this method is a parameterization of the dark matter distribution function in terms of integrals of motion, which can be analytically extended to infer properties of the global distribution if certain equilibrium conditions hold. The second feature of our method is a decomposition of the distribution function in moments of a model independent basis, with minimal reliance on the ansatz for its functional form. We illustrate our method using the Via Lactea II N-body simulation as well as an analytical model for the dark matter halo. We conclude that O(1000) events are necessary to measure deviations from the Standard Halo Model and constrain or measure the presence of anisotropies.


I. INTRODUCTION
Dark matter comprises 80% of the mass of the Universe and its presence has been inferred gravitationally in numerous different physical settings including galactic rotation curves, precision studies of the cosmic microwave background radiation, and through gravitational lensing of galaxies and clusters of galaxies. However, dark matter's non-gravitational interactions remain unknown and there is an active search underway to determine dark matter's identity.
Many well motivated proposals to extend the Standard Model accommodate naturally a dark matter candidate that is stable, electrically neutral and has the correct relic abundance [1]. These models generically predict that dark matter should interact with Standard Model particles with electroweak-sized cross sections. Currently, a large experimental effort is underway searching for interactions of dark matter with nuclei in underground experiments.
Typically, direct detection experiments attempt to infer that a dark matter particle scattered off a nucleus by detecting the nuclear recoil and measuring the energy deposited in the detector. While most presently running experiments were designed to only measure the energy of the nuclear recoil, more can be learned about dark matter if the direction of nuclear recoils is measured as well [2]. The anisotropy of dark matter-induced nuclear recoils partially arises from the motion of the Earth through galactic dark matter halo, but can also arise from intrinsic anisotropies in the dark matter halo [3]. Therefore, measuring the direction of the nuclear recoils is not only a powerful handle to discriminate signal from background but also provides a whole class of new measurements of the dark matter's velocity distribution in the solar neighborhood [4][5][6][7][8][9].
The majority of the upcoming directional detection experiments use time projection chambers (TPCs) to resolve the direction of the nuclear recoils, by means of drifting the ionized track originated from the nuclear recoils through an electric field and projecting it over a charge collection grid. Two of the dimensions are then inferred from this projection, while the third dimension is obtained by the drift time of the particles in the ionized track.
In order to resolve the direction of the track to good accuracy, these TPCs have to operate at low pressure and density, and therefore are limited in target mass, as well as in volume, since the drift length is limited by diffusion [13,16]. The limitation in exposure of directional detection experiments makes them currently not competitive with current experiments such as XENON100 [24] and CDMS [25], amongst many others.
Most likely, directional detection experiments will not be able to scale their exposures sufficiently fast to be the first to discover dark matter. However, they will provide unique measurements of dark matter properties, once dark matter is discovered through direct detection. The unique capabilities of directional detection experiments have been underemphasized. For instance, directional detection experiments will explore the kinematics of the interaction between dark matter and Standard Model nuclei more directly than direct detection experiments [26][27][28][29][30][31].
This article will explore how directional detection experiments can measure astrophysical properties of dark matter, such as its phase space distribution in the solar neighborhood.
Previous studies in this domain have been carried out in [32][33][34][35][36]  The organization of this article is as follows. In Sec. II we overview general aspects of directional detection in the context of post-discovery. In Sec. III we discuss Jeans theorem and how it motivates integrals of motion as powerful variables in terms of which the dark matter distribution function can be expressed. In Sec. IV we present a novel method to parameterize the functional form of the distribution function -instead of relying on ansatzes of its functional form, we perform a special function decomposition of the dark matter phase space distribution. The coefficients of this decomposition will be the measurable quantities in a directional experiment. We perform fits of these expansion coefficients in Sec. V by simulating data from a mock experiment, and discuss the level of precision in the measurement the velocity distribution function that can be achieved with plausible numbers of events at upcoming directional detection experiments. In Sec. VI we test our method using the Via Lactea II N-body simulation as a template of the dark-matter distribution in the solar neighborhood, by scattering the particles of the simulation at a detector placed at 8kpc from the center. We then use this mock signal to recover the input distribution, providing a closure/bias test for the procedure. Finally, in Sec. VII we discuss the implications for upcoming experiments.

II. GENERAL ASPECTS OF DIRECTIONALITY
In this section we discuss how the direction of nuclear recoils captures the underlying dynamics of the interaction between dark matter and nuclei. We will also show how anisotropies in the dark matter velocity distribution are manifest in rate anisotropies.

A. Directional dependence from form factors and kinematics
This subsection reviews the basic kinematics of dark matter direct detection. Consider a dark matter particle with mass m dm and initial velocity v dm scattering off of a nucleus at rest with mass m N . The angle between the nuclear recoil direction and v dm , η, is given by: Above, v min (E R ) is the minimum dark matter velocity required by energy and momentum conservation in order for the nucleus to be scattered with recoil energy E R . If the kinematics of the scattering is elastic, v min (E R ) is given by where µ is the reduced mass of the dark matter-nucleus system. Fig. 1  Notice that given v dm , the angle η and the nuclear recoil energy E R are not uniquely determined. In fact, given v dm , the recoil energy can range in the interval for elastic scattering. Given a flux of dark matter particles over the detector, the nuclear recoil rate is determined by Eq. 1 convolved with the dark matter velocity distribution and the scattering probability density over the energy recoil range in Eq. 3. If this probability distribution function (pdf) is flat in E R , the corresponding probability distribution over η will range from 0 to π/2, peaking at π/4, as illustrated in Fig. 2.
However, generically the scattering cross section will have a dependence on energy, modifying the probability that, given v dm , the nucleus will recoil with energy E R . For instance, the finite size of the nucleus will introduce a nuclear form factor dependence that suppresses the scattering probability at high recoil energies [37]. There may also be a form factor associated with the nuclear spin dependence of the interaction [38][39][40], or on the dark matter properties, for instance if the dark matter has a dipole moment or is a composite state [41][42][43][44][45][46][47][48]. Such form factors will alter not only the nuclear recoil energy distribution, but the nuclear recoil direction as well. v dm . A flat scattering probability over the kinematically allowed energy range will translate into a nuclear recoil directional pdf peaking at an angle of π/4 with respect to the incident dark matter direction. Form factors will introduce an energy dependence of the scattering probability, distorting not only the recoil spectrum but also the nuclear directional rate. Kinematic properties such as inelasticity will not distort the pdf's, but truncate them over the kinematically allowed range. A dark matter form factor depending on a positive power of the momentum transfer, such as q 2 , will have the opposite effect, suppressing the scattering probability at low energies and shifting the recoil direction rate towards small angles with respect to the dark matter incident direction.
Another interesting possibility is that the scattering is inelastic, and the dark matter particle up-scatters to an excited state with mass splitting δ relative to the ground state [49]. In this case v min (E R ) will be given by The presence of inelasticity will not alter the pdf's shape in E R and η for a given As mentioned in the previous subsection, the energy dependence of the dark matternucleus interaction affects the directional rate of nuclear recoils. In this section we discuss the additional effect of anisotropies in the velocity distribution.
So far we considered the effect of dark matter incidence with fixed velocity and direction.
That scenario may be realized if the DM distribution in the solar neighborhood is dominated by a stream with low velocity dispersion [50][51][52]. Another possibility is that the dominant dark matter component in the solar neighborhood is isotropic and at equilibrium. As a first approximation we can consider a thermal velocity distribution truncated at some escape velocity, v esc . If the dark matter halo and the solar system were co-rotating, this distribution would yield an isotropic nuclear recoil rate. However, if the solar system is moving through a static halo, the flux of dark matter particles will be boosted in the direction of the Earth's motion, creating interesting effects such as an annual modulation in the rate and a large anisotropy in the nuclear recoil direction, which will peak in the opposite direction to the Cygnus constellation. Anisotropic velocity distributions may generate even more interesting features in the directional rate of nuclear recoils, such as rings [53] and hot spots.
In order to illustrate those features, we begin by establishing formulas and conventions that will be used throughout this paper. Given a dark matter velocity distribution f ( v) in the galactic frame, the nuclear recoil rate as a function of energy and direction is given by [54][55][56][57]: where v ⊕ is the Earth's velocity in the galactic frame. N contains the dependence on the scattering cross-section at zero momentum-transfer, on the local dark matter density and on phase-space factors. F (E R ) parameterizes the cross-section at finite momentum transfer as a product of nuclear and dark matter form factors. Finally,v R is the direction of the nuclear recoil. We will ignore the Earth's orbit around the Sun and choose coordinates such v r ≡ ( sin θ cos φ , cos θ , sin θ sin φ ).
In those coordinates the z-direction is aligned with the axis of the disk and the x-direction points away from the center of the galaxy as depicted in Fig. 3. We will also neglect the eccentricity of the Sun's orbit around the center of the galaxy. In that limit the Earth's velocity is parallel to the y-direction.
Note from Eq. 5 that, if the velocity distribution f ( v) is isotropic in the galactic frame, then the only source of anisotropy in the rate will come from the Earth's boost, v ⊕ ·v R = v ⊕ cos θ. Therefore, in the coordinate system we adopted, isotropic velocity distributions should induce rates that have a dependence on θ but not on φ. That means that any φ dependence on the rate is a smoking gun signature of an anisotropic velocity distribution.  For illustration purposes, we compare the angular dependence of the nuclear recoil rate for an isotropic versus an anisotropic distribution. As pointed out, the predicted rate for the isotropic distribution is flat on φ, and has a dependence on θ that varies with the recoil energy. If the energy is low, the recoil direction will prefer high values of θ, given Eq. 1 and the fact that v min (E R ) decreases as E R decreases. On the other hand, increasing the recoil energy will shift the rate towards lower values of θ. That is illustrated pictorially in Fig. 4 low recoil energies exhibit "cone-like" features in the recoil direction that point away from the incoming dark matter flux. Increasing the recoil energy has the effect of collimating the "recoil cone", until the ring-like structure finally shrinks to a hot spot.
Now consider an anisotropic velocity distribution for which the tangential components of the velocities are suppressed relative to the radial component. In that case circular orbits are subdominant, and the preferred flux of dark matter particles is along the radial direction in the galactic frame. In the Earth's frame, however, the flux picks up a component along the Earth's motion. That generates two overlapping "recoil cones" and hot spots in the recoil rate. As in the isotropic case, increasing the recoil energy will shrink the recoil cones to two more pronounced hot spots that will be shifted in φ by π/2 relative to the low energy hot spots. That effect is illustrated in Fig.5.

III. INTEGRALS OF MOTION AND JEANS THEOREM
In the previous section we discussed how the particle and astrophysical properties of dark matter manifest themselves in directional nuclear recoil signatures. Since very little is known about those properties, educated guesses have to made in order to explore direct detection signals. However, in the event that a dark matter signal is discovered in direct detection, we would like to take the data at face value and interpret it free of theoretical prejudice whenever possible.
In this section we will develop the tools for the main goal of this work: to infer the dark matter phase space distribution in the solar neighborhood from directional data. The first step towards this goal is to choose the correct variables on which to express the phase space distribution function. Inquiring how the dark matter halo was formed in the first place will give us some guidance. For instance, since the Milky Way halo has not undergone recent mergers, we can assume that the majority of dark matter particles has had enough time to relax to a state of equilibrium, in which case Jeans Theorem should hold to a good approximation [58]: The Jeans theorem is the inspiration for the parameterization of the dark matter DF that will be used in this work. If a measurement of the local dark matter distribution function were our sole concern, any parameterization in terms of functions of the phase space coordinates could be adopted, even if such functions were not the true integrals of motion. But in scenarios where the equilibrium conditions of Jeans Theorem are valid, expressing the DF in terms of the true integrals of motion provides information beyond the scales of the solar neighborhood. In principle, we could analytically extend a measurement of the phase distribution in the solar system to any other position in the halo, and therefore infer the shape of the global DF through . . I n being integrals of motion of the halo.
Moreover, the strong version of Jeans Theorem allows us to assume that, for all practical purposes, the phase-space distribution will be a function of at most 3 integrals of motion.
Quasi-static haloes, for which Jeans theorem applies, always have energy, E, as an integral of motion. Other possibilities include the components of the angular momentum, L, and possibly non-classical integrals which have no analytical expression. Interestingly, the directional recoil rate is easily integrated if the distribution function depends only on energy.
Consider a generic DF f (E). According to Eq.5, the directional rate is given by where we have defined and the energy E is given in terms of the dark matter's velocity and the local galactic escape velocity v esc : The factor Θ(−E) excludes the possibility that the there are unbound particles (for which On Appendix A we show that, Hence, expanding f (E) in a Taylor series, using Eq. 12 and re-summing the expansion we finally obtain: Eq. 13 is a remarkable relation between the directional nuclear recoil rate and the underlying dark matter distribution. It makes it clear that, if the DF is isotropic in the solar neighborhood, measuring the dependence of the rate on v 2 esc /2 − ρ(E R , θ, φ) 2 /2 is equivalent to directly measuring the underlying DF, since the later is the derivative of the rate with More generically, however, the dark matter DF will be anisotropic and depend on other integrals of motion, for instance the magnitude of the angular momentum | L| or its components, such as L z . Unfortunately, no straightforward expression for the rate has been derived for a generic function f (E, | L|, L z ). However, we have been able to analytically integrate the rate for a certain class of functions (see Appendix A). Such functions have polynomial dependence in E and harmonic dependence in the components of the angular momentum. That will prove useful for the method we will propose in the next section, where we will perform a special function decomposition of the DF.

IV. SPECIAL FUNCTION DECOMPOSITION OF THE HALO DISTRIBUTION
Although the reconstruction of the dark matter distribution function from directional and non-directional detection signals has already been investigated for specific halo types [32][33][34][35][36], no model independent reconstruction study has been performed yet (see, however, halo independent proposals to interpret direct detection signals [59][60][61]). One of the main goals of this paper is to provide a parameterization of the DF that is as generic as possible.
The method we propose is to decompose the distribution function in an orthonormal basis of special functions, and fit the measured rate to the coefficients of this decomposition. If the basis is appropriately chosen, the higher order coefficients will be subdominant and we will be able to truncate the series to a finite and preferably manageable number of coefficients.
There is no optimal choice for the variables and for the basis of decomposition that applies for all possible cases, and our main purpose is to illustrate out method rather than explore all possible variables and bases. Therefore, in this paper we will choose the following integrals of motion on which to parameterize the local dark matter distribution function: • The energy E, which, as before, is given by • The z-component of the angular momentum, L z , which in the coordinate system we have adopted is given by In the solar neighborhood z = 0 and hence • The magnitude of the angular momentum, whose DF dependence will be incorporated Again, in the solar neighborhood that reduces to L t = r ⊕ v θ .
We will make a further assumption which is not generic but that will greatly facilitate the computation of the rate and reduce the number of decomposition coefficients in the fit.
This assumption is that the DF f (E, L t , L z ) is separable: Under this assumption, we will adopt the following choice of basis for f 1 (E), f 2 (L t ) and Above, f 1 (E) is expanded in shifted Legendre polinomials,P , and E lim is the lowest dark matter energy that an experiment with recoil threshold E th R is able to probe: In principle we could have chosen E lim to the lowest possible energy, −v 2 esc /2. However, since a realistic experiment will have no information about particles with energy less than E lim , there is no purpose in parameterizing the distribution function beyond what is measurable.
The functions f 2 (L t ) and f 3 (L y ) are decomposed in a Fourier series, where L max ≡ r ⊕ v esc .
A great practical advantage of the parameterization above is that the recoil rate can be analytically integrated. That way no numerical integration is necessary, reducing computational time to perform the fits. In Appendix A we derive and list the analytical forms for the rate on those basis functions.
In the following section we will illustrate how the fit is performed for a hypothetical experiment measuring the recoil rate for an anisotropic DF. We will choose as an example a Michie distribution [62] and generate mock up Monte Carlo data as the detector would measure it. We will then employ our method to recover the original distribution. In Sec. VI we will use the Via Lactea II simulation as the underlying dark matter distribution and simulate the rate at the detector. We will then employ our method to perform the fits.

V. DIRECTIONAL FITS TO A MICHIE DISTRIBUTION
In this section we illustrate our method for a hypothetical dark matter signal in a directional detector. Ideally we would like to have a large number of events to perform the fit. However, the current limits are quite stringent for spin-independent elastic scattering dark matter particles with mass of O(100) GeV. For instance, the latest XENON100 run with 4 kg-yrs of exposure observed 3 candidate events with an expected background of 1.8 ± 0.6 events [24]. That suggests that, in the heavy WIMP scenario where discovery is just around the corner, the detection of O(100) events would require at least hundreds of kg-yrs of exposure, which is still an unrealistic goal for the O(1)-kg detectors planned by current directional experiments. The most optimistic possibility is if dark matter is light, of order of a few GeV. In this mass range the limits are much less stringent [63,64], and the relatively light targets considered by directional experiments, such as CF 4 and CS 2 , are ideal.
Therefore we choose to generate this hypothetical signal for an elastically scattering dark matter particle with mass m dm = 6 GeV. We assume that the scattering is off of a CS 2 target, with recoil energy acceptance window E R ∈ [5 keV, 15 keV] and for simplicity we ignore the scattering off of carbon. We assume that the experiment operates in a zero background environment, but in principle background events could be straightforwardly incorporated in the fits if their directional shape is known. Finally, the results we will show in the following are obtained from fits without resolution effects taken into account. As a crude approximation of the detector angular resolution, we have divided the nuclear recoil direction in bins of size (10 • , 10 • ) and did not observe any significant alteration of our results.
In the results shown in the rest of this article, no binning has been used.
In order to illustrate the interesting effects of anisotropy, we choose a Michie distribution for the dark matter. In terms of integrals of motion, the Michie distribution is given by: Our choice of parameters is the following: v 0 = 280 km/s, v esc = 600 km/s, L 0 = r ⊕ v 0 and α = 1. We will assume in the fits that local escape velocity v esc as well as the dark matter mass m dm are know to sufficient precision. That is a reasonable assumption since directional experiments will only come into play after the recoil spectrum has been well measured by several non-directional experiments with different nuclear targets.
Note that the Michie distribution (18) trivially satisfies the separability condition (15) when expressed in terms of E, L t and L z , since | L| 2 = L 2 t + L 2 z . Then, with we can perform the special function decomposition as in Eq.16 and obtain the coefficients: The overall normalization of the coefficients above is unphysical since it can be reabsorbed in the cross section. Hence, the coefficients were scaled such that f 1 (E), f 2 (L t ) and f 3 (L z ) are normalized to unity. We have truncated the energy expansion in Legendre polynomials to the 5th order, and the Fourier decomposition of the components of the angular momentum to the 3rd harmonic. The truncation of the expansions at this orders approximates the true functions to an accuracy better than the percent level.
Knowing the true values of the decomposition coefficients, we can compare them to the values obtained from fits of the rate. We generated several Monte Carlo data ensembles, varying the number of observed events. We then performed unbinned log-likelihood fits in the 11-dimensional parameter space of decomposition coefficients using the publicly available code MultiNest v2.10, a Bayesian inference tool that uses a nested sampling algorithm to scan the parameter space [65,66].
The 1-σ range of the fitted coefficients can be extracted by approximating the loglikelihood function near its maximum by a quadratic function of the parameters of the fit: where C is the vector representing a given set of coefficients and σ is their standard deviation. In Fig. 7 we show the results of the fit for different numbers of events, varying from 100 to 10,000. We display the 68% C.L. contour intervals of the fitted coefficients, as well as the true values for comparison. As the errors are correlated, the 68% range of the decomposition coefficients is actually a multidimensional ellipsoid. The different plots in Fig. 7 are 2D projections of this ellipsoid in different planes of the parameter space.
As one can see, a few thousand events are necessary for a good measurement of the underlying dark matter distribution. Fig. 6 displays the functional forms of f 1 (E), f 2 (L t ) and f 3 (L z ) resulting from fits with 1000 events and 10,000 events. This shows that our method works remarkably well if the variables and decomposition basis are chosen judiciously. In the real world, however, an appropriate choice of decomposition basis may not be as trivial as in this hypothetical example. However, several checks of the data can be made before performing fits to test whether a given hypothesis is a good approximation. We will illustrate that in the next section, where we will use the Via Lactea II N-body simulation as a template for the dark matter distribution in our galaxy, and "scatter" the particles of this simulation is a detector placed at a location around 8 kpc from the center to try to infer the underlying phase space distribution.

VI. FITS TO THE VIA LACTEA II SIMULATION
No dark matter direct detection signal has been observed to this date to provide us with clues on the local distribution. Most of our current understanding of the range of possibilities for the local distribution comes from N-body simulations [3,50,51,67], which also provide one of the best test grounds to investigate the formation and dynamics of dark matter halos [68]. With about 10 9 particles with masses around 10 3 solar masses, Via Lactea II is currently one of the highest resolution N-body simulations and can be used as a model of the Milk Way halo [69].
In this section we will test the method proposed in this study using Via Lactea II as the underlying dark matter distribution. Ideally, we would like to choose a box located at R ⊕ = 8 kpc away from the center of the halo (which is the distance of the solar system from the center of the Milky Way) and with a volume (∆r) 3 R 3 ⊕ where ∆r is the length of a side of the box. Unfortunately, in Via Lactea II a region of this size typically does not enclose any particles. Our approach is then to choose a spherical shell of 8kpc-radius and 200pc-thickness, which encloses enough particles to generate a dark matter signal at a directional detector. Since we do not want to wash out anisotropies, we choose a fixed position r ⊕ at which to place the detector, and for each dark matter particle in the shell with coordinates v and r, we perform an affine map v( r) → v( r ⊕ ).
In the same spirit of the previous section, we would like to fit a signal generated from Via Lactea II to a distribution function in terms of E, L t , L z , which is separable as in Eq. 15, and decomposed as a Legendre series in E and a Fourier series in L t , L z , see Eq. 16. However, before attempting to perform those fits we should test the validity of those assumptions.
The first of those assumptions, that the energy and angular momentum are integrals of motion, can be shown to be a good approximation for positions not too far away from the center of the halo, where the gravitational potential can be reasonably approximated as being spherical. A quantitative measure of the sphericity of the potential is given by Above, ψ( r) is the exact gravitational potential at the position r, and ψ s (r) is the potential computed assuming spherical symmetry and evaluated using Gauss' law. χ( r) then parametrizes the deviation of the gravitational potential at position r from a spherically symmetric potential. and r = rẑ. One can see that for r < 200 kpc, the assumption of a spherically symmetric potential leads to an error O(6%). For r < 30 kpc, the error is reduced to O(3%). Hence the energy and angular momentum of dark matter particles within a 30 kpc radius should be approximately conserved.
The second of our assumptions, that of separability, can be evaluated through the function where which measures the correlation between the three variables E, L t and L z . If the dark matter distribution function verifies Eq. 15, G(E, L t , L z ) should be constant and equal to 1. Due to the lack of statistics, in Fig. 9 we plotḠ E ,Ḡ Ly ,Ḡ Lz , the averages of G over the two other variables, rather than the tridimensional behavior of the function. The fact thatḠ E ,Ḡ Lt ,Ḡ Lz are consistent with a constant value equal to one reassures us that our assumption of separability is well justified.

A. Local Fits to Via Lactea II
The range of validity of our hypotheses described in Sec. III were quantified in the previous section. In this section we perform fits to a directional detection signal generated from the distribution of VLII simulation, using particles contained within a shell located 8 kpc away from the center of the simulation, and with thickness of 200 pc. As mentioned previously, we fix the location for the detector at r ⊕ , and transform the velocities of all particles in the shell through an affine map v( r) → v( r ⊕ ).
As in Sec V, we set the mass of the dark matter particles to m dm = 6 GeV. Inspection of the VLII distributions allows us determine that the escape velocity at 8 kpc is v esc ≈ 545 km/s. We generate 10,000 events by elastically scattering the particles in the shell off of a CS 2 target, but ignore the scattering off of carbon. However, differently from Sec V, we are forced to lower the energy threshold of 5 keV. The reason for that is that with a 5 keV threshold, only particles in the tail of the energy distribution (v > ∼ 518 km/s) are kinematically allowed to scatter. Unfortunately there only about 10 5 particles inside the shell, and the tails of the distributions of those particles are not well populated. Hence, in order to avoid large errors resulting from fitting the tails, we lower the energy threshold to 2 keV, allowing particles with v > ∼ 330 km/s to scatter. In that range the distribution is better populated and the fits are more reliable. As in Sec V, detector resolution effects are not taken into account.
FIG. 10: θ − φ profile of the directional signal generated from VLII simulation. Notice that the hint of hot spots at φ = 0 [π] suggests the presence of anisotropies.
The θ − φ profile of the directional detection rate generated is shown in Fig. 10. The hint of hot spots at φ = 0 [π] suggest a suppression of circular orbits perpendicular to the plane of the Earth's orbit in the halo. At this point this is just a qualitative statement, which will be made quantitative below by the results of the fits. Note that the plane of the Earth's orbit was chosen arbitrarily given that VLII does not have a baryonic disk.
We perform unbinned log-likelihood fits to this directional signal as described in Secs. IV,V. The results of the fits are displayed in Fig. 11, where we plot E, L t and L z distributions for the particles that are kinematically allowed to generate nuclear recoils within the recoil energy acceptance window. We contrast the distributions extracted from the fit with the true underlying distributions to show that the fit is correctly capturing the local VLII phase space distribution within errors. The fits for the E, L t and L z dependent The most straightforward variable to check is E. Consider the distribution of dark matter particles at a given location r, binned in E, L t and L z with bin widths ∆E, ∆L t and ∆L z respectively. The number of particles in the bin centered around E, L t and L z is given by: where Det[J(E, L y , L z )] is the determinant of the Jacobian, given by: Since the bins with better statistics are those with lower values of E, and hence lower values of L t and L z , we can extract f 1 (E) from bins where L t and L z are centered around zero and E is centered around E i : The first factors in Eq. 34 are independent of E and can be normalized away, finally giving: In Fig. 12, we plot f 1 (E) using Eq. 35 for two different positions in the VLII simulation, Unfortunately, there is not enough statistics to extract f 2 (L t ) and f 3 (L z ) directly from the simulation using the same method used above for f 1 (E), since the population in bins with increasing values of L t and L z rapidly drop in number. The alternative is to fit the dark matter distributions at different radii using our parameterization in Eqs. 15 and 16.
We perform such fits for three different radii, r = 4.5 kpc, r = 8 kpc and r = 30 kpc, where as before we select particles within shells of 200 pc-thickness and radius r. Fig. 12 shows the results of the fits for f 1 (E), f 2 (L t ) and f 3 (L z ) and their corresponding error bands. The fact that there is a reasonable agreement over a wide range range of galactic radii is evidence that the parameterization motivated by Jeans theorem is valid and meaningful. From Fig. 12 it is also interesting to note that, although f 3 (L z ) is consistent with a flat distribution, f 2 (L t ) indicates that the VLII distribution is anisotropic, and the angular momentum component along the Earth's motion is suppressed. Another interesting fact is that f 1 (E) is not consistent with a single exponential, as is the case for Maxwellian or Michie distributions. Indeed, the high velocity tail of f 1 (E) is harder than Maxwellian, as pointed out in [67].

VII. DISCUSSION
In this paper we discussed important aspects of directional detection that will be relevant once dark matter is detected in underground experiments. Directional experiments will probe the dynamics and kinematics in the dark matter interactions with nuclei and allow us for directly measure the dark matter phase space distribution in the solar neighborhood.
It will be important to reduce our bias in extracting information from directional data.
With that in mind, we proposed to parametrize the local distribution function in terms of integrals of motion. If the dominant component of underlying dark matter distribution in the solar neighborhood is in equilibrium, expressing the DF in terms of integrals of motion will allow us to infer the DF at other positions in the halo. That can have a potentially large impact in studies of indirect signals, as well as implications for our understanding of the formation of the Milky Way halo. There is a chance that the local dark matter distribution could have a sizable component that is in a non-equilibrium configuration such as streams of dark matter arising from tidally disrupted dark matter sub-halos. These features can be readily identified in directional detection experiments and can be subtracted off [50][51][52].
We also proposed a decomposition of the distribution function in special functions of integrals of motion. That allows for a more model independent extraction of the true phase space distribution. We illustrated our method by generating Monte Carlo signals in a hypothetical directional detector and fitting the coefficients of our decomposition. This study was performed with an analytical model of the halo, as well as using Via-Lactea II simulation. In the later, we also tested the results of our fit by extrapolating the fitted local distribution to other positions in the halo and comparing to the true distribution. We found that for radii less than ∼ 30 kpc our extrapolation worked to a very good approximation.
The number of signal events necessary to detect non-standard features of the dark matter distribution obviously depend on the degree of such effects. For Via Lactea II we found that O(1000) events were necessary to detect presence anisotropies and departures from a Maxwellian distribution.
The parameterization and decomposition of the dark matter phase space distribution we have proposed may also have applications beyond directional detection. For instance, it could also allow a better understanding of the dark matter haloes obtained in N-body simulations.
In particular, reconstructing the dark matter phase space distribution would allow us to compare the different existing N-body simulations to each other, studying for example how the anisotropies of the dark matter distribution function depend on the different parameters of the simulation, especially the presence of a baryonic disk.
In conclusion, this article has shown the power of directional detection experiments in the post-dark matter-discovery era. By using the results from directional detection experiments and N-body simulations, a better understanding of dark matter in the Milky Way can be obtained.
Usingv R = (x R , y R , z R ) ≡ (sin θ cos φ, cos θ, sin θ sin φ), and integrating over v x gives The domain of integration for the variables v y and v z is then set by the condition The limits of integration on v z are then v z± (v y ) = 1 where The condition ∆ > 0 sets the limits for v y : v y± = y R ρ(E R , θ) ± (1 − y 2 R )(v 2 esc − ρ 2 (E R , θ)).
With the limits of integration defined, the only remaining constraint is v 2 esc − ρ 2 (E R , θ)) > 0, which requires that for a given recoil energy E R and recoil angle θ, the minimum velocity of the incoming dark matter particles in the galactic frame must be smaller than the dark matter escape velocity.
We then do the variable change The new boundaries are then given by ±V 0 z (V y ) and ±V 0 y with where The differential rate then becomes For a distribution function expressed in terms of energy E and the components of the angular momentum L t and L z , we can recover the dependence on V y and V z using Eqs. A15 to A18 allow us to integrate the rate analytically for several special forms of f (E, L t , L z ), such as E cos(nπL z ) cos(mπL t ), E cos(nπL z ) sin(mπL t ), (A20) E sin(nπL z ) cos(mπL t ), E sin(nπL z ) sin(mπL t ).
For instance, for f (E, L t , L z ) as in Eq. A19, the differential rate is given by cos mπr ⊕ |x r | x 2 r + z 2 r (V z + v 0 z ) cos nπr ⊕ (V y + v 0 y ) dV z dV y .
Remembering thatv we can rewrite the expression above in a frame independent form: where u T is an arbitrary vector in the tangential plane.
Finally, when there is no dependence on L the rate takes the form: In particular, for a generic f (E),