Discretising the velocity distribution for directional dark matter experiments

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.


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. The resulting recoil spectrum will be peaked in the opposing direction, towards the constellation of Cygnus (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][4][5][6][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] and even DNA [24].
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 [25]. However, such a simple description is unlikely to be an accurate description of the Milky Way halo [26][27][28][29][30]. In addition, N -body simulations raise the possibility of more complicated distributions, including debris flows [31,32], tidal streams [33,34] or a dark disk [35][36][37]. 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 χ and interaction cross section σ) 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 [38][39][40][41][42].
Astrophysics-independent approaches to non-directional experiments have received significant attention in the past few years, including so-called 'halo-independent' methods [43][44][45][46][47][48][49][50][51][52][53] and attempts at parametrising the 1-dimensional speed distribution [38,40,41,54]. There have also been several attempts to construct a suitable parametrisation for the 3-dimensional velocity distribution, relevant for directional experiments [55][56][57][58]. 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. [57,58]) 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 Sec. 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. [38,40]). Instead, we investigate the magnitude of the error in the event rate introduced by this angular discretisation for two different benchmark velocity distributions. We leave the choice of a suitable parametrisation of the remaining 1-dimensional functions of v to future work. 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 Sec. 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 Sec. 3 by a discussion of previous attempts to overcome these uncertainties. In Sec. 4, we present the discretised velocity distribution which is the focus of this work. In Sec. 5 and 6, 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 Sec. 7. 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.

Directional formalism
The directional rate per unit detector mass for recoils of energy E R can be written as [59] dR whereq = (sin θ cos φ, sin θ sin φ, cos φ) 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 σ 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 DMproton rate. The local DM density is written as ρ 0 , the DM mass as m χ and the reduced DM-proton mass as µ χp = m χ m p /(m χ + m p ). 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 spindependent (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 [60,61] where J is the total nuclear spin, S p,n 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) [62] or values motivated by specific models (e.g. a n /a p = ±1) [63]. 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 = √ 2m N E R , for a nucleus of mass m N : 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 [61]. The spin structure functions (and the proton and neutron spin matrix elements) can be calculated using nuclear shell models (see e.g. Refs. [64][65][66][67][68][69]). Here, we focus on Fluorine targets, assuming a p = a n , with spin structure functions and matrix elements taken from Ref. [70]. 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 transformf where v min is the minimum DM speed required to excite a nuclear recoil of energy E R : Geometrically, the Radon transform corresponds to an integral over f (v) on a plane perpendicular toq, which has a perpendicular distance v 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 directionq. 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 ρ(r) ∝ r −2 , the resulting velocity distribution has a Maxwell-Boltzmann form in the referenceframe of the Galaxy, [59] with velocity dispersion σ v . The corresponding Radon transform also takes the form of a Gaussian,f Using the properties of the Radon transform we can obtain the corresponding result in the Earth frame by performing a Galilean transformation v → v + v e , with v e the velocity of the Earth with respect to the rest frame of the halo, in which casê The SHM velocity distribution in the laboratory frame is then (2.10) Within the SHM, one typically assumes values of v e ≈ 220 km s −1 and σ v = v e / √ 2 ≈ 156 km s −1 . However, there are still uncertainties on these values of at least 10% [71][72][73]. Introducing a cut off in the distribution at the Galactic escape speed v esc introduces further uncertainty into the model [74].
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 [26] or with more realistic density profiles [27] have been suggested, as well as analytic parametrisations which should provide more realistic behaviour at low and high speeds [28]. Self-consistent distribution functions reconstructed from the potential of the Milky Way have also been obtained [29,30].
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 [75][76][77]. 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 [32].
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 [35,36]. The result is a dark disk which corotates with the stellar disk, though with a smaller velocity dispersion, typically σ DD v ∼ 50 km s −1 . Such a dark disk could comprise a large fraction (up to 50%) of the total DM density, though more recent results [37] 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 [78,79]. However, local structures could still be significant, especially streams of DM from the tidal disruption of Milky Way satellites, such as Sagittarius [33,34,80].
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.

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 [55]. The possibility of including additional structures, such as streams and dark disks, has also been considered [56]. 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 [57] 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 [81]. 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 [58]. 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 χ or σ 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 : 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.

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: Over each bin in θ , 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 θ . We note that in this work we will only be considering the azimuthally-averaged event rate (i.e. integrated over the recoil angle φ). In this case, we emphasise that no assumptions about the φ -dependence of f (v) are required, as the integral over φ 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 φ without loss of generality. However, this analysis could be extended to consider bins in the φ angle, with an additional discretisation in φ 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 [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 should 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 Fig. 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 φ . These approximations are obtained by averaging the full velocity distribution over each bin in θ : 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 Sec. 5. Upon discretisation, the distribution which is initially focused in one direction (towards θ = 0 • ) now becomes constant over θ in each bin. For the SHM (left column), there is a predominantly forward-going component (cos θ > 0), as well as a smaller backwards component (cos θ < 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 (θ = 90 • ), 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): whereq = (sin θ cos φ, sin θ sin φ, cos φ). 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 off 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 θ or cos θ). However, for the reader interested in simply calculatingf 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 eachf 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: where β = v min /v. 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 andq. We then bin these events in angle and compare with the calculated distribution off j (v min ). The results are in perfect agreement with the results expected from Appendix B.

Comparison with exact results: Radon transform
We now wish to compare these approximate IRTs with the IRTs obtained from the full (nondiscretised) velocity distribution. To do this, we select a benchmark velocity distribution (such as the SHM) and calculate the f k (v) of Eq. 4.1 by averaging over cos θ in each bin, as in Eq. 4.2 and Fig. 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 θ.
We fix the peak of the underlying speed distribution to be aligned in the forward direction θ = 0 • . However, we discuss the consequences for 'misaligned' distributions in Sec. 7. We consider two example distributions, which have already been illustrated in Fig. 1. The first is the canonical Standard Halo Model (SHM), with parameters v e = 220 km s −1 and σ 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 k (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 σ 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 [33,34,80]. This stream should therefore be considered a 'worst-case' scenario which is difficult to approximate accurately. Figure 2 shows this 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). 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 (j = 1) Radon transform, the error induced by the discretisation is less than 20%, while for the backward (j = 2) case, the error is typically larger (50 − 100%), though the absolute value is smaller.

N=2 discretisation
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 Fig. 1. The full stream distribution (top right) is strongly focused in the forward (θ = 0 • ) direction, meaning that particles must scatter through almost 90 • to produce a recoil is the backward direction. This is kinematically allowed only for a small fraction of the population. Due to the angular averaging, however, the N = 2 velocity distribution (middle right) has a significant population of DM particles with velocities at right angles (θ = 90 • ) 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 compares the approximate and exact IRTs for the N = 3 discretisation illustrated on the bottom row of Fig. 1. In the case of the SHM, the agreement between the two functions is improved compared to the N = 2 case. In the forward direction (left panel), the fractional error is reduced to at most 10%. For the stream, there is a slight improvement in the approximate result. In particular, scattering in the backwards direction (right panel) is substantially reduced compared to the N = 2 case. Again, the reason for this is clear from Fig. 1: the population of DM particles with transverse velocities (θ = 90 • ) is reduced compared to the N = 2 case due to the finer discretisation. 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 −1 to 500 km s −1 . This is because almost all particles are moving in the forward direction with speed v = 500 km s −1 . The kinematics of the scattering requires that v ·q = v min , meaning that scattering in the forward direction (θ = 0 • ) is only allowed for v min ≈ 500 km s −1 . Scattering away from the forward direction, but still included in the first angular bin (θ < 60 • ), is kinematically allowed for v min > cos(60 • ) × 500 km s −1 ≈ 250 km s −1 . For the discretised velocity distribution, there is a population of DM particles initially directed away from the forward direction (θ = 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.

N=3 discretisation
Finally, looking at 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 Fig. 1), it is still kinematically possible for particles on the edge of the forward bin (θ ≈ 60 • ) to scatter into the backward direction (θ > 120 • ). Thus, the approximate calculation predicts a small but non-zero IRT for the j = 3 case.

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 Fig. 4. 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 Fig. 1 that most of the distribution is distributed over an angle θ 10 • , 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.
Finally, 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 k = 1 forward-directed bin. 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.

Comparison with exact results: event rate
In the previous section, we examined the discretisation error in the integrated Radon transform (IRT) which arises from the use of the discretised velocity distribution. However, direc-  tional direct detection experiments do not directly measure the Radon Transform, but rather the differential recoil rate dR/dE R dΩ q of Eq. 2.1. Analogously to the previous sections, we define the directionally-integrated recoil spectrum: The integrated recoil spectrum is suppressed by the form factor and is therefore related to the IRT by where v min (E R ) is defined in Eq. 2.6. In analogy to Sec. 5, we can now compare the directionally-integrated spectrum obtained using the discretised velocity distribution (approximate) and using the full Radon transform (exact). Figure 5 shows an example of this comparison, for the case of N = 2 discretisation. This is similar to Fig. 2, but now showing the integrated recoil spectrum, rather than the integrated Radon transform. 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  Figure 5. Exact (solid) and approximate (dashed) forward and backward differential event rates, dR 1 /dE R and dR 2 /dE R , defined in Eq. 6.1, for a Fluorine target, DM mass m χ = 50 GeV and for N = 2. Results are shown for the SHM (blue) and stream (red) distribution functions. Shown as a shaded band are energies below 20 keV, a typical directional detection threshold [83]. DM mass of m χ = 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. Energies below 20 keV -a typical threshold energy for directional detection experiments [83] -are shown by a grey shaded band.

N=2 discretisation
Compared to Fig. 2, the results of Fig. 5 show a slight suppression of the rate at high energies, relative to the rate at low energies, due to the effects of the form factor. We note that the range of energies below threshold is also where the discrepancy between the approximate and exact results is largest. This is relatively generic; as was discussed in Sec. 5, the discretisation typically leads to a large discrepancy in scattering events at low v min , corresponding to large-angle scattering. The result is that the experiment is only sensitive to recoil events in the range of energies where the agreement between the approximate and exact rates is closest.
We note also from Fig. 5 that, as in the case of Fig. 2, the approximate rate and exact rate are much closer in the case of the SHM than in the case of the stream distribution. However, we note that the two distribution functions result in substantially different event spectra (in terms of both shape and normalisation). The result is that the discretisation error between the exact and approximate rates for the stream is typically much smaller than the difference between the SHM and stream distributions to begin with. This demonstrates why such an approximate, discretised approach is justified. The error in the event rate induced by assuming the incorrect distribution function is much larger than the possible error induced by the discretisation.

N=3 discretisation and the folded spectrum
We briefly give one final example of the discretised event rate. Specifically, we choose the N = 3 case, analogous to the results of Fig. 3. This 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 [84], 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 θ and another with direction − cos θ. They therefore probe the so-called 'folded' recoil spectrum: For the N = 2 case, the folded spectrum is entirely isotropic, as the forward and backward IRTs differ only in the sign of cos θ. However, in the N = 3 case, the transverse IRT, given byf 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 Fig. 3) are folded together, to obtain the longitudinal Radon transform Analogously, we can define the longitudinal and transverse event spectra, dR L,T /dE R , derived from the corresponding Radon transforms. In Fig. 6, we show the approximate and exact recoil spectra for the N = 3 discretisation, while in Fig. 7 we show the corresponding longitudinal and transverse (folded) spectra, relevant when sense discrimination is not possible. The largest difference between the two sets of spectra is that the longitudinal rate (left panel of Fig. 7) shows slightly better agreement between the approximate and exact spectra in the case of the SHM (compared to the left and right panels of Fig. 6). In addition, the largest discrepancy now appears at very low recoil energy, below the energy threshold of the experiment. The discretisation procedure leads to an angular smearing of the distribution function, leading to more scattering events in the backwards direction and fewer in the forward direction, as previously discussed. The folding of the event rates causes this effect to average out, improving the agreement between the exact and approximate rates. In the case of the stream, this 'leakage' of the events into the backwards direction 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 addition, that it can in fact reduce the associated discretisation error.

Forward-backward asymmetry
In order to better quantify how large the error induced by discretisation could be, we can compare the number of events obtained using the exact and approximate IRTs. Because the N = 1 case (full average over all angles) is exact, the discretised velocity distribution will always recover the correct total number of scattering events. We have checked this explicitly up to N = 10 by summing the number of events in each angular bin. Instead of comparing the total number of events, then, we can instead compare the forward-backward asymmetry in the number of events, A FB . In terms of the number of forward and backward events between the threshold and maximum energies E th and E max ,  Figure 6. Exact (solid) and approximate (dashed) directionally-integrated differential event rates, dR j /dE R , defined in Eq. 6.1, for a Fluorine target, DM mass m χ = 50 GeV and for N = 3. Results are shown for the SHM (blue) and stream (red) distribution functions. Shown as a shaded band are energies below 20 keV, a typical directional detection threshold [83].  . Exact (solid) and approximate (dashed) longitudinal and transverse differential event rates, dR L,T /dE R , defined in Eqs. 6.4 and 6.5 and the associated text, for a Fluorine target, DM mass m χ = 50 GeV and for N = 3. Results are shown for the SHM (blue) and stream (red) distribution functions. Shown as a shaded band are energies below 20 keV, a typical directional detection threshold [83].
the asymmetry is given by By comparing the values of A 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 8 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). We again consider a Fluorine target with E th = 20 keV, E max = 100 keV and m χ = 50 GeV. In the case of the stream (red), the event rate is strongly asymmetric (A FB ≈ 1) as the velocity distribution is highly focused in the forward direction. By comparison, the SHM (blue) has a smaller asymmetry (A FB ≈ 0.9) due to its wider velocity dispersion.
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 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 forwardbackward asymmetry. Even with only 4 angular bins, the error in A FB for both velocity distributions is less than 10%. This is smaller than the difference in A FB between the two distribution functions, indicating that the discretised velocity distribution can potentially be useful in discriminating between possible distributions, even for a relatively small number of bins.

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).
This means that we should be able to parametrise the speed distributions in each angular bin f k (v) and account for astrophysical uncertainties without introducing a large error in the recoil spectrum. For comparison, we have also considered a more extreme and highlydirectional stream distribution. The disagreement between the true recoil spectrum and that obtained from the discretised velocity distribution is larger in this case, though increasing the number of bins N improves the agreement. However, there is a significant difference between the event spectra produced by the SHM and stream distributions, which is typically larger than the error induced by the discretisation of f (v).
In Sec. 6, we have attempted to analyse how closely the discretised velocity distribution can reproduce the recoil spectrum which is actually measured in directional detectors. To this end, we examined both the directionally-integrated recoil spectra and the forward-backward asymmetry in the event rate. However, translating from the Radon transform to the event rate is not necessarily straightforward in a realistic experiment. Conversion from v min to E R requires knowledge of the DM mass m χ and attempting to fit both the velocity distribution and mass simultaneously can lead to spurious results [38,40]. Furthermore, for comparison with experimental data, we must take into account the fact that experiments have finite angular resolution, typically in the range 20 • -80 • , 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 underlying velocity distribution (i.e. v e aligned along θ = 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 Sec. 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 θ = 0 • would also lead to a distribution more isotropic in θ . 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.
For a realistic analysis, all of these issues should be taken into account. Based on mock (or future) data, the optimal direction for θ = 0 • 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. 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.
Once a suitable parametrisation for the radial functions f k (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 (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 (v). We have therefore focused here on the discretisation of f (v) and left more involved investigations for future work.

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: In Appendices A and B, we have provided recipes for calculating the corresponding directionallyintegrated Radon transforms (IRTs), which appear in the directional event rate. Alternative methods for parametrising the DM velocity distribution (such as Refs. [57,58], 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. For the standardly-assumed SHM, the IRT is well recovered even for a small number of bins N . The typical fractional error in the Radon transform is around 20% for N = 2, falling to 10% for N = 3 and decreasing further with increasing N . For comparison, we also consider an extreme, highly directional stream distribution. In this case, the discretisation error is much larger, with an estimated N = 18 angular bins required to accurately map the full velocity distribution. However, the difference in event spectra between the SHM and stream distributions is substantial and typically much larger than the error caused by using a discretised velocity distribution. This further emphasises the importance of developing a suitable parametrisation for f (v) to mitigate uncertainties in the directional rate.
In order to better quantify the error induced by this discretisation, 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 ≥ 4 bins. 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 ≥ 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. [38,40,42]), 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.

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 We write the coordinates of the velocity and recoil momentum as v = v sin θ cos φ , sin θ sin φ , cos θ q = (sin θ cos φ, sin θ sin φ, cos θ) , (A. 2) and consider here only the azimuthally integrated Radon transform. We write this aŝ f (v min , cos θ), which is given bŷ We expand the δ−function explicitly in terms of the angular coordinates: We rewrite the argument of the δ−function as a function φ: Here, we sum over those values of φ i satisfying g(φ i ) = 0: where we have also defined β = v min /v. The solutions for φ ∈ [0, 2π] are: We note that these solutions exist only for β ∈ [0, 1] (or equivalently v > v min ) and for α ∈ [−1, 1], otherwise Eq. A.6 cannot be satisfied. If these constraints are satisfied, there exist exactly 2 solutions for a given value of φ and therefore 2 δ-functions in Eq. A.5.
For the derivative of g(φ) we obtain Substituting the values of φ 1,2,3,4 , we see that Each of the two δ-functions therefore contributes the same amount to the integral, regardless of the value of φ . We can now perform the integral over φ: where C(α) = 1 for α ∈ [−1, 1] and vanishes otherwise. We note that D is independent of the azimuthal angle φ . The constraint α ∈ [−1, 1] is satisfied for cos θ ∈ [x − , x + ], where We therefore obtain f (v min , cos θ) = 2 x + x − because I does not depend on φ . We note from this that the azimuthally-integrated Radon transform is unaffected by the exact φ -dependence of f (v). Instead, the azimuthallyintegrated 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 φ without loss of generality. We will consider a velocity distribution which is discretised into N angular pieces: (A.14) We would then ultimately like to calculate the directionally-integrated Radon transformf j , integrated over the same bins in cos θ as for the velocity distribution, Here, we have introduced the notation a k = cos(kπ/N ) and defined f (v, cos θ ) = 2πf (v, cos θ , φ ). We divide the integration range for cos θ into different regions: . . .
where k ± are defined such that: These regions are chosen such that f (v, cos θ ) is independent of cos θ within each region. We note that if k + = k − , there is only one region: cos θ ∈ [x − , x + ]. Using these definitions we can therefore rewrite the term in brackets in Eq. A.15 as x + x − f (v, cos θ )I(β, cos θ, cos θ ) d cos θ x + a k + I(β, cos θ, cos θ ) d cos θ . (A.18) We can now explicitly perform the integral over cos θ , where we have defined We note that F (x ± , cos θ) = ± π 2 , (A. 21) for all cos θ. We also note that F (y, cos θ) can be integrated analytically: x F (y, cos θ) cos θ where it is understood that J(x, y) is a function of β and where (A.23) It is now possible to also complete the integration over cos θ. Assuming that k ± do not change over a range cos θ ∈ [b i , b i+1 ], we can combine Eqs. A.18 and A.22 to obtain However, we must take into account the fact that the values of x ± and therefore of k ± depend on β and cos θ. We now discuss how to divide up the integration regions of v and cos θ such that we can apply Eq. A.24.
As an illustration, we show in Fig. 9 the values of x + (x − ) as solid (dashed) lines as a function of cos θ, for a fixed value of β. 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 θ. As an example, we show with horizontal and vertical dotted lines the edges of the discretised regions in cos θ and cos θ for N = 3. That is, the dotted lines show the values cos(nπ/N ) for n = 0, 1, 2, 3. If we wish to evaluatef 1 (v min ), we need to integrate over the shaded region in Fig. 9.
In Fig. 10, we show x ± for different values of β. We give two examples, for two different ranges of β: β < 1/2 (black) and β > √ 3/2 (blue). The values of k ± (and therefore the relevant f k (v)) will clearly depend on the value of β, so we must be careful to account for this properly.
Calculation of the values of k ± and the corresponding bin edges in cos θ ({b i }) is involved, though not technically difficult. We therefore sketch the procedure here. We begin with the definition of x + :  Figure 10. As Fig. 9, but for two different values of β. The black lines show x ± for a value of β in the range β < 1/2 while blue lines show x ± for a value in the range β > √ 3/2.
By straight-forward differentiation, we observe that this function has a maximum at cos θ = β.
For concreteness, we consider β in the range, We must consider three separate cases, depending on j (i.e. which cos θ bin we are considering): (i) j ≤ (n − 1) In this case, cos θ > β, meaning that ∂x + /∂ cos θ < 0, so x + is monotonically decreasing in cos θ.
(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 θ bin by using the definition in Eq. A.25 . This allows us to determine the value of k + (i.e. the cos θ 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 θ bin edges exactly once, when cos θ = cos θ j + , satisfying x + (cos θ j + ) = cos(lπ/N ) , (A. 27) for some l = 1, 2, ..., N . This means that for a given j, there are two distinct regions in cos θ (cos θ < cos θ j + and cos θ > cos θ 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 θ 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 θ ∈ [b i , b i+1 ]. It is clear then that we must subdivide each cos θ bin into 3 smaller bins. The edges of these sub-bins will be cos θ j + and cos θ j − , at which the value of either k + or k − changes. Finally, we need to determine which of these crossings cos θ j ± occurs first. The crossings occur at the same value of cos θ if cos θ j + = cos θ j − , which occurs for β = cos ((n − 1/2)π/N ). For each value of n, we therefore need to distinguish two regimes, depending on which of the crossings cos θ j ± is larger. It is helpful then to further subdivide the range of β, writing We can now determine which of cos θ 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 θ ({b i }) are now well determined, along with the values of k ± 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, where we remind the reader that a k = cos(kπ/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 ± and cos θ j ± (calculated as described here) in Appendix B.

(B.9)
Finally, the angular function J is given by: (B.10) where we define β = v min /v and t = (1 − x 2 )(1 − β 2 ) − (y − β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.