<figure><img src="image_1.png" /></figure>

HAL
open science

# Discretising the velocity distribution for directional dark matter experiments

Bradley J. Kavanagh

► To cite this version:

Bradley J. Kavanagh. Discretising the velocity distribution for directional dark matter experiments. Journal of Cosmology and Astroparticle Physics, 2015, 2015 (07), pp.019. 10.1088/1475-7516/2015/07/019. cea-01335007

HAL Id: cea-01335007

https://cea.hal.science/cea-01335007

Submitted on 5 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Discretising the velocity distribution for directional dark matter experiments

Bradley J. Kavanagh

Institut de Physique Théorique, CNRS, URA 2306 & CEA/Saclay,
F-91191 Gif-sur-Yvette, France

School of Physics & Astronomy, University of Nottingham,
University Park, Nottingham, NG7 2RD, U.K.

E-mail: bradley.kavanagh@cea.fr

Received March 2, 2015

Revised May 13, 2015

Accepted June 16, 2015

Published July 13, 2015

**Abstract.** Dark matter (DM) direct detection experiments which are directionally-sensitive may be the only method of probing the full velocity distribution function (VDF) of the Galactic DM halo. We present an angular basis for the DM VDF which can be used to parametrise the distribution in order to mitigate astrophysical uncertainties in future directional experiments and extract information about the DM halo. This basis consists of discretising the VDF in a series of angular bins, with the VDF being only a function of the DM speed $v$ within each bin. In contrast to other methods, such as spherical harmonic expansions, the use of this basis allows us to guarantee that the resulting VDF is everywhere positive and therefore physical. We present a recipe for calculating the event rates corresponding to the discrete VDF for an arbitrary number of angular bins $N$ and investigate the discretisation error which is introduced in this way. For smooth, Standard Halo Model-like distribution functions, only $N = 3$ angular bins are required to achieve an accuracy of around 10–30% in the number of events in each bin. Shortly after confirmation of the DM origin of the signal with around 50 events, this accuracy should be sufficient to allow the discretised velocity distribution to be employed reliably. For more extreme VDFs (such as streams), the discretisation error is typically much larger, but can be improved with increasing $N$. This method paves the way towards an astrophysics-independent analysis framework for the directional detection of dark matter.

**Keywords:** dark matter experiments, dark matter detectors, dark matter theory

**ArXiv ePrint:** 1502.04224

<figure><img src="image_3.png" /></figure>

Article funded by SCOAP³. Content from this work may be used
under the terms of the Creative Commons Attribution 3.0 License.

Any further distribution of this work must maintain attribution to the author(s)
and the title of the work, journal citation and DOI.

doi:10.1088/1475-7516/2015/07/019
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Directional formalism</td><td>3</td></tr><tr><td>3 Problems with parametrising the velocity distribution</td><td>5</td></tr><tr><td>4 A discretised velocity distribution</td><td>6</td></tr><tr><td>5 Comparison with exact results</td><td>10</td></tr><tr><td>5.1 N=2 discretisation</td><td>11</td></tr><tr><td>5.2 N=3 discretisation</td><td>13</td></tr><tr><td>5.2.1 The folded rate</td><td>15</td></tr><tr><td>5.3 N=5 discretisation</td><td>16</td></tr><tr><td>5.4 Forward-backward asymmetry</td><td>17</td></tr><tr><td>6 Discussion</td><td>19</td></tr><tr><td>7 Conclusions</td><td>21</td></tr><tr><td>A Calculating the Radon transform</td><td>22</td></tr><tr><td>B The directionally-integrated Radon transform (IRT)</td><td>28</td></tr></table>

## 1 Introduction

The measurement of the directionality of nuclear recoils in dedicated low-background detectors has long been proposed as a method to detect particle dark matter (DM) [1]. The motion of the Sun through the Galactic DM halo generates a so-called 'wind' of DM particles, which would appear to originate from the direction of the constellation of Cygnus. The resulting recoil spectrum will be peaked in the opposing direction (or for high mass DM, in a ring around this direction [2]). Such a directional signal would be a smoking gun for DM against the typically isotropic distribution of residual backgrounds [3-7]. Directional detection may also provide a way to circumvent the neutrino floor [8] and to probe the local structure of the DM velocity distribution [9, 10].

A number of directional experiments are currently in the prototype stage. The predominant technology is the low-pressure time-projection chamber (TPC), in which nuclear recoils leave an O(1 mm) track of ionised electrons, which can be imaged by drifting the electrons (or electron-transport gas) to an anode grid. This technology is currently utilised by the DRIFT [11, 12], MIMAC [13, 14], DMTPC [15, 16], NEWAGE [17, 18] and D3 [19] collaborations and the problem of reconstructing the recoil direction using this technology has been much studied [20, 21]. More recent proposals for directional detectors include nuclear emulsions [22], electronic scattering in crystals [23], DNA-based methods [24] and columnar recombination in Xenon targets [25, 26].

Directional experiments such as these may be the only possibility for probing the full 3-dimensional velocity distribution of DM, $f(v)$, which is a priori unknown. The standard
analysis of direct detection experiments assumes the so-called Standard Halo Model (SHM), a simplified model for the DM halo as an isotropic, isothermal sphere of particles [27]. However, such a simple description is unlikely to be an accurate description of the Milky Way halo [28–32]. In addition, N-body simulations raise the possibility of more complicated distributions, including debris flows [33, 34], tidal streams [35, 36] or a dark disk [37–39]. The impact of these astrophysical uncertainties on directional signals has previously been studied [2, 9, 10]. While non-standard astrophysics should not severely affect the ‘smoking gun’ directional signature of DM recoils [9], it may still pose a potential problem for future data. In particular, trying to reconstruct the DM parameters (such as mass $m_\chi$ and interaction cross section $\sigma$) from a potential signal requires some assumptions to be made about the form of $f(v)$. As with non-directional detection techniques, poor assumptions about the astrophysics of DM can lead to biased reconstructions of these parameters [40–44].

Astrophysics-independent approaches to non-directional experiments have received significant attention in the past few years, including so-called ‘halo-independent’ methods [45–55] and attempts at parametrising the 1-dimensional speed distribution [40, 42, 43, 56]. There have also been several attempts to construct a suitable parametrisation for the 3-dimensional velocity distribution, relevant for directional experiments [57–60]. This is significantly more difficult than in the non-directional case, owing to the presence of the angular components of $f(v)$. A single 1-dimensional function is therefore no longer sufficient to describe the full distribution function and a suitable angular basis must be found. Several bases have been suggested (e.g. refs. [59, 60]) but all typically suffer from the same problem, which has not previously been addressed in the literature: they allow $f(v)$ to take negative — and therefore unphysical — values, which can lead to spurious results. We discuss this problem in more detail in section 3.

In this work, we present a framework which allows the velocity distribution to be parametrised in a model-independent way, while guaranteeing that it is everywhere positive. This consists of an angular discretisation of the velocity distribution, decomposing it into a series of 1-dimensional functions, each of which is constant over a given bin in the angular coordinates. This is motivated in the first instance by its ability to describe the simplest directional signal, a forward-backward asymmetry in the scattering rate. However, we extend the discretisation to an arbitrary number of bins $N$ in the angular variables, such that in the limit of large $N$, the full 3-dimensional distribution can be recovered. For arbitrary $N$, we also demonstrate how to calculate the Radon transform (the function which encodes the angular dependence of the scattering rate) and therefore how to compute the full scattering rate.

Here, we focus on the angular discretisation itself, and do not attempt a full reconstruction of the DM velocity distribution and particle physics parameters (as in e.g. [40, 42]). Instead, we investigate the magnitude of the error introduced by this angular discretisation for two different benchmark velocity distributions. To do this, we compare the event rate within each angular bin obtained using the full and discrete velocity distributions. Qualitatively, we discuss and compare the energy dependence of the event rate within each bin. Quantitatively, we can calculate the number of signal events in each angular bin for the full and discrete distributions, and compare the result. This allows us to evaluate whether (and in which scenarios) fitting to data using such an angular discretisation gives accurate results.

In order to perform this analysis, we assume that the velocity distribution is known exactly. In the analysis of real experiments, a parametrisation of $f(v)$ within each bin must be chosen and fit to the data. We leave the choice of a suitable parametrisation of the remaining
1-dimensional functions of v to future work, meaning that the results presented here represent a lower limit on the error induced using the discretised velocity distribution. However, we lay out a framework which can be used in the analysis of future data to extract coarse-grained information about the DM velocity distribution in a way which was not previously possible.

In section 2, we present the directional detection formalism, including calculation of the scattering cross section and the Radon transform, and a brief discussion of the associated astrophysical uncertainties. This is followed in section 3 by a discussion of previous attempts to overcome these uncertainties. In section 4, we present the discretised velocity distribution which is the focus of this work. In section 5, we investigate how this discretised approximation compares to the exact result, both at the level of the Radon transform and the event rate. Finally, we discuss the significance of these results and plans for future work in section 6. The procedure for calculating the Radon transform from the discretised distribution is presented in appendix B, while we relegate the full derivation of this formula to appendix A.

## 2 Directional formalism

The directional rate per unit detector mass for recoils of energy $E_R$ can be written as [61]

$$ \frac{dR}{dE_R d\Omega_q} = \frac{\rho_0}{4\pi\mu_{\chi p}^2 m_{\chi}} \rho_T C_N F^2(E_R) \hat{f}(v_{\min}, \hat{q}), \quad (2.1) $$

where $\hat{q} = (\sin \vartheta \cos \varphi, \sin \vartheta \sin \varphi, \cos \varphi)$ is a unit vector pointing in the direction of the recoil. We have written the rate $R$ in terms of the DM-proton cross section at zero-momentum transfer $\sigma_p$, which may be spin-independent or spin-dependent; the form factor $F^2(E_R)$, which describes the loss of coherence due to the finite size of the nucleus; and an enhancement factor $C_N$, which describes the enhancement of the rate for a nucleus $N$ relative to the DM-proton rate. The local DM density is written as $\rho_0$, the DM mass as $m_{\chi}$ and the reduced DM-proton mass as $\mu_{\chi p} = m_{\chi}m_{\rho}/(m_{\chi} + m_{\rho})$.

For the gaseous targets used in low-pressure TPCs, relatively light nuclei are typically used. For example, the DRIFT experiment [11, 12] uses CF$_4$ as the target gas, with spin-dependent (SD) scattering from $^{19}$F nuclei expected to be the dominant interaction with DM particles. For SD interactions, the enhancement factor can be written as [62, 63]

$$ C_N = \frac{4J+1}{3J} |\langle S_p \rangle + a_n/a_p \langle S_n \rangle|^2, \quad (2.2) $$

where $J$ is the total nuclear spin, $\langle S_{p,n} \rangle$ is the expectation value of the proton and neutron spin in the nuclear ground state with maximal magnetic quantum number, and $a_n/a_p$ is the ratio of the DM-neutron and DM-proton couplings. The ratio $a_n/a_p$ depends on the specific model of DM under consideration, although simplifying assumptions are often used, including pure-nucleon couplings ($a_p = 0$ or $a_n = 0$) [64] or values motivated by specific models (e.g. $a_n/a_p = \pm 1$) [65].

The SD form factor can be written in terms of the spin structure functions of the nucleus $S_{ij}(q)$, which depend on the momentum transfer $q = \sqrt{2m_N E_R}$, for a nucleus of mass $m_N$:

$$ F_{\text{SD}}^2(q) = S(q)/S(0), \quad (2.3) $$

with

$$ S(q) = a_0^2 S_{00}(q) + a_0 a_1 S_{01}(q) + a_1^2 S_{11}(q). \quad (2.4) $$
The isoscalar and isovector couplings are related to the proton and neutron couplings by $a_0 = a_p + a_n$ and $a_1 = a_p - a_n$ [63]. The spin structure functions (and the proton and neutron spin matrix elements) can be calculated using nuclear shell models (see e.g. refs. [66–71]). Here, we focus on Fluorine targets, assuming $a_p = a_n$, with spin structure functions and matrix elements taken from ref. [72]. The spin structure functions lead to a roughly exponential suppression of the recoil spectrum as a function of recoil energy.

The DM velocity distribution enters into the scattering rate through the Radon transform

$$ \hat{f}(v_{\text{min}}, \hat{q}) = \int_{\mathbb{R}^3} f(\mathbf{v}) \delta(\mathbf{v} \cdot \hat{q} - v_{\text{min}}) \, d^3\mathbf{v}, \qquad (2.5) $$

where $v_{\text{min}}$ is the minimum DM speed required to excite a nuclear recoil of energy $E_R$:

$$ v_{\text{min}} = v_{\text{min}}(E_R) = \frac{s \overline{m_N E_R}}{2 \mu_{\chi N}^2} \qquad (2.6) $$

Geometrically, the Radon transform corresponds to an integral over $f(v)$ on a plane perpendicular to $\hat{q}$, which has a perpendicular distance $v_{\text{min}}$ from the origin. Physically, this corresponds to integrating over all DM velocities which satisfy the kinematic constraints for exciting a nuclear recoil of energy $E_R$ in direction $\hat{q}$.

The typical assumption for the form of the DM velocity distribution is the Standard Halo Model (SHM). For a spherically symmetric, isothermal DM halo, with density profile $\rho(r) \propto r^{-2}$, the resulting velocity distribution has a Maxwell-Boltzmann form in the reference-frame of the Galaxy, [61]

$$ f(\mathbf{v}) = \frac{1}{(2\pi\sigma_v^2)^{3/2}} \exp \left[ -\frac{\mathbf{v}^2}{2\sigma_v^2} \right], \qquad (2.7) $$

with velocity dispersion $\sigma_v$. The corresponding Radon transform also takes the form of a Gaussian,

$$ \hat{f}(v_{\text{min}}, \hat{q}) = \frac{1}{(2\pi\sigma_v^2)^{1/2}} \exp \left[ -\frac{v_{\text{min}}^2}{2\sigma_v^2} \right]. \qquad (2.8) $$

Using the properties of the Radon transform we can obtain the corresponding result in the Earth frame by performing a Galilean transformation $\mathbf{v} \rightarrow \mathbf{v} + \mathbf{v}_e$, with $\mathbf{v}_e$ the velocity of the Earth with respect to the rest frame of the halo, in which case

$$ \hat{f}(v_{\text{min}}, \hat{q}) \rightarrow \hat{f}(v_{\text{min}} + \mathbf{v}_e \cdot \hat{q}, \hat{q}). \qquad (2.9) $$

The SHM velocity distribution in the laboratory frame is then

$$ \hat{f}(v_{\text{min}}, \hat{q}) = \frac{1}{(2\pi\sigma_v^2)^{1/2}} \exp \left[ -\frac{(v_{\text{min}} + \mathbf{v}_e \cdot \hat{q})^2}{2\sigma_v^2} \right]. \qquad (2.10) $$

Within the SHM, one typically assumes values of $v_e \approx 220 \text{ km s}^{-1}$ and $\sigma_v = v_e/\sqrt{2} \approx 156 \text{ km s}^{-1}$. However, there are still uncertainties on these values of at least 10% [73–75]. Introducing a cut off in the distribution at the Galactic escape speed $v_{\text{esc}}$ introduces further uncertainty into the model [76].

Despite its widespread use, the SHM is unlikely to be an accurate representation of the DM halo. Observations and N-body simulations indicate that the halo should deviate from a $1/r^2$ profile and may not be spherically symmetric. As a result alternative models have been proposed. Speed distributions associated with triaxial halos [28] or with more realistic density
profiles [29] have been suggested, as well as analytic parametrisations which should provide more realistic behaviour at low and high speeds [30]. Self-consistent distribution functions reconstructed from the potential of the Milky Way have also been obtained [31, 32].

It is also possible to extract the speed distribution from N-body simulations. Such distribution functions tend to peak at lower speeds than the SHM and have a more populated high speed tail [77–79]. There are also indications that DM substructure may be significant, causing ‘bumps’ in the speed distribution, or that DM which has not completely phase-mixed — so-called ‘debris flows’ — may have a contribution [34].

Another result obtained from simulations is the possibility of a dark disk. When baryons are included in simulations of galaxy formation, this can result in DM subhalos being preferentially dragged into the disk plane where they are tidally stripped [37, 38]. The result is a dark disk which corotates with the stellar disk, though with a smaller velocity dispersion, typically $\sigma_v^{DD} \sim 50 \text{ km s}^{-1}$. Such a dark disk could comprise a large fraction (up to 50%) of the total DM density, though more recent results [39] suggest a smaller dark disk contribution (around 10%), depending on the merger history of the Milky Way.

Finally, direct detection experiments probe the DM halo on sub-milliparsec scales and there is therefore the possibility that DM sub-structures could dominate the local distribution. Analyses of N-body simulations suggest that no individual subhalos should dominate the local density [80, 81]. However, local structures could still be significant, especially streams of DM from the tidal disruption of Milky Way satellites, such as Sagittarius [35, 36, 82].

Such a wide range of possibilities for the form of $f(v)$ leads to substantial uncertainties on the predicted event rate for directional detectors. This in turn can lead to potential bias in the reconstruction of other DM parameters, if these uncertainties are not properly accounted for.

### 3 Problems with parametrising the velocity distribution

With promising developments in directional detector technology, it is interesting to ask what information about the velocity distribution we could, in principle, extract from a directional signal. Early attempts to address this question involved extending the SHM to allow anisotropic velocity dispersions and fitting these new parameters to mock data [57]. The possibility of including additional structures, such as streams and dark disks, has also been considered [58]. However, such methods are likely to be accurate only if the underlying velocity distribution is well-described by the chosen model.

Alves, Hendri and Wacker [59] investigated the more general possibility of describing $f(v)$ in terms of a series of special functions of integrals of motion (energy and angular momentum). These can then be fit to data, with around 1000 events required to distinguish between the SHM and a Via Lactea II distribution function [83]. However, the special, separable form of the velocity distribution requires that the dark matter halo is in equilibrium. More troubling is that the resulting velocity distribution is not guaranteed to be everywhere positive and therefore not all combinations of parameters correspond to physical distribution functions.

A more general parametrisation for the velocity distribution was recently proposed by Lee [60]. In this approach, the velocity distribution is decomposed into products of Fourier-Bessel functions and spherical harmonics. This is completely general and does not require assumptions about the halo being in equilibrium. Lee also gives an analytic expression for the Radon transform of the Fourier-Bessel basis, making this approach computationally efficient. However, this approach can also produce negative-valued, non-physical distribution functions, as in the case of the method of Alves et al.
In fact, any decomposition in terms of spherical harmonics leads to this problem, because the spherical harmonic basis functions can have negative values. It is unclear how this issue will affect parameter reconstruction. It may lead to parameter estimates (e.g. for $m_\chi$ or $\sigma_p$) which are unphysical, as they require an unphysical distribution function to fit the data well. Without some criteria which determines which coefficients of the spherical harmonics lead to strictly positive distribution functions, it may not be possible to reject such parameter points. We may attempt to numerically test each parametrised distribution function for negative values but for a real function of three parameters $f(v) = f(v_x, v_y, v_z)$ this would require a very large number of evaluations, which may not be computationally feasible (and does not absolutely guarantee positive-definiteness). In addition, physical distributions may occupy only a small fraction of the total space of parameters making parameter sampling and reconstruction difficult.

Even if positive-definiteness could be ensured, it is not clear how to interpret the parameters reconstructed in this way. Spherical harmonic approximations of typical velocity distributions such as the SHM (obtained by integrating out the coefficients of the basis functions) tend to produce distributions which do contain negative values. This is especially true when only a small number of basis functions is used. However, the 'true' spherical harmonic coefficients obtained in this way cannot be obtained by fitting to data (because we would reject those which lead to negative values). This may again lead to a bias in the reconstructed DM parameters, because the 'true' values do not lie in the allowed parameter space. Such potential problems have not previously been noted in the literature.

In order to fit to data, then, it is necessary to decompose $f(v)$ into a series of angular components $A^i$:

$$f(v) = f(v, \cos \vartheta^j, \varphi^j) = f^1(v)A^1(\vartheta^j, \varphi^j) + f^2(v)A^2(\vartheta^j, \varphi^j) + f^3(v)A^3(\vartheta^j, \varphi^j) + \dots\dots\dots\dots\dots (3.1)$$

We then truncate the series at some order, leaving only a finite number of 1-dimensional functions $f^i(v)$ which are unknown. This reduces the problem of attempting to fit a function of the 3-dimensional variable $v$ to the problem of parametrising a series of 1-dimensional functions, which is much more tractable. Of course, we should be careful that this truncation preserves enough angular information to still provide a good approximation to $f(v)$. However, as more data becomes available, we can add more terms to the series to capture more angular features in the distribution.

As we have discussed, the spherical harmonic basis may not be an appropriate choice for this decomposition. In the next section, I will present an alternative decomposition which can guarantee that the velocity distribution is everywhere positive and therefore represents a promising and general method for extracting information from directional experiments.

## 4 A discretised velocity distribution

In order to ensure that the velocity distribution is everywhere positive, we propose that the velocity distribution be discretised into $N$ angular components:

$$ \begin{align} & f^1(v) \quad \text{for } \vartheta^j \in [0, \pi/N], \nonumber \\ & \Box f^2(v) \quad \text{for } \vartheta^j \in [\pi/N, 2\pi/N], \nonumber \\ & \Box \dots \nonumber \\ & f^N(v) \quad \text{for } \vartheta^j \in [(N-1)\pi/N, \pi]. \tag{4.1} \end{align} $$
Over each bin in $\vartheta^j$, $f(v)$ has no angular dependence and depends only on a single function of the DM speed. The resulting velocity distribution will be everywhere positive, as long as a suitable parametrisation for the $f^k(v)$ is chosen which is itself everywhere positive.

We consider for simplicity only a discretisation in $\cos \vartheta^j$. We note that in this work we will only be considering the azimuthally-averaged event rate (i.e. integrated over the recoil angle $\varphi$). In this case, we emphasise that no assumptions about the $\varphi^j$-dependence of $f(v)$ are required, as the integral over $\varphi^j$ can be performed exactly, and the azimuthally-averaged rate depends only on the azimuthally-averaged velocity distribution (see appendix A). We can therefore take the velocity distribution to be independent of $\varphi^j$ without loss of generality. However, this analysis could be extended to consider bins in the $\varphi$ angle, with an additional discretisation in $\varphi^j$ if required.

The motivation for this discretised description is that the simplest signal (beyond an isotropic $N=1$ signal) which can be observed with a directional detector is an asymmetry between the event rates in, say, the forward and backward directions. Shortly after the confirmation of a dark matter signal at a directional detector, the number of events may still be quite small (for example, the roughly 10 events required to distinguish from an isotropic background [5], or 30 events required to confirm the peak recoil direction [7]). In this small statistics scenario, constraining a large number of free functions is not feasible. However, if we discretise $f(v)$ into $N=2$ angular components, it may be possible to extract some meaningful directional information with only a small number of events. With larger numbers of events, $N$ can be increased to allow more directional information to be extracted.

We show in figure 1 some examples of this discretised velocity distribution. We show the SHM velocity distribution (top left), as well as the $N=2$ (middle left) and $N=3$ (bottom left) discretised approximations, in all cases integrated over $\varphi^j$. These approximations are obtained by averaging the full velocity distribution over each bin in $\vartheta^j$:

$$f^k(v) = \frac{1}{\cos((k-1)\pi/N) - \cos(k\pi/N)} \int_{\cos(k\pi/N)}^{\cos((k-1)\pi/N)} f(v) \, d\cos\vartheta^j. \quad (4.2)$$

For comparison, the results for a stream distribution are also shown in the right column. We describe these two distribution functions in more detail at the start of section 5.

Upon discretisation, the distribution which is initially focused in one direction (towards $\vartheta^j = 0^\circ$) now becomes constant over $\vartheta^j$ in each bin. For the SHM (left column), there is a predominantly forward-going component ($\cos \vartheta^j > 0$), as well as a smaller backwards component ($\cos \vartheta^j < 0$). For the stream (right column), which has a much narrower dispersion and higher peak speed, the $N=2$ discretised distribution is almost entirely in the forward direction. However, we note that due to the discretisation, there is now a sizeable population of particles with transverse velocities ($\vartheta^j = 90^\circ$), which was not the case in the full distribution. In the $N=3$ discretisation (bottom row), the velocity distributions become more focused and begin to recover more of the directionality of the full distributions. In the limit of large $N$, the full distribution can be recovered.

In order to determine the corresponding event spectrum, we must calculate the Radon transform from this discretised $f(v)$. We introduce the integrated Radon transform (IRT), integrated over the same angular bins as for $f(v)$:

$$\hat{f}^j(v_{\min}) = \int_{\varphi=0}^{2\pi} \int_{\cos(j\pi/N)}^{\cos((j-1)\pi/N)} \hat{f}(v_{\min}, \hat{q}) \, d\cos\vartheta \, d\varphi, \quad (4.3)$$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The full velocity distribution (top) as well as $N = 2$ (middle) and $N = 3$ (bottom) discretised approximations for two examples: the SHM with $v_e = 220$ km s$^{-1}$, $\sigma_v = 156$ km s$^{-1}$ (left column) and a stream with $v_e = 500$ km s$^{-1}$, $\sigma_v = 20$ km s$^{-1}$ (right column). In each case, we have integrated over the $\phi'$ direction and only show $f(v, \cos\theta')$. The vector $v_e$ is aligned along $\theta' = 0$. The same colour scale is used in each plot. In the case of the full stream distribution (top right), large values have been truncated for easier comparison with the other plots.</figcaption></figure>
where $\hat{q} = (\sin \vartheta \cos \varphi, \sin \vartheta \sin \varphi, \cos \varphi)^{\mathrm{T}}$.<fnref n="1" /> In principle, we can define the IRT over a set of bins different from those used to define the discretised distribution function. However, using the same bins in both cases typically simplifies the calculation of the IRT. A notable exception would be to consider fewer bins for the Radon transform than for the distribution function (in order reduce possible discretisation error). In this scenario, calculation of the IRT would be even simpler.

The IRT arises from the calculation of the directional recoil rate, integrated over a given angular bin, which could then be compared with the number of events observed in that bin. We consider the IRT for two reasons. First, it allows us to perform all of the relevant angular integrals in the calculation of the Radon transform, eliminating the need for computationally-intensive numerical integrals over the angular variables. Secondly, the loss of angular information in discretising the velocity distribution means that the full Radon transform of this discretised distribution is unlikely to provide a good fit to the data on an event-by-event basis. Instead, considering the IRT (and correspondingly binned data) should reduce the error which is introduced by using a discretised approximation to the velocity distribution. This in turn allows us to parametrise the v-dependence of each angular bin and mitigate uncertainties in the velocity distribution.

The full details of the derivation of $\hat{f}^j(v_{\min})$ for arbitrary $N$ are included in appendix A. We include this derivation as a pedagogical tool in the interests of anyone who wishes to modify or extend the approach presented here (for example, by considering different definitions for the bins in $\cos \vartheta^j$ or $\cos \vartheta$). However, for the reader interested in simply calculating $\hat{f}^j(v_{\min})$, given a set of $f^k(v)$, the algorithm is given in appendix B. We note that all of the angular integrals involved can be performed analytically, meaning that computing the discretised Radon transform reduces to performing a series of 1-dimensional integrals over $v$ (one integral for each $\hat{f}^j$), which can be performed numerically. A python code which implements the full calculation of appendix B is available from the author.

We note that in the case $N = 1$, the IRT simply reduces to an integral over all angles, exactly recovering the velocity integral which appears in non-directional experiments. Also, in the $N = 2$ case, the IRT takes the following simple form:

$$ \hat{f}^1(v_{\min}) = 4\pi \int_{v_{\min}}^{\infty} \nu \pi f^1(\nu) + \tan^{-1} \frac{\sqrt{1-\beta^2}!}{\beta} f^2(\nu) - f^1(\nu) \, d\nu, \quad (4.4) $$

$$ \hat{f}^2(v_{\min}) = 4\pi \int_{v_{\min}}^{\infty} \nu \pi f^2(\nu) + \tan^{-1} \frac{\sqrt{1-\beta^2}!}{\beta} f^1(\nu) - f^2(\nu) \, d\nu, \quad (4.5) $$

where $\beta = v_{\min}/\nu$. We have also checked using Monte Carlo calculations that our calculations give the correct forms for the IRT in the cases $N = 2$ and $N = 3$. In order to do this, we randomly sample particles from the discretised distributions. For each particle, we simulate a random scattering event and store the values of $v_{\min}$ and $\hat{q}$. We then bin these events in angle and compare with the calculated distribution of $\hat{f}^j(v_{\min})$. The results are in perfect agreement with the results expected from appendix B.

<footnote n="1">¹The integration over $\phi$ essentially reduces the experiment to what is referred to as 1-dimensional directional detection, in which only information about $\theta$ is available. However, it has been shown that this does not significantly reduce the discovery power of the detector when compared with 3-dimensional detection [84].</footnote>
## 5 Comparison with exact results

We now wish to compare these approximate IRTs with the IRTs obtained from the full (non-discretised) velocity distribution. To do this, we select a benchmark velocity distribution (such as the SHM) and calculate the $f^*(v)$ of eq. (4.1) by averaging over $\cos \vartheta^j$ in each bin, as in eq. (4.2) and figure 1. We then insert these into the algorithm presented in appendix B to obtain the corresponding IRTs. We refer to these as the approximate IRTs. For comparison, we use the full Radon transform of eq. (2.8) to obtain the exact IRTs by integrating over $\cos \vartheta$. We fix the peak of the underlying speed distribution to be aligned in the forward direction $\vartheta^j = 0^\circ$. However, we discuss the consequences for ‘misaligned’ distributions in section 6.

We discuss qualitatively the differences in shape between the exact and approximate IRTs, as any such differences could potentially lead to biases in the reconstruction of particle physics parameters and the velocity distribution. However, directional direct detection experiments do not directly measure the Radon Transform, but rather the differential recoil rate $dR/dE_R d\Omega_q$ of eq. (2.1). Analogously to the previous section, we define the directionally-integrated recoil spectrum:

$$ \frac{dR^j}{dE_R} = \int_{\varphi=0}^{2\pi} \int_{\cos(j\pi/N)}^{\cos((j-1)\pi/N)} \frac{dR}{dE_R d\Omega_q} d\cos\vartheta d\varphi. \qquad (5.1) $$

The integrated recoil spectrum is then related to the IRT by

$$ \frac{dR^j}{dE_R} \propto \hat{f}^j(v_{\min}(E_R)) F^2(E_R), \qquad (5.2) $$

where $v_{\min}(E_R)$ is defined in eq. (2.6). The form factor $F^2(E_R)$ leads to a roughly exponential suppression of the rate (relative to the IRT) with increasing $E_R$. Because of this simple proportionality relationship between the IRT and the rate, we do not discuss the shape of the rates in detail.

Instead, we perform a more quantitative comparison between the full and discretised event rates. Assuming a particular idealised detector, we compare the expected number of events $N_j$ obtained within each bin using the exact and approximate IRTs presented in the previous section:<fnref n="2" />

$$ N_j = \int_{E_{\min}}^{E_{\max}} \frac{dR^j}{dE_R} dE_R. \qquad (5.3) $$

In converting from $v_{\min}$ to recoil energy $E_R$, it is necessary to fix the values of the DM and nuclear masses. We therefore take as an example a Fluorine target and DM mass of $m_\chi = 50$ GeV. Direct detection experiments do not probe down to arbitrarily low energies and typically an energy threshold is set below which either the signal may be dominated by unrejected backgrounds or the direction of the recoil cannot be reconstructed. We consider a typical threshold energy of 20 keV [85].

We fix the normalisation of the event rate by requiring a total of 50 signal events across all bins. As previously noted, the anisotropy of a DM signal can be confirmed with around 10 signal events [5], while the median recoil direction can be confirmed to match that expected due to the Earth's motion with of order 30 events [7]. It is reasonable then to expect that

<footnote n="2">²Because the $N = 1$ case (full average over all angles) is exact, the discretised velocity distribution will always recover the correct total number of signal events. We have checked this explicitly up to $N = 10$ by summing the number of events in each angular bin.</footnote>
the approach presented here might be employed once around 50 signal events are observed. In addition, we include a background of 1 event, distributed isotropically. We discuss the impact of a larger number of signal events in section 6.

The error between the exact and approximate results can then be compared with the typical Poissonian uncertainty on the expected number of events in each bin $\Delta N_j \sim \sqrt{N_j}$. If the discrepancy between the exact and approximate results is smaller than this Poisson error in the number of events, we would expect the two methods to be in agreement within statistical uncertainties. This means that the data should be approximately as well-fit by the discretised velocity distribution as by the full distribution. This in turn indicates that the discretised velocity distribution can be used, along with a suitable parametrisation, to fit and extract information about the underlying $f(v)$.

We consider two example distributions, which have already been illustrated in figure 1. The first is the canonical Standard Halo Model (SHM), with parameters $v_e = 220$ km s$^{-1}$ and $\sigma_v = 156$ km s$^{-1}$. We use this not only because it is so often studied in the literature, but also because it has a relatively smooth, simple structure and is not too strongly peaked. If the method presented here is to be useful, it should give accurate results for such a simple underlying distribution. For simplicity we do not truncate the SHM at the Galactic escape speed. The inclusion of such a cut-off would reduce the high-speed DM population, reducing the directionality of the signal and therefore improving the results presented here. In the reconstruction of real data, an escape speed cut-off could be included easily in the chosen parametrisation for $f^*(v)$.

As a comparison, we also consider a stream distribution, modeled using eq. (2.10), but with parameters $v_e = 500$ km s$^{-1}$ and $\sigma_v = 20$ km s$^{-1}$. This leads to a sharply peaked, strongly directional distribution function, which allows us to compare the SHM with a more extreme case. We note in particular that these parameters lead to a sharper distribution (and therefore a more directional rate) than values typically assumed, for example, for the Sagittarius stream [35, 36, 82]. This stream should therefore be considered a 'worst-case' scenario which is difficult to approximate accurately.

## 5.1 N=2 discretisation

Figure 2 shows the comparison between the exact and approximate IRTs for the $N = 2$ discretisation. The exact IRT is shown as a solid line, while the approximate IRT obtained from the discretisation is shown as a dashed line. We show the results for both the SHM (blue) and stream distributions (red). Values of $v_{min}$ in the shaded region lie below the assumed 20 keV energy threshold for a Fluorine target and 50 GeV DM particle. In the case of the SHM, the shape of the Radon transform is well reconstructed in both the forward (left panel) and backward (right panel) directions, even with the crude $N = 2$ discretisation. For the forward IRT ($j = 1$), the ratio between the exact and approximate results is less than 20%, while for the backward ($j = 2$) case, the error is typically larger (50–100%), though the absolute value is smaller.

For the stream distribution, the differences between the exact and approximate results are more significant. The results using the discretised velocity distribution underestimate the forward IRT at low $v_{min}$, while overestimating the backward IRT for all values of $v_{min}$. The reason for this can be seen by examining the velocity distributions illustrated in figure 1. The full stream distribution (top right) is strongly focused in the forward ($\theta^J = 0^\circ$) direction, meaning that particles must scatter through almost 90° to produce a recoil in the backward direction. This is kinematically allowed only for a small fraction of the population. Due to the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Exact (solid) and approximate (dashed) integrated Radon transforms, $f^j$, defined in eq. (4.3), for $N=2$. Results are shown for the SHM (blue) and stream (red) distribution functions. The approximate Radon transforms are obtained by discretising the full velocity distribution into $N$ angular bins. The vector $v_e$ is aligned along $\theta' = 0$. The shaded grey region lies below the energy threshold of 20 keV for a DM mass of 50 GeV.</figcaption></figure>

angular averaging, however, the $N=2$ velocity distribution (middle right) has a significant population of DM particles with velocities at right angles ($\theta' = 90^\circ$) to the forward direction. Thus, the discretised velocity distribution has a greater chance of producing scatters in the backward direction. Overall, the discretised distribution is less focused in the forward direction than the full distribution, resulting in a reduced asymmetry between the forward and backward scattering rates.

Figure 3 shows the comparison between the number of events $N_j$ in each of the angular bins, defined in eq. (5.3), for the case of $N=2$ discretisation. The exact and approximate calculations are shown as solid and hatched bars respectively, while the SHM and stream distributions are shown separately in the left and right panels. In addition, the 'exact' event numbers in each bin are assigned an error bar given by the Poissonian standard deviation $\sqrt{N_j}$.

In the case of the SHM (left), the number of events in the forward direction ($j=1$) is roughly in agreement (within Poisson uncertainties) when calculated using the exact and approximate IRTs. However, in the backward direction ($j=2$), the approximate calculation overestimates the number of events by a factor of 3. As shown in the left hand panel of figure 2, in the forward direction, the greatest discrepancy between the IRTs occurs at low $v_{\min}$. When we consider events only above the energy threshold (i.e. above the shaded grey region), this discrepancy is therefore minimal. In the right hand panel of figure 2, however, the fractional discrepancy grows as a function of $v_{\min}$, and is greatest above the energy threshold of the experiment. This is due to the enhancement in DM particles travelling in the backward direction $\theta' = 180^\circ$ when using the discretised velocity distribution (see figure 1) and leads to the large discrepancy in the number of backward-going events.

This problem is even greater for the stream distribution (right), for which neither the forward nor backward event numbers agree within the statistical errors. In the backwards direction, the exact calculation predicts much less than 1 event, while the approximate calculation overestimates this by a factor of ~30. For such a highly directional distribution, such a small number of bins is clearly not suitable.

We note briefly that in the case of $N=2$ discretisation, the approximate calculation of event numbers gives results which are almost indistinguishable for the SHM and stream distributions. This is because the discretised forms of the SHM and stream distributions
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison of the number of events in each angular bin (eq. (5.3)) for N=2 bins, for a Fluorine target with 20 keV threshold and $m_\chi = 50$ GeV. Event numbers are calculated using the exact (solid bars) and approximate (hatched bars) IRTs. The SHM distribution is shown in the left panel, while the stream distribution is shown on the right. A total of 50 signal events are assumed (as well as 1 isotropically distributed background event). The error bars on the 'exact' number of events are the Poissonian errors $\Delta N_j \sim \sqrt{N_j}$.</figcaption></figure>

differ most substantially at low $v_{\min}$. This can be seen in the middle row of figure 1; eliminating those particles which are below the threshold speed of ~300 km s⁻¹ gives two similar distributions. Given that there is poor agreement between the number of expected events for the *N* = 2 discretisation, and that the discretised SHM and stream distributions are largely indistinguishable, it is clear that such a coarse angular approximation is not suitable for fitting to data.

## 5.2 N=3 discretisation

Figure 4 compares the approximate and exact IRTs for the *N* = 3 discretisation illustrated on the bottom row of figure 1. In the case of the SHM, the agreement between the two functions is improved compared to the *N* = 2 case. The slope of the event rate is degenerate with the DM mass $m_\chi$, so any error in the shape of the IRTs could translate to a bias in the reconstructed value of $m_\chi$. In the forward direction (left panel), the fractional error is reduced to at most 10% and the shape of the IRT appears to be agree closely in all three directional bins. This suggests that any such bias should be minimal.

However, there is still a significant difference between the shapes of the exact and approximate IRTs for the stream distribution. For example, in the forward direction (left panel), the exact calculation predicts scattering events occurring only with a narrow range of $v_{\min}$, from 250 km s⁻¹ to 500 km s⁻¹. This is because almost all particles are moving in the forward direction with speed $v = 500$ km s⁻¹. The kinematics of the scattering requires that $v \cdot \hat{q} = v_{\min}$, meaning that scattering in the forward direction ($\theta = 0°$) is only allowed for $v_{\min} \approx 500$ km s⁻¹. Scattering away from the forward direction, but still included in the first angular bin ($\theta < 60°$), is kinematically allowed for $v_{\min} > \cos(60°) \times 500$ km s⁻¹ $\approx 250$ km s⁻¹. For the discretised velocity distribution, there is a population of DM particles initially directed away from the forward direction ($\theta = 0°$) which can scatter through larger angles and still produce scattering events within the first angular bin. These large-angle scattering events mean that values of $v_{\min}$ down to zero are now kinematically allowed.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** As figure 2, for N = 3 angular bins. In the right panel, the exact Radon transform for the stream is indistinguishable from zero.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** As figure 3, for N = 3 angular bins.</figcaption></figure>

Considering the *j* = 3 angular bin (right panel), we note that the exact stream distribution predicts no scattering events in the backward direction. This is because the particles would have to scatter through an angle much larger than 90° to produce events in the *j* = 3 bin, which is not kinematically allowed. However, for the discretised distribution (bottom right panel of figure 1), it is still kinematically possible for particles on the edge of the forward bin (ϑ^j ≈ 60°) to scatter into the backward direction (ϑ > 120°). Thus, the approximate calculation still predicts a small but non-zero IRT for the *j* = 3 case.

Figure 5 shows the event numbers expected in each of the *N* = 3 angular bins calculated using the exact and approximate IRTs. For the SHM (left panel), in the forward direction (*j* = 1), the approximate calculation underestimates the rate by around 13%, while in the transverse direction (*j* = 2), the approximate calculation overestimates the rate by around 30%. Finally, in the backwards direction (*j* = 3) the rate is overestimated by a factor of 2. However, these discrepancies are roughly the same magnitude as the typical Poisson error in each bin for a sample of 50 events. In contrast to the *N* = 2 case then, the *N* = 3 case can provide a reasonable approximation to the event rate for the SHM.

For the stream distribution (right panel), there are significant discrepancies, which remain larger than the typical statistical error. In particular, there is a substantial leakage of events from the forward direction (*j* = 1) into the transverse direction (*j* = 2) in the approximate calculation. This can be seen in the middle panel of figure 4, in which the exact IRT is
zero above the experimental threshold, while the approximate IRT is substantially non-zero. This in turn arises due to the 'smearing' of the distribution away from the forward direction, as previously discussed. In contrast, in the backward direction (j = 3), the approximate calculation now agrees with the exact calculation, giving zero events. As shown in the right panel of figure 5, the discrepancy between the exact and approximate IRTs arises only at low $v_{min}$ (below the energy threshold of the experiment), where scattering into the backwards bin is kinematically allowed.

### 5.2.1 The folded rate

The $N = 3$ case is important for the scenario where head-tail discrimination is not possible when reconstructing recoil tracks. Head-tail discrimination of tracks has previously been demonstrated [86], but may not be possible with 100% accuracy in near-future detectors [21]. Within the formalism considered here, detectors without sense discrimination cannot distinguish between a recoil with direction $cos \theta$ and another with direction $-cos \theta$. They therefore probe the so-called 'folded' recoil spectrum:

$$ \frac{dR}{dE_R d|cos \theta|} = \frac{dR}{dE_R d cos \theta} + \frac{dR}{dE_R d(-cos \theta)}. \quad (5.4) $$

For the $N = 2$ case, the folded spectrum is entirely isotropic, as the forward and backward IRTs differ only in the sign of $cos \theta$. However, in the $N = 3$ case, the transverse IRT, given by

$$ \hat{f}^{\mathrm{T}}(v_{\mathrm{min}}) = \hat{f}^{2}(v_{\mathrm{min}}) = \int_{\varphi=0}^{2\pi} \int_{-1/2}^{1/2} \hat{f}(v_{\mathrm{min}}, \cos \vartheta) d \cos \vartheta d\varphi, \quad (5.5) $$

is invariant under $\hat{f}(v_{min}, cos \vartheta) \rightarrow \hat{f}(v_{min}, -cos \vartheta)$. That is, the transverse event rate 'folds' back onto itself. Thus, even without sense discrimination, directional experiments will still be sensitive to this transverse scattering rate. By comparison, if the forward and backward directions cannot be distinguished, the remaining two IRTs (the left and right panels in figure 4) are folded together, to obtain the longitudinal Radon transform

$$ \hat{f}^{\mathrm{L}}(v_{\mathrm{min}}) = \int_{\varphi=0}^{2\pi} \int_{-1/2}^{1/2} f(v_{\mathrm{min}}, \cos \vartheta, \varphi) d \cos \vartheta d\varphi + \int_{\varphi=0}^{2\pi} \int_{1/2}^{1} \hat{f}(v_{\mathrm{min}}, \cos \vartheta) d \cos \vartheta d\varphi. \quad (5.6) $$

Analogously, we can define the longitudinal and transverse event spectra, $dR^{LT}/dE_R$, and event numbers $N^{LT}$, derived from the corresponding Radon transforms.

Figure 6 shows the event numbers for this folded rate in the longitudinal and transverse bins. For the SHM (left panel), there is a slight improvement in the agreement between the exact and approximate results compared with figure 5. The backward bin, which was previously overestimated by a factor of 2, has been removed and merged with the forward bin to form the longitudinal bin. Because the approximate result previously underestimated the true number of events in the forward bin, this combination of the bins leads to marginally improved agreement, though the effect is small in this example.

In the case of the stream (right panel), the number of events in each bin appears almost indistinguishable from those of figure 5. This is because the leakage of the events into the backwards direction when using the approximate distribution is minimal and therefore the effect of folding is also minimal. However, we have demonstrated that this method can still be employed when sense discrimination of recoils is not possible and, in general, it should lead to a reduction in the associated discretisation error.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** As figure 3, for the folded event rate with $N = 3$ underlying bins, producing one longitudinal and one transverse directional bin (see eq. (5.4) and associated text).</figcaption></figure>

## 5.3 N=5 discretisation

We now consider a more finely discretised velocity distribution, namely $N = 5$. The comparison of the exact and approximate IRTs is shown in figure 7. The results for the SHM, showing typical discrepancies at the 5% level, indicate that the approximate Radon transform does in fact tend to the true Radon transform in the limit of large $N$. The results for the stream again show much poorer agreement. We expect that the agreement will not improve significantly until the angular size of each bin is close to the angular extent of the underlying distribution function. For the stream, we can see by eye in figure 1 that most of the distribution is distributed over an angle $\theta^j \approx 4 \, 10^\circ$, meaning that roughly $N = 18$ bins are required to prevent the smearing of the distribution function visible for $N = 2$ and $N = 3$. In spite of this, some structures (such as the peak in the forward direction of the upper left panel) are better reconstructed than in the $N = 3$ case, indicating that additional information is still gained by adding more bins. In principle, there is no obstacle to increasing the number of bins up to (and beyond) $N = 18$, as the formalism of appendix B is applicable for arbitrary $N$. However, the appropriate number of bins would typically be dictated by the amount of data available, with a large amount of data required to justify fitting the number of parameters associated with $N = 18$ bins.

We note that for the bin focused in the backward direction, $j = 5$ (bottom right), the exact and approximate IRTs for the stream are indistinguishable and very close to zero. As in the $N = 3$ case, scattering in the backward direction is not possible for the full stream distribution because of the strong directionality. For the discretised distribution, most of the distribution is focused in the forward-directed bin ($k = 1$ in eq. (4.2)). In contrast to the $N = 3$ result, however, this is sufficiently directional that scattering into the backward direction (into the $j = 5$ bin) is no longer kinematically allowed, resulting in significantly closer agreement in the backward direction.

In figure 8, we show the event numbers in each bin for the $N = 5$ discretisation. For the SHM (left panel), there is now close agreement between the exact and approximate results, with less than 10% discrepancy for the forward bins $j = 1, 2$ and an error between 20% and 40% in the transverse and backward bins $j = 3, 4, 5$. For all bins, the agreement between the exact and approximate results is significantly smaller than the Poisson uncertainty for the example of 50 signal events. This closely matches the expectation from figure 7, which shows close agreement between the exact and approximate IRTs.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Exact (solid) and approximate (dashed) integrated Radon transforms, $\hat{f}^j$, defined in eq. (4.3), for $N=5$. Results are shown for the SHM (blue) and stream (red) distribution functions. The approximate Radon transforms are obtained by discretising the full velocity distribution into $N$ angular bins. The vector $v_e$ is aligned along $\theta' = 0$. In the lower left panel, the exact IRT for the stream is indistinguishable from zero, while in the lower right panel, both the approximate and exact IRTs are indistinguishable from zero.</figcaption></figure>

For the stream (right panel), there is an improved fit between the exact and approximate results, with bins $j = 1, 2, 4, 5$ showing agreement within the statistical uncertainty. Leakage of events between the two forward bins leads to a slight overestimation of the $j = 2$ rate. However, the biggest problem remains leakage of events into the transverse bin $j = 3$ which gives a significant overestimation of the number of events when using the approximate IRT. Once again, this is seen easily in the upper right panel of figure 8, where the approximate IRT is non-zero above the energy threshold of the experiment.

We note that while there remains a notable discretisation error for the stream distribution, it is now possible to distinguish the results for the SHM and stream distribution when using the approximate IRTs, in contrast to $N = 2$ discretisation. For example, in the $j = 3$ bin, the approximate calculation predicts around ~8 events for the SHM, compared to ~2 events for the stream. This difference is larger than the error induced by the discretisation (a discrepancy of ~2 events for both the SHM and stream). This indicates that while the $N = 5$ discretisation may not be able to accurately reproduce the event rate for a stream distribution, it may still be useful for distinguishing between different distributions. For example, if the $N = 5$ discretisation were used to fit to data and the results indicated a smaller number of events than expected in the $j = 3$ bin, this could point to deviations from the SHM.

## 5.4 Forward-backward asymmetry

Finally, we calculate the forward-backward asymmetry in the number of events, $A_{FB}$, and compare the exact and approximate results. In terms of the number of forward and backward
<figure><img src="image_11.png" /><figcaption>Figure 8. As figure 3, for N=5 angular bins.</figcaption></figure>

events between the threshold and maximum energies $E_{\text{th}}$ and $E_{\text{max}}$,

$$N_F = \int_{E_{\text{max}}} \int_{2\pi} \int_{1} \frac{dR}{dE_R d\Omega_q} d\cos\vartheta d\phi dE_R \quad (5.7)$$

$$N_B = \int_{E_{\text{max}}}^{E_{\text{th}}} \int_{2\pi}^{\varphi=0} \int_{0}^{0} \frac{dR}{dE_R d\Omega_q} d\cos\vartheta d\phi dE_R, \quad (5.8)$$

the asymmetry is given by

$$A_{\text{FB}} = \frac{N_F - N_B}{N_F + N_B}. \quad (5.9)$$

By comparing the values of $A_{\text{FB}}$ obtained from the full and discretised distributions, we can obtain a measure of how well the discretised distribution captures the directionality of the signal, as well as showing how this improves with increasing $N$. In addition, this forward-backward asymmetry would potentially be the first directional signal measured with directional detectors and so determining how well it can be recovered is of substantial importance. Using the current formalism, we can consider only even values of $N$, for which each angular bin lies entirely in either the forward or backward direction. In principle, the calculation can also be extended to include odd values of $N$ using the framework of appendix A.

Figure 9 shows the forward-backward asymmetry obtained from the discretised velocity distribution for $N = 2, 4, 6, 8, 10$ (filled squares), compared to the true forward-backward asymmetry, obtained from the full velocity distribution (solid lines). In the case of the stream (red), the event rate is strongly asymmetric ($A_{\text{FB}} \approx 1$) as the velocity distribution is highly focused in the forward direction. By comparison, the SHM (blue) has a smaller asymmetry ($A_{\text{FB}} \approx 0.9$) due to its wider velocity dispersion. Dotted lines show the $1\sigma$ statistical uncertainties on the measured value of $A_{\text{FB}}$ assuming 50 signal events.

For both the SHM and stream, the asymmetry is significantly underestimated when using only the simple forward-backward ($N = 2$) discretisation. This matches the analysis of the previous sections, which demonstrated that the discretised velocity distribution tends to lead to a leakage of events from the forward to the backward direction. However, for both distributions, the asymmetry obtained from the discretised distribution rapidly approaches the true value with increasing $N$. For the stream distribution, the approximation converges
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Forward-backward asymmetry of the event rate $A_{FB}$, defined in eq. (5.9), using the full velocity distribution (solid lines) and the discretised velocity distribution (filled squares). Results are shown for the SHM and stream distributions in blue and red respectively. Dotted lines show the $1\sigma$ statistical uncertainties on the measured value of $A_{FB}$ assuming 50 signal events.</figcaption></figure>

more rapidly than for the SHM. This is because of the strong directionality of the stream velocity distribution, so only relatively few bins are required to capture a simple forward-backward asymmetry. Even with only 4 angular bins, the error in $A_{FB}$ for both velocity distributions is less than 10% and lies within the $1\sigma$ statistical uncertainty expected from a sample of 50 signal events.

# 6 Discussion

We have demonstrated that for smooth SHM-like distributions, the discretised velocity distribution allows us to obtain a good approximation to the integrated Radon transform (IRT) with relatively few angular bins. Though the $N = 2$ discretisation is too simple an approximation, for 3 or more bins the discretised distribution leads to a number of events in each bin which matches the true number within statistical uncertainties (assuming 50 signal events). This means that we should be able to parametrise the speed distributions in each angular bin $f^*(v)$ and account for astrophysical uncertainties without introducing a large error in the number of observed events. For comparison, we have also considered a more extreme and highly-directional stream distribution. The disagreement between the true recoil spectrum and that obtained from the discretised velocity distribution is significantly larger in this case. Though increasing the number of bins $N$ improves the agreement, a rather large number of bins ($N \sim 18$) may be needed to accurately describe the stream distribution.

The results we have presented here have assumed an idealised future directional detector. However, application of these results is not necessarily as straightforward in a realistic experiment. Conversion from $v_{min}$ to $E_R$ requires knowledge of the DM mass $m_\chi$ and attempting to fit both the velocity distribution and mass simultaneously can lead to spurious results [40, 42]. Furthermore, for comparison with experimental data, we must take into account the fact that experiments have finite angular resolution, typically in the range $20^\circ-80^\circ$,
with higher resolution at high energy [21]. Finally, throughout this work, we have assumed
that the basis for discretising the velocity distribution is aligned with the peak of the under-
lying velocity distribution (i.e. $v_e$ aligned along $\theta^J = 0$). However, this is not known a priori
and a reliable method of selecting the orientation of the angular basis is required.

In spite of these open issues, the work presented here is general and conservative. We
have demonstrated that the stream distribution is more poorly fit than the SHM distribution
with few angular bins. However, the stream is an extreme example, as discussed in section 5,
and it is unlikely the discretisation error will be any greater than observed in that case.
Other possibilities for the velocity distribution would tend to be less directional than the
stream example, leading to better agreement between the exact and approximate rates. For
a dark disk, the low value of $v_e$ means that the velocity distribution appears almost isotropic
in the lab frame, which will only reduce the angular discretisation error. A misalignment
between $v_e$ and $\theta^J = 0^\circ$ would also lead to a distribution more isotropic in $\theta^J$. Similarly,
including finite angular resolution will smear the velocity distribution (in a style similar
to the discretisation), reducing the directional asymmetry and therefore the discretisation
error. We have also demonstrated that this method is applicable in detectors where sense
recognition is not possible.

For a realistic analysis, all of these issues should be taken into account. Based on mock
(or future) data, the optimal direction for $\theta^J = 0^\circ$ should be selected (potentially based on
the observed median recoil direction [5, 7]). An appropriate number of bins $N$ should then be
chosen, which may be influenced by the angular resolution of the experiment and the number
of observed events. In some scenarios, such as for strongly peaked signals, deviations from
the smooth SHM may be obvious from the data. In this case, a definition of the angular bins
different from that used here (including differently-sized bins) may be optimal to reduce the
discretisation error. Using the results of appendix A, such a scenario can be straightforwardly
accommodated, although we leave the issue of choosing the number and size of the bins to
future studies.

In order to make more quantitative statements, we have focused on a specific scenario,
in which 50 signal events (and 1 background event) are observed. This has allowed us to
compare the typical Poisson uncertainties on the number of observed events with the error
induced by using the discretised velocity distribution. Though a signal consisting of 50 events
is somewhat arbitrary, it is representative of the lower limit on the number of events required
before such an analysis would be reasonable. As the number of signal events increases, the
typical statistical uncertainty decreases. Eventually, this uncertainty becomes smaller than
the error between the approximate and exact results. Increasing the number of angular bins
would then reduce the discretisation error and increase the statistical uncertainty on the
number of events in each bin, reducing any bias which may be induced by the discretisation.
This therefore provides a natural scheme for deciding how many angular bins are required,
based on the number of signal events observed.

Once the number of bins $N$ has been fixed and a suitable parametrisation for the radial
functions $f^k(\mathbf{v})$ has been chosen, the full parameter space (of both DM particle physics
parameters and astrophysics parameters) should be fit to the data. We note in particular
that by summing the radial functions $f^k(\mathbf{v})$, the 1-dimensional speed distribution is recovered,
meaning that non-directional experiments can easily be included in this framework. As
emphasised in this section, the fitting process will be highly dependent on the underlying
DM and experimental parameters, as well as on the chosen parametrisation for $f^k(\mathbf{v})$. We
have therefore focused here on the discretisation of $f(\mathbf{v})$ and left more involved investigations
for future work.
## 7 Conclusions

In this work, we have presented an angular basis which can be used to parametrise the DM velocity distribution. This involves discretising the velocity distribution into N angular bins:

$$f(v) = f(v, \cos \vartheta^j, \varphi^j) = \begin{cases} f^1(v) & \text{for } \vartheta^j \in [0, \pi/N], \\ \Box f^2(v) & \text{for } \vartheta^j \in [\pi/N, 2\pi/N], \\ \cdot & \\ \Box f^k(v) & \text{for } \vartheta^j \in [(k-1)\pi/N, k\pi/N], \\ \cdot & \\ \Box f^N(v) & \text{for } \vartheta^j \in [(N-1)\pi/N, \pi]. \end{cases} \quad (7.1)$$

In appendices A and B, we have provided recipes for calculating the corresponding directionally-integrated Radon transforms (IRTs), which appear in the directional event rate. Alternative methods for parametrising the DM velocity distribution (such as refs. [59, 60], and in particular those based on spherical harmonics) may lead to negative values $f(v)$, leading to unphysical distribution functions. It is not clear what affect this might have on parameter reconstruction or how this problem can be mitigated. The advantage of the basis presented here is that it guarantees that the resulting distribution function is everywhere positive and therefore physical.

We have investigated the possible size of discretisation errors in the directional event rate by comparing the IRT in each bin obtained from the discretised and full velocity distributions. The simplest possible discretisation ($N = 2$) is unsuitable for fitting to data, as it substantially underestimates the forward scattering rate, while underestimating the backwards scattering rate. However, starting from $N = 3$, the standardly-assumed SHM can be well fit by the discretised velocity distribution. We have shown that once the DM nature of the signal is confirmed (with around 50 events), the $N = 3$ discretisation allows us to calculate the number of events in each angular bin and obtain good agreement with the true event numbers (within statistical uncertainties). Increasing the number of bins to $N = 5$ improves this agreement, giving discretisation errors in the range 10–50% depending which angular bin is considered.

For comparison, we also consider an extreme, highly directional stream distribution. In this case, the discretisation error is much larger and exceeds the statistical uncertainties in the $N = 5$ case. An estimated $N = 18$ angular bins would be required to accurately map the full velocity distribution. However, increasing $N$ does reduce the discrepancy between the exact and approximate results. In addition, with as few as $N = 5$ bins, it is possible to distinguish between the SHM and stream distribution. This suggests that even though the discretised distribution gives a poor approximation to the stream, it may be useful in distinguishing different distributions, even for relatively few bins.

Finally, we have considered the forward-backward asymmetry $A_{FB}$ in the event rate, calculated using both the full velocity distribution and the discretised velocity distribution. The latter calculation rapidly converges to the correct value with increasing $N$. For both the SHM and stream distributions, the error in $A_{FB}$ is below 10% for $N \ge 4$ bins, within the statistical error expected for 50 signal events. We have also demonstrated that the method presented here can still be used when head-tail discrimination of the recoil tracks is not possible, in which case $N \ge 3$ bins are required.
In this work, we have only considered discretising the angular component of the DM velocity distribution, leaving the N radial functions $f^k(v)$ fixed to their 'true' values. In future work, it will be necessary to combine the discretisation presented here with a parametrisation of these radial functions (such as the parametrisations of refs. [40, 42, 44]), in order to reconstruct the velocity distribution (and other DM parameters) from mock data. In addition, several questions are still to be addressed, including how to decide on the optimal 'forward' direction for the discretisation, and how realistic angular resolution would impact the results presented here. However, these additional considerations are only expected to improve the agreement between the true and approximate event rates. This work represents an initial step towards a parametrisation of the full DM velocity distribution, which would allow future data from directional experiments to be analysed without astrophysical uncertainties and which would potentially allow the velocity distribution itself to be probed and reconstructed.

## Acknowledgments

The author thanks Anne M. Green, Julien Billard and Samuel K. Lee for helpful discussions on the subject of directional detection. The author also thanks Anne M. Green and Ciaran A.J. O'Hare for critical comments on this manuscript. The author is supported by the European Research Council (ERC) under the EU Seventh Framework Programme (FP7/2007–2013)/ERC Starting Grant (agreement No. 278234 – 'NewDark' project).

## A Calculating the Radon transform

In this appendix, we derive the Radon transform corresponding to the discretised velocity distribution. The final result can be found in appendix B. The Radon transform of $f(v)$ is defined by:

$$ \hat{f}(v_{\min}, \hat{q}) = \int_{\mathbb{R}^3} f(v) \delta(v \cdot \hat{q} - v_{\min}) \, d^3 v. \qquad (A.1) $$

We write the coordinates of the velocity and recoil momentum as

$$ \begin{gathered} \mathbf{v} = v \sin \vartheta^j \cos \varphi^j, \sin \vartheta^j \sin \varphi^j, \cos \vartheta^j \\ \hat{\mathbf{q}} = (\sin \vartheta \cos \varphi, \sin \vartheta \sin \varphi, \cos \vartheta), \end{gathered} \qquad (A.2) $$

and consider here only the azimuthally integrated Radon transform. We write this as $\hat{f}(v_{\min}, \cos \vartheta)$, which is given by

$$ \begin{align} \hat{f}(v_{\min}, \cos \vartheta) &= \int_{2\pi} f(v_{\min}, \cos \vartheta, \varphi) \, d\varphi \nonumber \\ &= \int_{0}^{2\pi} \int_{0}^{\pi} f(v) \delta(v \cdot \hat{q} - v_{\min}) \, d^3 v \, d\varphi \nonumber \\ &= \int_{0}^{2\pi} \int_{0}^{\pi} f(v) \delta(v \cdot \hat{q} - v_{\min}) \, d\varphi \, d^3 v \tag{A.3} \\ &= \int_{0}^{2\pi} \int_{0}^{\pi} f(v) D(v_{\min}, \cos \vartheta, v) \, d^3 v. \nonumber \\ &\equiv \end{align} $$

We expand the $\delta$-function explicitly in terms of the angular coordinates:

$$ \delta(v \cdot \hat{q} - v_{\min}) = \frac{1}{v} \delta \left( \sin \vartheta \sin \vartheta^j \cos(\varphi - \varphi^j) + \cos \vartheta \cos \vartheta^j - v_{\min}/v \right) \equiv \frac{1}{v} \delta(g(\varphi)). \qquad (A.4) $$
We rewrite the argument of the $\delta$-function as a function $\varphi$:

$$ \delta(g(\varphi)) = \sum_{i} \frac{\delta(\varphi - \varphi_i)}{|g^j(\varphi_i)|} \qquad (A.5) $$

Here, we sum over those values of $\varphi_i$ satisfying $g(\varphi_i) = 0$:

$$ \cos(\varphi_i - \varphi^j) = \frac{\beta - \cos \vartheta \cos \vartheta^j}{\sin \vartheta \sin \vartheta^j} \equiv \alpha, \qquad (A.6) $$

where we have also defined $\beta = v_{\min}/v$. The solutions for $\varphi \in [0, 2\pi]$ are:

$$ \begin{aligned} \varphi_1 &= \varphi^j + \cos^{-1} \alpha, & \text{for } \varphi^j &\in (0, 2\pi) \text{, } \cos^{-1} \alpha \in \mathbb{R} \\ \varphi_2 &= \varphi^j + 2\pi \cos^{-1} \alpha, & \text{for } \varphi^j &\in (0, \cos^{-1} \alpha) \\ \varphi_3 &= \varphi^j + \cos^{-1} \alpha, & \text{for } \varphi^j &\in (2\pi, \cos^{-1} \alpha, 2\pi) \\ \varphi_4 &= \varphi^j - \cos^{-1} \alpha, & \text{for } \varphi^j &\in \cos^{-1} \alpha, 2\pi \end{aligned} \qquad (A.7) $$

We note that these solutions exist only for $\beta \in [0, 1]$ (or equivalently $v > v_{\min}$) and for $\alpha \in [-1, 1]$, otherwise eq. (A.6) cannot be satisfied. If these constraints are satisfied, there exist exactly 2 solutions for a given value of $\varphi^j$ and therefore 2 $\delta$-functions in eq. (A.5).

For the derivative of $g(\varphi)$ we obtain

$$ g'(\varphi) = - \sin \vartheta \sin \vartheta^j \sin(\varphi - \varphi^j). \qquad (A.8) $$

Substituting the values of $\varphi_{1,2,3,4}$, we see that

$$ |g^j(\varphi_{1,2,3,4})| = (\sin \vartheta \sin \vartheta^j)^2 - (\beta - \cos \vartheta \cos \vartheta^j)^2. \qquad (A.9) $$

Each of the two $\delta$-functions therefore contributes the same amount to the integral, regardless of the value of $\varphi^j$. We can now perform the integral over $\varphi$:

$$ \begin{aligned} D(v_{\min}, \cos \vartheta, \mathbf{v}) &= \frac{1}{v} \int_0^{2\pi} \delta(g(\varphi)) \, d\varphi \\ &= -\frac{2C(\alpha)}{v \left(\sin \vartheta \sin \vartheta^j\right)^2 - \left(\beta - \cos \vartheta \cos \vartheta^j\right)^2} \Theta(v - v_{\min}) \qquad (A.10) \\ &= \frac{2}{v} C(\alpha) I(\beta, \cos \vartheta, \cos \vartheta^j) \Theta(v - v_{\min}) \end{aligned} $$

where $C(\alpha) = 1$ for $\alpha \in [-1, 1]$ and vanishes otherwise. We note that $D$ is independent of the azimuthal angle $\varphi^j$. The constraint $\alpha \in [-1, 1]$ is satisfied for $\cos \vartheta^j \in [x_-, x_+]$, where

$$ x_{\pm} = \beta \cos \vartheta \pm \sqrt{1 - \beta^2 \sin \vartheta}. \qquad (A.11) $$

We therefore obtain

$$ \hat{f}(v_{\min}, \cos \vartheta) = 2 \int_{x_{-}}^{x_{+}} \int_{v_{\min}}^{\infty} f(v, \cos \vartheta^j) I(\beta, \cos \vartheta, \cos \vartheta^j) v \, dv \, d \cos \vartheta^j, \qquad (A.12) $$
where we have performed the $\varphi^j$ integral over the velocity distribution, and defined

$$f(v, \cos \vartheta^j) = \int_0^{2\pi} f(v, \cos \vartheta^j, \varphi^j) d\varphi^j, \quad (A.13)$$

because $l$ does not depend on $\varphi^j$. We note from this that the azimuthally-integrated Radon transform is unaffected by the exact $\varphi^j$-dependence of $f(v)$. Instead, the azimuthally-integrated Radon transform depends only on the azimuthally-integrated velocity distribution. For the framework considered here, then, we can assume that $f(v)$ is independent of $\varphi^j$ without loss of generality.

We will consider a velocity distribution which is discretised into $N$ angular pieces:

$$ \begin{array}{ll} f^1(v) & \text{for } \vartheta^j \in [0, \pi/N], \\ \Box f^2(v) & \text{for } \vartheta^j \in [\pi/N, 2\pi/N], \\ \hline f(v) = f(v, \cos \vartheta^j, \varphi^j) = & \cdot \\ \Box f^k(v) & \text{for } \vartheta^j \in [(k-1)\pi/N, k\pi/N], \\ \hline \Box f^N(v) & \text{for } \vartheta^j \in [(N-1)\pi/N, \pi]. \end{array} \quad (A.14) $$

We would then ultimately like to calculate the directionally-integrated Radon transform $\hat{f}^j$, integrated over the same bins in $\cos \vartheta$ as for the velocity distribution,

$$ \begin{aligned} \hat{f}^j(v_{\min}) &= \int \cos((j-1)\pi/N) f(v_{\min}, \cos \vartheta) d\cos \vartheta & (A.15) \\ &= 2 \int_{a_j}^{\cos(\pi/N)} \int_{v_{\min}}^{x_{-1}} \int_{x_{-1}}^{x_{+}} f(v, \cos \vartheta^j) I(\beta, \cos \vartheta, \cos \vartheta^j) d\cos \vartheta^j \, dv \, dv \, d\cos \vartheta, \end{aligned} $$

Here, we have introduced the notation $a_k = \cos(k\pi/N)$ and defined $f(v, \cos \vartheta^j) = 2\pi f(v, \cos \vartheta^j, \varphi^j)$. We divide the integration range for $\cos \vartheta^j$ into different regions:

$$ \begin{aligned} \cos \vartheta^j &\in [x_{-}, a_{k_{-1}}] \\ \cos \vartheta^j &\in [a_{k_{-1}}, a_{k_{-2}}] \\ \cos \vartheta^j &\in [a_{k_{+1}}, a_{k_{+}}] \\ \cos \vartheta^j &\in [a_{k_{+}}, x_{+}], \end{aligned} \quad (A.16) $$

where $k_{\pm}$ are defined such that:

$$ x_{\pm} \in [a_{k_{\pm}}, a_{k_{\pm}-1}]. \quad (A.17) $$

These regions are chosen such that $f(v, \cos \vartheta^j)$ is independent of $\cos \vartheta^j$ within each region. We note that if $k_+ = k_-$, there is only one region: $\cos \vartheta^j \in [x_-, x+]$. Using these definitions
we can therefore rewrite the term in brackets in eq. (A.15) as

$$
\begin{align}
& \int_{x_{-}}^{x_{+}} \left( f(v, \cos \vartheta)^{J} I(\beta, \cos \vartheta, \cos \vartheta)^{J} d \cos \vartheta^{J} \right. \nonumber \\
& \left. \qquad - \int_{x_{-}}^{x_{+}} \left( k_{-}(v) \int_{x_{-}}^{x_{+}} I(\beta, \cos \vartheta, \cos \vartheta')^{J} d \cos \vartheta' \right) \right. \nonumber \\
& \left. \qquad = 2\pi f \sum_{k=0}^{k_{-1}} \int_{x_{-}}^{x_{+}} \left( \begin{aligned}[t]
& 2\pi f^{k}(v) \int_{a_{k-1}}^{a_{k}} I(\beta, \cos \vartheta, \cos \vartheta')^{J} d \cos \vartheta' \\
& + 2\pi f^{k+1}(v) \int_{a_{k+1}}^{a_{k}} I(\beta, \cos \vartheta, \cos \vartheta')^{J} d \cos \vartheta'
\end{aligned} \right) \right. \tag{A.18} \\
& \qquad \left. + \sum_{k=k_{-1}+1}^{k_{+}} \int_{x_{+}}^{x_{-}} \left( \begin{aligned}[t]
& 2\pi f^{k}(v) \int_{a_{k}}^{a_{k+1}} I(\beta, \cos \vartheta, \cos \vartheta')^{J} d \cos \vartheta' \\
& + 2\pi f^{k+1}(v) \int_{a_{k+1}}^{a_{k}} I(\beta, \cos \vartheta, \cos \vartheta')^{J} d \cos \vartheta'
\end{aligned} \right) \right). \nonumber
\end{align}
$$

We can now explicitly perform the integral over cos $\vartheta^J$,

$$
\begin{align}
& \int_{a_{k-1}}^{a_k} I(\beta, \cos \vartheta, \cos \vartheta^J) d \cos \vartheta^J \nonumber \\
&= \int_{\text{a}_k}^{\text{a}_{k-1}} \frac{d \cos \vartheta^J}{(\sin \vartheta \sin \vartheta^J)^2 - (\beta - \cos \vartheta \cos \vartheta^J)^2} \nonumber \\
&\stackrel{\text{!!, a}_k}{=} \sin^{-1} \frac{\sqrt{1 - \beta^2 \sin \vartheta}}{\vartheta - \beta \cos \vartheta} \quad \text{a}_k \nonumber \\
&= F(\text{a}_{k-1}, \cos \vartheta) - F(\text{a}_k, \cos \vartheta), \tag{A.19}
\end{align}
$$

where we have defined

$$
F(y, \cos \vartheta) = \sin^{-1} \frac{y - \beta \cos \vartheta}{\sqrt{1 - \beta^2 \sin^2 \vartheta}}. \quad (A.20)
$$

We note that

$$
F(x_{\pm}, \cos \vartheta) = \pm \frac{\pi}{2}, \qquad (A.21)
$$

for all $\cos \vartheta$. We also note that $F(y, \cos \vartheta)$ can be integrated analytically:

$$
\begin{align}
& \int_x F(y, \cos \vartheta) \cos \vartheta \nonumber \\
&= x \sin^{-1} \frac{\sqrt{y - \beta x}}{1 - x^2 \sqrt{1 - \beta^2}}! \nonumber \\
&\quad + y \tan^{-1} \frac{x - \sqrt{y\beta}}{t} \tag{A.22} \\
&\quad + 2 \tan^{-1} \frac{1 - y^2}{-2 \sqrt{\beta^2 - \beta x^2}} \frac{y\beta(1+x)}{y - \beta x \sqrt{1 - \beta^2}} \nonumber \\
&\quad - 2 \tan^{-1} \frac{1 - y^2 - \beta^2 + x - y\beta(1-x)}{(y+\beta)t} \nonumber \\
&\equiv J(x, y) + C \tag{A.23}
\end{align}
$$

where it is understood that $J(x, y)$ is a function of $\beta$ and where

$$
t = (1 - x^2)(1 - \beta^2) - (y - \beta x)^2. \tag{A.23}
$$
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Integration limits for eq. (A.15) as a function of cos $ϑ$, for a fixed value of $\theta$. The limits $x_+$ and $x_-$ are shown as solid and dashed lines respectively. The dotted vertical and horizontal lines mark the values $\cos(n\pi/N)$ for $n = 0, 1, 2, 3$ and $N = 3$. If we wish to calculate $\hat{f}^1(v_{\min})$, we must perform the angular integral over the shaded region.</figcaption></figure>

It is now possible to also complete the integration over $\cos \vartheta$. Assuming that $k_{\pm}$ do not change over a range $\cos \vartheta \in [b_i, b_{i+1}]$, we can combine eqs. (A.18) and (A.22) to obtain

$$ 
\begin{aligned}
& \int_{b_{i+1}}^{b_i} \int_{x_{-i}}^{x_i} f(v, \cos \vartheta^j) I(\theta, \cos \vartheta, \cos \vartheta^j) d \cos \vartheta^j d \cos \vartheta \\
&= 2\pi f^{k-}(v) \frac{\pi}{2} (b_{i+1} - b_i) + J(b_{i+1}, a_{k-1}) - J(b_i, a_{k-1}) \\
& \quad + \sum_{k=k-1}^{k=1} 2\pi f^k(v) (J(b_{i+1}, a_{k-1}) - J(b_i, a_{k-1}) - J(b_{i+1}, a_k) + J(b_i, a_k)) \\
& \quad + 2\pi f^{k+}(v) \frac{\pi}{2} (b_{i+1} - b_i) - J(b_{i+1}, a_{k+}) - J(b_i, a_{k+})
\end{aligned}
\quad (A.24) $$

However, we must take into account the fact that the values of $x_{\pm}$ and therefore of $k_{\pm}$ depend on $\theta$ and $\cos \vartheta$. We now discuss how to divide up the integration regions of $v$ and $\cos \vartheta$ such that we can apply eq. (A.24).

As an illustration, we show in figure 10 the values of $x_+(x_-)$ as solid (dashed) lines as a function of $\cos \vartheta$, for a fixed value of $\theta$. In evaluating eq. (A.15), we wish to integrate over the region enclosed by the solid and dashed lines, for the relevant range of $\cos \vartheta$. As an example, we show with horizontal and vertical dotted lines the edges of the discretised regions in $\cos \vartheta$ and $\cos \vartheta^j$ for $N = 3$. That is, the dotted lines show the values $\cos(n\pi/N)$ for $n = 0, 1, 2, 3$. If we wish to evaluate $\hat{f}^1(v_{\min})$, we need to integrate over the shaded region in figure 10.

In figure 11, we show $x_{\pm}$ for different values of $\theta$. We give two examples, for two different ranges of $\theta$: $\theta < 1/2$ (black) and $\theta > 3/2$ (blue). The values of $k_{\pm}$ (and therefore the relevant $f^k(v)$) will clearly depend on the value of $\theta$, so we must be careful to account for this properly.

Calculation of the values of $k_{\pm}$ and the corresponding bin edges in $\cos \vartheta$ ($\{b_i\}$) is involved, though not technically difficult. We therefore sketch the procedure here. We begin with the definition of $x_+$:

$$ x_+ = x_+(\cos \vartheta) = \theta \cos \vartheta + \sqrt{1 - \theta^2 \sin \vartheta}. \quad (A.25) $$
<figure><img src="image_14.png" /><figcaption>**Figure 11.** As figure 10, but for two different values of $\theta$. The black lines show $x_{\pm}$ for a value of $\theta$ in the range $\theta < 1/2$ while blue lines show $x_{\pm}$ for a value in the range $\theta > 3/2$.</figcaption></figure>

By straight-forward differentiation, we observe that this function has a maximum at $\cos \vartheta = \theta$. For concreteness, we consider $\theta$ in the range,

$$ \theta \in [\cos(n\pi/N), \cos((n-1)\pi/N)] \quad \text{for } n = 1, 2, \dots, N/2. \qquad (A.26) $$

We must consider three separate cases, depending on $j$ (i.e. which $\cos \vartheta$ bin we are considering):

(i) $j \le (n-1)$

In this case, $\cos \vartheta > \theta$, meaning that $\partial x_+/\partial \cos \vartheta < 0$, so $x_+$ is monotonically decreasing in $\cos \vartheta$.

(ii) $j \ge (n+1)$

In this case, $\cos \vartheta < \theta$, meaning that $\partial x_+/\partial \cos \vartheta > 0$, so $x_+$ is monotonically increasing in $\cos \vartheta$.

(iii) $j = n$

In this case, the maximum of $x_+$ lies in the $j^{th}$ bin.

For each value of $j$, we can now determine the maximum and minimum value of $x_+$ within that $\cos \vartheta$ bin by using the definition in eq. (A.25). This allows us to determine the value of $k_+$ (i.e. the $\cos \vartheta^j$ bin into which $x_+$ falls). In fact, it can be shown that for each value of $j$ (apart from $j = n$), $x_+$ crosses one of the $\cos \vartheta^j$ bin edges exactly once, when $\cos \vartheta = \cos \vartheta^j_+$, satisfying

$$ x_+(\cos \vartheta^j_+) = \cos(n\pi/N), \qquad (A.27) $$

for some $l = 1, 2, \dots, N$. This means that for a given $j$, there are two distinct regions in $\cos \vartheta$ ($\cos \vartheta < \cos \vartheta^j_+$ and $\cos \vartheta > \cos \vartheta^j_+$) each with a different value of $k_+$.

We can perform a similar analysis for $x_-$ to obtain the values of $k_-$ for a given $n$ and $j$, as well as the values of $\cos \vartheta^j_-$, where the curve $x_-$ crosses from one bin to the next. However, we can only apply eq. (A.24) if both $k_+$ and $k_-$ do not change over the range $\cos \vartheta \in [b_i, b_{i+1}]$.
It is clear then that we must subdivide each $\cos \vartheta$ bin into 3 smaller bins. The edges of these sub-bins will be $\cos \vartheta^j_+$ and $\cos \vartheta^j_-$, at which the value of either $k_+$ or $k_-$ changes.

Finally, we need to determine which of these crossings $\cos \vartheta_{\pm}^j$ occurs first. The crossings occur at the same value of $\cos \vartheta$ if $\cos \vartheta^j_+ = \cos \vartheta^j_-$ which occurs for $\beta = \cos((n - 1/2)\pi/N)$. For each value of $n$, we therefore need to distinguish two regimes, depending on which of the crossings $\cos \vartheta_{\pm}^j$ is larger. It is helpful then to further subdivide the range of $\beta$, writing

$$ \beta \in [\cos(m\pi/(2N)), \cos((m - 1)\pi/(2N))] \quad \text{for } m = 1, 2, \dots, N, \quad (A.28) $$

from which we obtain

$$ n = \begin{cases} (m+1)/2 & \text{for } m \text{ odd} \\ m/2 & \text{for } m \text{ even.} \end{cases} \quad (A.29) $$

We can now determine which of $\cos \vartheta_{\pm}^j$ is greater, depending on whether $m$ is even or odd. For a given value of $m$ then, the values of the sub-bin edges in $\cos \vartheta$ ($\{b_i\}$) are now well determined, along with the values of $k_{\pm}$ within each sub-bin. This allows us to apply eq. (A.24) to each of these sub-bins.

Finally, we can perform the integral over $v$ described in eq. (A.15). However, we note that we should sub-divide the range of integration, depending on the value of $m$. That is, the integral over $v$ becomes,

$$ \int_{v_{\min}}^{\infty} \int_{m=1}^{\infty} \int_{v_{\min}/a_{(m-1)/2}}^{v_{\min}/a_{m/2}} dv \, dv, \quad (A.30) $$

where we remind the reader that $a_k = \cos(k\pi/N)$. For each term in this sum, the value of $m$ does not change within the range of integration, so we can apply the results which have been obtained so far. We give the full expression for the IRT, as well as the values of $k_{\pm}$ and $\cos \vartheta_{\pm}^j$ (calculated as described here) in appendix B.

## B The directionally-integrated Radon transform (IRT)

The discretised velocity distribution, with $N$ angular bins, is defined as

$$ \begin{align} & f^1(v) && \text{for } \vartheta^j \in [0, \pi/N], \nonumber \\ & \Box f^2(v) && \text{for } \vartheta^j \in [\pi/N, 2\pi/N], \nonumber \\ & \Box^k f(v) && \text{for } \vartheta^j \in [(k-1)\pi/N, k\pi/N], \nonumber \\ & \Box^{\bullet} f^N(v) && \text{for } \vartheta^j \in [(N-1)\pi/N, \pi]. \tag{B.1} \end{align} $$

The integrated Radon transform (IRT), integrated over the same angular bins, is defined as

$$ \hat{f}^j(v_{\min}) = \int_{\varphi=0}^{2\pi} \int_{a_j-1}^{a_j} \hat{f}(v_{\min}, \hat{q}) \, d\cos\vartheta \, d\varphi, \quad (B.2) $$
where have defined the shorthand $a_k = \cos(k\pi/N)$. The IRT is given in terms of the distribution functions $f^k(v)$ as

$$ \hat{f}^j(v_{\min}) = 4\pi \sum_{m=1}^{\infty} \int_{v_{\min}/a_{m/2}}^{v_{\min}/a_{(m-1)/2}} \sum_{i=0,1,2} \frac{\pi}{2} (f_{k_+}^{k_i}(v) + f_{k_-}^{k_i}(v)) (b_{i+1} - b_i) \\ + (1 - \delta_{k_+ k_-}) f_{k_+}^{k_i}(v) \mathcal{J}(b_{i+1}, a_{k_-^{i-1}}) - \mathcal{J}(b_i, a_{k_-^{i-1}}) \\ - (1 - \delta_{k_- k_+}) f_{k_-}^{k_i}(v) \mathcal{J}(b_{i+1}, a_{k_+^{i}}) - \mathcal{J}(b_i, a_{k_+^{i}}) \\ + \sum_{k=k_+^{i}+1}^{k=k_-^{i}+1} f^k(v) \mathcal{J}(b_{i+1}, a_{k-1}) - \mathcal{J}(b_i, a_{k-1}) \\ - \mathcal{J}(b_{i+1}, a_k) + \mathcal{J}(b_i, a_k) . \tag{B.3} $$

The function $J(x, y)$ depends implicitly on the speed $v$ and is defined below. For each value of $j$, we integrate over $N$ ranges in $v$, labeled by the index $m$. We also sum over three ‘bins’, labeled by the index $i = 0, 1, 2$, which are distinguished by their values of the integers $k^i$. We note that the values of $k_{\pm}$ will depend on the values of $j$ and $m$ under consideration, as well as the index $i$. We write $k_{\pm} = (k_{\pm}^0, k_{\pm}^1, k_{\pm}^2)$ and give the values below.

For odd $m$, we define $n = (m + 1)/2$ and bin edges

$$ \begin{aligned} b_0 &= \cos(j\pi/N) \\ b_1 &= \cos(\gamma + (j + n - 1)\pi/N) \\ b_2 &= \cos(\gamma + (j - n)\pi/N) \\ b_3 &= \cos((j - 1)\pi/N), \end{aligned} \tag{B.4} $$

with $\gamma = \cos^{-1}(v_{\min}/v)$. The corresponding $k_{\pm}$ values are:

$$ \mathbf{k}_{+} = \begin{cases} (n-j, n-j, n-j+1) & \text{for } j \le (n-1) \\ (1, 1, 1) & \text{for } j = n \\ (j-n+1, j-n+1, j-n) & \text{for } j \ge (n+1), \end{cases} \tag{B.5} $$

and

$$ \begin{aligned} \mathbf{k}_{-} &= (N, N, N) & \text{for } j = (N + 1 - n) \tag{B.6} \\ (2N - n - j + 1, 2N - n - j + 2, 2N - n - j + 2) & \text{for } j \ge (N + 2 - n). \end{aligned} $$

For even $m$, we write $n = m/2$ and bin edges

$$ \begin{aligned} b_0 &= \cos(j\pi/N) \\ b_1 &= \cos(\gamma + (j - n)\pi/N) \\ b_2 &= \cos(\gamma + (j + n - 1)\pi/N) \\ b_3 &= \cos((j - 1)\pi/N). \end{aligned} \tag{B.7} $$
The corresponding $k_{\pm}$ values are:

$$ 
\begin{array}{lll}
(n-j, n-j+1, n-j+1) & \text{for } j \le (n-1) & \\
\mathbf{k}_{+} = (1, 1, 1) & \text{for } j = n & \text{(B.8)} \\
(j-n+1, j-n, j-n) & \text{for } j \ge (n+1), &
\end{array}
$$

and

$$ 
\begin{array}{lll}
(n+j, n+j, n+j-1) & \text{for } j \le (N-n) & \\
\mathbf{k}_{-} = (N, N, N) & \text{for } j = (N+1-n) & \text{(B.9)} \\
(2N-n-j+1, 2N-n-j+1, 2N-n-j+2) & \text{for } j \ge (N+2-n). &
\end{array}
$$

Finally, the angular function $J$ is given by:

$$ 
\begin{aligned}
J(x, y) = & x \sin^{-1} \frac{\sqrt{y - \beta x}}{1 - x^2 \sqrt{1 - \beta^2}} \\
& + y \tan^{-1} \frac{x - \sqrt{y - \beta x}}{t} \\
& + 2 \tan^{-1} \frac{-1 \frac{1 - y^2 (y \beta^2 - \beta^4) x^2}{t} y \beta (1 + x)}{\beta} \\
& - 2 \tan^{-1} \frac{-1 \frac{1 - y^2 - \beta^2 + x \sqrt{y \beta (1-x)}}{(y+\beta) \bar{t}}}{\bar{t}},
\end{aligned}
\quad \text{(B.10)}
$$

where we define $\beta = v_{\min}/v$ and $t = (1 - x^2)(1 - \beta^2) - (y - \beta x)^2$. We note that the angular integrals in the Radon transform lead to this analytic form for arbitrary $N$, meaning that only 1-dimensional integrals over $v$ need to be performed. A python implementation of this algorithm is available on request from the author.

## References

[1] D.N. Spergel, *The Motion of the Earth and the Detection of Wimps*, Phys. Rev. D 37 (1988) 1353 [INSPIRE].

[2] N. Bozorgnia, G.B. Gelmini and P. Gondolo, *Ring-like features in directional dark matter detection*, JCAP 06 (2012) 037 [arXiv:1111.6361] [INSPIRE].

[3] C.J. Copi, J. Heo and L.M. Krauss, *Directional sensitivity, WIMP detection and the galactic halo*, Phys. Lett. B 461 (1999) 43 [hep-ph/9904499] [INSPIRE].

[4] C.J. Copi and L.M. Krauss, *Comparing WIMP interaction rate detectors with annual modulation detectors*, Phys. Rev. D 67 (2003) 103507 [astro-ph/0208010] [INSPIRE].

[5] B. Morgan, A.M. Green and N.J.C. Spooner, *Directional statistics for realistic weakly interacting massive particle direct detection experiments*, Phys. Rev. D 71 (2005) 103507 [astro-ph/0408047] [INSPIRE].

[6] J. Billard, F. Mayet, J.F. Macías-Pérez and D. Santos, *Directional detection as a strategy to discover galactic Dark Matter*, Phys. Lett. B 691 (2010) 156 [arXiv:0911.4086] [INSPIRE].

[7] A.M. Green and B. Morgan, *The median recoil direction as a WIMP directional detection signal*, Phys. Rev. D 81 (2010) 061301 [arXiv:1002.2717] [INSPIRE].

[8] P. Grothaus, M. Fairbairn and J. Monroe, *Directional Dark Matter Detection Beyond the Neutrino Bound*, Phys. Rev. D 90 (2014) 055018 [arXiv:1406.5047] [INSPIRE].
[9] J. Billard, Q. Riffard, F. Mayet and D. Santos, *Is a co-rotating Dark Disk a threat to Dark Matter Directional Detection?*, Phys. Lett. B 718 (2013) 1171 [arXiv:1207.1050] [INSPIRE].

[10] C.A.J. O'Hare and A.M. Green, *Directional detection of dark matter streams*, Phys. Rev. D 90 (2014) 123511 [arXiv:1410.2749] [INSPIRE].

[11] E. Daw, A. Dorofeev, J.R. Fox, J.L. Gauvreau, C. Ghag et al., *The DRIFT Directional Dark Matter Experiments*, EAS Publ. Ser. 53 (2012) 11 [arXiv:1110.0222] [INSPIRE].

[12] DRIFT collaboration, J.B.R. Battat et al., *First background-free limit from a directional dark matter experiment: results from a fully fiducialised DRIFT detector*, arXiv:1410.7821 [INSPIRE].

[13] Q. Riffard, J. Billard, G. Bosson, O. Bourrion, O. Guillaudin et al., *Dark Matter directional detection with MIMAC*, arXiv:1306.4173 [INSPIRE].

[14] D. Santos, G. Bosson, J.L. Bouly, O. Bourrion, C. Fourel et al., *MIMAC: Mlcro-tpc MAtrix of Chambers for dark matter directional detection*, J. Phys. Conf. Ser. 469 (2013) 012002 [arXiv:1311.0616] [INSPIRE].

[15] DMTPC collaboration, J. Monroe, *Status and Prospects of the DMTPC Directional Dark Matter Experiment*, AIP Conf. Proc. 1441 (2012) 515 [arXiv:1111.0220] [INSPIRE].

[16] DMTPC collaboration, J.B.R. Battat, *Updates from the Dark Matter Time Projection Chamber Group (DMTPC)*, J. Phys. Conf. Ser. 469 (2013) 012001 [INSPIRE].

[17] K. Miuchi, H. Nishimura, K. Hattori, N. Higashi, C. Ida et al., *First underground results with NEWAGE-0.3a direction-sensitive dark matter detector*, Phys. Lett. B 686 (2010) 11 [arXiv:1002.1794] [INSPIRE].

[18] K. Nakamura, K. Miuchi, S. Iwaki, H. Kubo, T. Mizumoto et al., *NEWAGE*, J. Phys. Conf. Ser. 375 (2012) 012013 [arXiv:1109.3099] [INSPIRE].

[19] S. Vahsen et al., *The directional dark matter detector (D³)*, EAS Publ. Ser. 53 (2012) 43.

[20] S. Burgos, J. Forbes, C. Ghag, M. Gold, V.A. Kudryavtsev et al., *Track Reconstruction and Performance of DRIFT Directional Dark Matter Detectors using Alpha Particles*, Nucl. Instrum. Meth. A 584 (2008) 114 [arXiv:0707.1758] [INSPIRE].

[21] J. Billard, F. Mayet and D. Santos, *Three-dimensional track reconstruction for directional Dark Matter detection*, JCAP 04 (2012) 006 [arXiv:1202.3372] [INSPIRE].

[22] T. Naka, M. Kimura, M. Nakamura, O. Sato, T. Nakano et al., *R&D Status of Nuclear Emulsion For Directional Dark Matter Search*, EAS Publ. Ser. 53 (2012) 51 [arXiv:1109.4485] [INSPIRE].

[23] R. Essig, J. Mardon and T. Volansky, *Direct Detection of Sub-GeV Dark Matter*, Phys. Rev. D 85 (2012) 076007 [arXiv:1108.5383] [INSPIRE].

[24] A. Drukier, K. Freese, D. Spergel, C. Cantor, G. Church and T. Sano, *New Dark Matter Detectors using DNA or RNA for Nanometer Tracking*, arXiv:1206.6809 [INSPIRE].

[25] D.R. Nygren, *Columnar recombination: a tool for nuclear recoil directional sensitivity in a xenon-based direct detection WIMP search*, J. Phys. Conf. Ser. 460 (2013) 012006 [INSPIRE].

[26] G. Mohlabeng, K. Kong, J. Li, A. Para and J. Yoo, *Dark Matter Directionality Revisited with a High Pressure Xenon Gas Detector*, arXiv:1503.03937 [INSPIRE].

[27] A.M. Green, *Astrophysical uncertainties on direct detection experiments*, Mod. Phys. Lett. A 27 (2012) 1230004 [arXiv:1112.0524] [INSPIRE].

[28] N.W. Evans, C.M. Carollo and P.T. de Zeeuw, *Triaxial haloes and particle dark matter detection*, Mon. Not. Roy. Astron. Soc. 318 (2000) 1131 [astro-ph/0008156] [INSPIRE].
[29] L.M. Widrow, Distribution functions for cuspy dark matter density profiles, Astrophys. J. Suppl. 131 (2000) 39.

[30] M. Lisanti, L.E. Strigari, J.G. Wacker and R.H. Wechsler, The Dark Matter at the End of the Galaxy, Phys. Rev. D 83 (2011) 023519 [arXiv:1010.4300] [INSPIRE].

[31] P. Bhattacharjee, S. Chaudhury, S. Kundu and S. Majumdar, Sizing-up the WIMPs of Milky Way: Deriving the velocity distribution of Galactic Dark Matter particles from the rotation curve data, Phys. Rev. D 87 (2013) 083525 [arXiv:1210.2328] [INSPIRE].

[32] M. Fornasa and A.M. Green, Self-consistent phase-space distribution function for the anisotropic dark matter halo of the Milky Way, Phys. Rev. D 89 (2014) 063531 [arXiv:1311.5477] [INSPIRE].

[33] M. Lisanti and D.N. Spergel, Dark Matter Debris Flows in the Milky Way, Phys. Dark Univ. 1 (2012) 155 [arXiv:1105.4166] [INSPIRE].

[34] M. Kuhlen, M. Lisanti and D.N. Spergel, Direct Detection of Dark Matter Debris Flows, Phys. Rev. D 86 (2012) 063505 [arXiv:1202.0007] [INSPIRE].

[35] K. Freese, P. Gondolo, H.J. Newberg and M. Lewis, The effects of the Sagittarius dwarf tidal stream on dark matter detectors, Phys. Rev. Lett. 92 (2004) 111301 [astro-ph/0310334] [INSPIRE].

[36] K. Freese, P. Gondolo and H.J. Newberg, Detectability of weakly interacting massive particles in the Sagittarius dwarf tidal stream, Phys. Rev. D 71 (2005) 043516 [astro-ph/0309279] [INSPIRE].

[37] J.I. Read, L. Mayer, A.M. Brooks, F. Governato and G. Lake, A dark matter disc in three cosmological simulations of Milky Way mass galaxies, Mon. Not. Roy. Astron. Soc. 397 (2009) 44 [arXiv:0902.0009] [INSPIRE].

[38] J.I. Read, V. Debattista, O. Agertz, L. Mayer, A.M. Brooks et al., A dark matter disc in the Milky Way, AIP Conf. Proc. 1240 (2010) 391 [arXiv:0901.2938] [INSPIRE].

[39] A. Pillepich, M. Kuhlen, J. Guedes and P. Madau, The Distribution of Dark Matter in the Milky Way's Disk, Astrophys. J. 784 (2014) 161 [arXiv:1308.1703] [INSPIRE].

[40] A.H.G. Peter, WIMP astronomy and particle physics with liquid-noble and cryogenic direct-detection experiments, Phys. Rev. D 83 (2011) 125029 [arXiv:1103.5145] [INSPIRE].

[41] M. Fairbairn, T. Douce and J. Swift, Quantifying Astrophysical Uncertainties on Dark Matter Direct Detection Results, Astropart. Phys. 47 (2013) 45 [arXiv:1206.2693] [INSPIRE].

[42] B.J. Kavanagh and A.M. Green, Improved determination of the WIMP mass from direct detection data, Phys. Rev. D 86 (2012) 065027 [arXiv:1207.2039] [INSPIRE].

[43] B.J. Kavanagh and A.M. Green, Model independent determination of the dark matter mass from direct detection experiments, Phys. Rev. Lett. 111 (2013) 031302 [arXiv:1303.6868] [INSPIRE].

[44] B.J. Kavanagh, Parametrizing the local dark matter speed distribution: a detailed analysis, Phys. Rev. D 89 (2014) 085026 [arXiv:1312.1852] [INSPIRE].

[45] P.J. Fox, J. Liu and N. Weiner, Integrating Out Astrophysical Uncertainties, Phys. Rev. D 83 (2011) 103514 [arXiv:1011.1915] [INSPIRE].

[46] P.J. Fox, G.D. Kribs and T.M.P. Tait, Interpreting Dark Matter Direct Detection Independently of the Local Velocity and Density Distribution, Phys. Rev. D 83 (2011) 034007 [arXiv:1011.1910] [INSPIRE].

[47] M.T. Frandsen, F. Kahlhoefer, C. McCabe, S. Sarkar and K. Schmidt-Hoberg, Resolving astrophysical uncertainties in dark matter direct detection, JCAP 01 (2012) 024 [arXiv:1111.0292] [INSPIRE].
[48] P. Gondolo and G.B. Gelmini, Halo independent comparison of direct dark matter detection data, JCAP 12 (2012) 015 [arXiv:1202.6359] [INSPIRE].

[49] E. Del Nobile, G.B. Gelmini, P. Gondolo and J.-H. Huh, Halo-independent analysis of direct detection data for light WIMPs, JCAP 10 (2013) 026 [arXiv:1304.6183] [INSPIRE].

[50] N. Bozorgnia, J. Herrero-Garcia, T. Schwetz and J. Zupan, Halo-independent methods for inelastic dark matter scattering, JCAP 07 (2013) 049 [arXiv:1305.3575] [INSPIRE].

[51] E. Del Nobile, G. Gelmini, P. Gondolo and J.-H. Huh, Generalized Halo Independent Comparison of Direct Dark Matter Detection Data, JCAP 10 (2013) 048 [arXiv:1306.5273] [INSPIRE].

[52] B. Feldstein and F. Kahlhoefer, A new halo-independent approach to dark matter direct detection analysis, JCAP 08 (2014) 065 [arXiv:1403.4606] [INSPIRE].

[53] B. Feldstein and F. Kahlhoefer, Quantifying (dis)agreement between direct detection experiments in a halo-independent way, JCAP 12 (2014) 052 [arXiv:1409.5446] [INSPIRE].

[54] P.J. Fox, Y. Kahn and M. McCullough, Taking Halo-Independent Dark Matter Methods Out of the Bin, JCAP 10 (2014) 076 [arXiv:1403.6830] [INSPIRE].

[55] J.F. Cherry, M.T. Frandsen and I.M. Shoemaker, Halo Independent Direct Detection of Momentum-Dependent Dark Matter, JCAP 10 (2014) 022 [arXiv:1405.1420] [INSPIRE].

[56] L.E. Strigari and R. Trotta, Reconstructing WIMP Properties in Direct Detection Experiments Including Galactic Dark Matter Distribution Uncertainties, JCAP 11 (2009) 019 [arXiv:0906.5361] [INSPIRE].

[57] J. Billard, F. Mayet and D. Santos, Markov Chain Monte Carlo analysis to constrain Dark Matter properties with directional detection, Phys. Rev. D 83 (2011) 075002 [arXiv:1012.3960] [INSPIRE].

[58] S.K. Lee and A.H.G. Peter, Probing the Local Velocity Distribution of WIMP Dark Matter with Directional Detectors, JCAP 04 (2012) 029 [arXiv:1202.5035] [INSPIRE].

[59] D.S.M. Alves, S.E. Hedri and J.G. Wacker, Dark Matter in 3D, arXiv:1204.5487 [INSPIRE].

[60] S.K. Lee, Harmonics in the Dark-Matter Sky: Directional Detection in the Fourier-Bessel Basis, JCAP 03 (2014) 047 [arXiv:1401.6179] [INSPIRE].

[61] P. Gondolo, Recoil momentum spectrum in directional dark matter detectors, Phys. Rev. D 66 (2002) 103513 [hep-ph/0209110] [INSPIRE].

[62] D.G. Cerdeño and A.M. Green, Direct detection of WIMPs, in Particle Dark Matter, Gianfranco Bertone ed., Cambridge University Press, Cambridge U.K. (2010), pp. 347–369 [arXiv:1002.1912] [INSPIRE].

[63] M. Cannoni, Reanalysis of nuclear spin matrix elements for dark matter spin-dependent scattering, Phys. Rev. D 87 (2013) 075014 [arXiv:1211.6050] [INSPIRE].

[64] D.R. Tovey, R.J. Gaitskell, P. Gondolo, Y.A. Ramachers and L. Roszkowski, A New model independent method for extracting spin dependent (cross-section) limits from dark matter searches, Phys. Lett. B 488 (2000) 17 [hep-ph/0005041] [INSPIRE].

[65] D. Albornoz Vásquez, G. Bélanger, J. Billard and F. Mayet, Probing neutralino dark matter in the MSSM and the NMSSM with directional detection, Phys. Rev. D 85 (2012) 055023 [arXiv:1201.6150] [INSPIRE].

[66] J.R. Ellis and R.A. Flores, Realistic Predictions for the Detection of Supersymmetric Dark Matter, Nucl. Phys. B 307 (1988) 883 [INSPIRE].

[67] J. Engel and P. Vogel, Spin-dependent cross-sections of weakly interacting massive particles on nuclei, Phys. Rev. D 40 (1989) 3132 [INSPIRE].

[68] F. Iachello, L.M. Krauss and G. Maino, Spin-Dependent Scattering of Weakly Interacting Massive Particles in Heavy Nuclei, Phys. Lett. B 254 (1991) 220 [INSPIRE].
[69] M.T. Ressell, M.B. Aufderheide, S.D. Bloom, K. Griest, G.J. Mathews and D.A. Resler, Nuclear shell model calculations of neutralino-nucleus cross sections for $^{29}$Si and $^{73}$Ge, Phys. Rev. D 48 (1993) 5519 [hep-ph/9307228].

[70] V. Dimitrov, J. Engel and S. Pittel, Scattering of weakly interacting massive particles from $^{73}$Ge, Phys. Rev. D 51 (1995) 291 [hep-ph/9408246] [INSPIRE].

[71] J. Menéndez, D. Gazit and A. Schwenk, Spin-dependent WIMP scattering off nuclei, Phys. Rev. D 86 (2012) 103511 [arXiv:1208.1094] [INSPIRE].

[72] P.C. Divari and J.D. Vergados, Renormalization of the Spin-dependent WIMP scattering off nuclei, arXiv:1301.1457 [INSPIRE].

[73] M. Feast and P. Whitelock, Galactic kinematics of cepheids from hipparcos proper motions, Mon. Not. Roy. Astron. Soc. 291 (1997) 683 [astro-ph/9706293] [INSPIRE].

[74] R. Schönrich, Galactic Rotation and Solar Motion from Stellar Kinematics, Mon. Not. Roy. Astron. Soc. 427 (2012) 274 [arXiv:1207.3079] [INSPIRE].

[75] J. Bovy, C.A. Prieto, T.C. Beers, D. Bizyaev, L.N. da Costa et al., The Milky Way's circular velocity curve between 4 and 14 kpc from APOGEE data, Astrophys. J. 759 (2012) 131 [arXiv:1209.0759] [INSPIRE].

[76] T. Piffl, C. Scannapieco, J. Binney, M. Steinmetz, R.-D. Scholz et al., The RAVE survey: the Galactic escape speed and the mass of the Milky Way, Astron. Astrophys. 562 (2014) A91 [arXiv:1309.4293] [INSPIRE].

[77] M. Vogelsberger, A. Helmi, V. Springel, S.D.M. White, J. Wang et al., Phase-space structure in the local dark matter distribution and its signature in direct detection experiments, Mon. Not. Roy. Astron. Soc. 395 (2009) 797 [arXiv:0812.0362] [INSPIRE].

[78] M. Kuhlen, N. Weiner, J. Diemand, P. Madau, B. Moore et al., Dark Matter Direct Detection with Non-Maxwellian Velocity Structure, JCAP 02 (2010) 030 [arXiv:0912.2358] [INSPIRE].

[79] Y.-Y. Mao, L.E. Strigari, R.H. Wechsler, H.-Y. Wu and O. Hahn, Halo-to-Halo Similarity and Scatter in the Velocity Distribution of Dark Matter, Astrophys. J. 764 (2013) 35 [arXiv:1210.2721] [INSPIRE].

[80] A. Helmi, S.D.M. White and V. Springel, The Phase space structure of a dark matter halo: Implications for dark matter direct detection experiments, Phys. Rev. D 66 (2002) 063502 [astro-ph/0201289] [INSPIRE].

[81] M. Vogelsberger, S.D.M. White, A. Helmi and V. Springel, The fine-grained phase-space structure of cold dark matter halos, Mon. Not. Roy. Astron. Soc. 385 (2008) 236 [arXiv:0711.1105] [INSPIRE].

[82] C. Savage, K. Freese and P. Gondolo, Annual Modulation of Dark Matter in the Presence of Streams, Phys. Rev. D 74 (2006) 043531 [astro-ph/0607121] [INSPIRE].

[83] M. Kuhlen, J. Diemand, P. Madau and M. Zemp, The Via Lactea INCITE Simulation: Galactic Dark Matter Substructure at High Resolution, J. Phys. Conf. Ser. 125 (2008) 012008 [arXiv:0810.3614] [INSPIRE].

[84] J. Billard, Comparing readout strategies to directly detect dark matter, Phys. Rev. D 91 (2015) 023513 [arXiv:1411.5946] [INSPIRE].

[85] E. Daw, J.R. Fox, J.L. Gauvreau, C. Ghag, L.J. Harmon et al., Spin-Dependent Limits from the DRIFT-IIId Directional Dark Matter Detector, Astropart. Phys. 35 (2012) 397 [arXiv:1010.3027] [INSPIRE].

[86] S. Burgos, E. Daw, J. Forbes, C. Ghag, M. Gold et al., First measurement of the Head-Tail directional nuclear recoil signature at energies relevant to WIMP dark matter searches, Astropart. Phys. 31 (2009) 261 [arXiv:0809.1831] [INSPIRE].