---
abstract: 'In a recent paper, a new parametrization for the dark matter (DM) speed distribution $f(v)$ was proposed for use in the analysis of data from direct detection experiments. This parametrization involves expressing the *logarithm* of the speed distribution as a polynomial in the speed $v$. We present here a more detailed analysis of the properties of this parametrization. We show that the method leads to statistically unbiased mass reconstructions and exact coverage of credible intervals. The method performs well over a wide range of DM masses, even when finite energy resolution and backgrounds are taken into account. We also show how to select the appropriate number of basis functions for the parametrization. Finally, we look at how the speed distribution itself can be reconstructed, and how the method can be used to determine if the data are consistent with some test distribution. In summary, we show that this parametrization performs consistently well over a wide range of input parameters and over large numbers of statistical ensembles and can therefore reliably be used to reconstruct both the DM mass and speed distribution from direct detection data.'
author:
- 'Bradley J. Kavanagh'
bibliography:
- 'Model\_Indep.bib'
title: 'Parametrizing the local dark matter speed distribution: a detailed analysis'
---

Introduction
============

The dark matter (DM) paradigm has enjoyed much success in explaining a wide range of astronomical observations (for a review, see e.g. Ref. [@Bertone:2005]). As yet no conclusive evidence has been provided for the identity of particle DM, though there are a variety of candidates, including the supersymmetric neutralino [@Jungman:1996], sterile neutrinos [@Dodelson:1994], axions [@Duffy:2009] and the lightest Kaluza-Klein particle [@Kolb:1984]. Here, we focus on the search for particles which belong to the generic class of Weakly Interacting Massive Particles (WIMPs). Direct detection experiments [@Goodman:1985; @Drukier:1986] aim to measure the energies of nuclear recoils induced by WIMP DM in the Galactic halo. Under standard assumptions about the DM halo, this data can be used to extract the WIMP mass and interaction cross section, allowing us to check for consistency with other search channels (such as indirect detection [@Lavalle:2012] and collider experiments [@Battaglia:2010]) and to probe underlying models of DM.

Direct detection experiments are traditionally analyzed within the framework of the Standard Halo Model (SHM), in which WIMPs are assumed to have a Maxwell-Boltzmann speed distribution in the Galactic frame. The impact of uncertainties in the WIMP speed distribution has been much studied (see e.g. Refs. [@Green:2010; @Peter:2011; @Fairbairn:2012]), leading to the conclusion that such uncertainties may introduce a bias into any reconstruction of the WIMP mass from direct detection data. As yet, the speed distribution is unknown, while a number of proposals have been put forward for its form, including analytic parametrizations (e.g. Ref. [@Lisanti:2010]) and distributions reconstructed from the potential of the Milky Way [@Bhattacharjee:2012] or from N-body simulations [@Vogelsberger:2009; @Kuhlen:2010; @Kuhlen:2012; @Mao:2012]. Recent results from N-body simulations which attempt to include the effects of baryons on structure formation also report the possible presence of a dark disk in the Milky Way [@Read:2009; @Read:2010; @Kuhlen:2013]. With such a wide range of possibilities, we should take an agnostic approach to the speed distribution, not only to avoid introducing bias into the analysis of data, but also with the hope of measuring the speed distribution and thereby probing the formation history of the Milky Way.

Several methods of evading these uncertainties have been proposed. These include simultaneously fitting the parameters of the SHM and dark matter properties [@Strigari:2009; @Peter:2009], fitting to empirical forms of the speed distribution (e.g. Ref. [@Pato:2011]) and fitting to a self-consistent distribution function [@Pato:2013]. However, these methods typically require that the speed distribution can be well fitted by a particular functional form. More model-independent methods, such as fitting the moments of the speed distribution [@Drees:2007; @Drees:2008] or using a step-function speed distribution [@Peter:2011], have also been presented. However, these methods can still introduce a bias into the measurement of the WIMP mass and perform less well with the inclusion of realistic experimental energy thresholds.

In a recent paper [@Kavanagh:2013a] (hereafter referred to as Paper 1), a new parametrization of the speed distribution was presented, which allowed the WIMP mass to be extracted from hypothetical direct detection data without prior knowledge of the speed distribution itself. Paper 1 demonstrated this for a WIMP of mass 50 GeV, using several underlying distribution functions. In the present paper, we extend this analysis to a wider range of masses. We also aim to demonstrate the statistical properties of the method and show how realistic experimental parameters affect its performance. Finally, we will also elaborate on some of the technical details of the method and assess its ability to reconstruct the underlying WIMP speed distribution.

Section \[sec:DDRate\] of this paper explains the direct detection event rate formalism and presents the parametrization of the speed distribution introduced in Paper 1. In Sec. \[sec:ParameterRecon\], the methodology for testing the parametrization is outlined. In Section \[sec:Parametrization\], we consider the choice and number of basis functions for the method. We then study the performance of the method as a function of input WIMP mass (Sec. \[sec:mass\]) and when Poisson fluctations in the data are taken into account (Sec. \[sec:stats\]). In Sec. \[sec:Recon\], we demonstrate how the speed distribution can be extracted from this parametrization and examine whether or not different distribution functions can be distinguished. Finally, we summarize in Sec. \[sec:Conclusions\] the main results of this paper.

Direct detection event rate {#sec:DDRate}
===========================

Dark matter direct detection experiments aim to measure the energies $E$ of nuclear recoils induced by interactions with WIMPs in the Galactic halo. Calculation of the event rate at such detectors has been much studied (e.g. Refs. [@Goodman:1985; @Drukier:1986; @Lewin:1996; @Jungman:1996]). For a target nucleus with nucleon number $A$, interacting with a WIMP of mass $m_\chi$, the event rate per unit detector mass is given by: $$\label{eq:Rate}
\frac{\textrm{d}R}{\textrm{d}E} = \frac{\rho_0 \sigma_p}{2 m_\chi \mu_{\chi p}^2} A^2 F^2(E) \eta(v_\textrm{min})\,,$$ where $\rho_0$ is the local dark matter mass density, $\sigma_p$ is the WIMP-proton spin-independent cross section and the reduced mass is defined as $\mu_{A B} = m_A m_B/(m_A + m_B)$. The Helm form factor $F^2(E)$ [@Helm:1956] describes the loss of coherence of spin-independent scattering due to the finite size of the nucleus. A wide range of possible interactions have been considered in the literature, including inelastic [@Smith:2001], isospin-violating [@Kurylov:2003] and more general non-relativistic interactions [@Fan:2010; @Fitzpatrick:2012; @Fitzpatrick:2013]. We focus here on the impact of the WIMP speed distribution on the direct detection event rate. We therefore restrict ourselves to considering only spin-independent scattering, which is expected to dominate over the spin-dependent contribution for heavy nuclei, due to the $A^2$ enhancement in the rate.

Information about the WIMP velocity distribution $f(\textbf{v})$ is encoded in the function $\eta$, sometimes referred to as the mean inverse speed, $$\label{eq:eta}
\eta(v_\textrm{min}) = \int_{v > v_\textrm{min}} \frac{f(\textbf{v})}{v} \, \textrm{d}^3\textbf{v}\,,$$ where $\textbf{v}$ is the WIMP velocity in the reference frame of the detector. The integration is performed only over those WIMPs with sufficient speed to induce a nuclear recoil of energy $E$. The minimum required speed for a target nucleus of mass $m_N$ is $$\label{eq:v_min}
v_\textrm{min}(E) = \sqrt{\frac{m_N E}{2\mu_{\chi N}^2}}\,.$$

We distinguish between the directionally averaged velocity distribution $$f(v) = \oint f(\textbf{v}) \, \textrm{d}\Omega_{\textbf{v}}\,,$$ and the 1-dimensional speed distribution $$f_1(v) = \oint f(\textbf{v}) v^2 \textrm{d}\Omega_{\textbf{v}}\,.$$ The distribution function should in principle be time-dependent, due to the motion of the Earth around the Sun. However, this is expected to be a percent-level effect (for a review, see e.g. Ref. [@Freese:2013]) and we therefore assume that $f_1(v)$ is time independent in the present work.

We consider several benchmark speed distributions in this work, including the SHM and the SHM with the addition of a moderate dark disk which accounts for 23% of the total WIMP density [@Kuhlen:2013]. We model the speed distributions as combinations of Gaussian functions in the Earth frame $$\label{eq:gaussian}
g(\textbf{v}) = N \exp\left(-\frac{(\textbf{v} - \textbf{v}_\textrm{lag})^2}{2\sigma_v^2}\right) \Theta(v_\textrm{esc} - |\textbf{v} - \textbf{v}_\textrm{lag}|)\,,$$ where $\textbf{v}_\textrm{lag}$ specifies the peak velocity of the distribution in the Earth frame and $\sigma_v$ the velocity dispersion. We truncate the distribution above the escape speed $v_\textrm{esc}$ in the Galactic frame and the factor $N$ is required to satisfy the normalization condition (Eq. \[eq:normalization\]). We use the value $v_\textrm{esc} = 544 {\,\textrm{km s}^{-1}}$, which lies within the 90% confidence limits obtained from the RAVE survey [@RAVE:2007; @RAVE:2013]. In addition, we also use the speed distribution of Lisanti et al. [@Lisanti:2010], which has the following form in the Earth’s frame: $$\label{eq:lisanti}
f(\textbf{v}) = N \left[\exp\left(\frac{v_\textrm{esc}^2 - |\textbf{v} - \textbf{v}_0|^2}{k v_0^2}\right) -1\right]^k \Theta(v_\textrm{esc} - |\textbf{v} - \textbf{v}_0|)\,.$$ We use the parameter values $k = 2$ and $v_0 = 220 {\,\textrm{km s}^{-1}}$ in this work. We summarize in Tab. \[tab:distributions\] the different speed distributions considered. We also plot several of these in Fig. \[fig:Ensemble\_distributions\] for reference.

[m[3cm]{}|ccc]{} Speed distribution benchmark & Fraction & $v_\textrm{lag} / {\,\textrm{km s}^{-1}}$ & $\sigma_v / {\,\textrm{km s}^{-1}}$\
SHM & 1 & 220 & 156\
& 0.77 & 220 & 156\
& 0.23 & 50 & 50\
Stream & 1 & 400 & 20\
& 0.97 & 220 & 156\
& 0.03 & 500 & 20\
& 0.5 & 200 & 20\
& 0.5 & 400 & 20\
Lisanti et al. & & $v_0 = 220 {\,\textrm{km s}^{-1}}$ & $k = 2$

![Several of the benchmark speed distributions used in this work. They are defined in Eqs. \[eq:gaussian\] and \[eq:lisanti\] with parameters from Tab. \[tab:distributions\]. These distributions are the SHM (solid blue), SHM+DD (dashed green), Lisanti et al. (dot-dashed red) and the stream (dotted magenta).[]{data-label="fig:Ensemble_distributions"}](SpeedDistributions-Ensemble.pdf){width="49.00000%"}

In Paper 1, a parametrization for the WIMP speed distribution was introduced, for use in the analysis of direct detection data. The parametrization of Paper 1 has the form: $$\label{eq:parametrization}
f_1(v) = v^2 \exp\left\{ -\sum_{k =0}^{N-1} a_k P_k(v)\right\}\,$$ where $P_k(v)$ is some basis of polynomial functions of v. We fit the coefficients $\left\{a_1, ..., a_{N-1}\right\}$ using data, and fix $a_o$ by normalization $$\label{eq:normalization}
a_0 = \ln\left(\int_{0}^\infty v^2 \exp\left\{ -\sum_{k = 1}^{N-1} a_k P_k(v)\right\} \, \textrm{d}v\right)\,.$$ This form of parametrization ensures that the distribution function $f_1(v)$ is everywhere positive and can be used to fit an arbitrary underlying directionally-averaged distribution function (given a sufficiently large number of polynomial basis functions). We explore in Sec. \[sec:Parametrization\] which basis functions should be used in the parametrization, as well as how many basis functions are required.

Parameter Reconstruction {#sec:ParameterRecon}
========================

In order to assess the performance of the parametrization method, we attempt to reconstruct the WIMP mass $m_\chi$ and polynomial coefficients $\left\{a_1, ..., a_{N-1} \right\}$ using the nested sampling software <span style="font-variant:small-caps;">MultiNest</span> [@MultiNest1; @MultiNest2; @MultiNest3]. We also include the WIMP-proton spin-independent cross section $\sigma_p$ as a free parameter. However, we are forced to treat the cross section as a nuisance parameter. As has previously been noted [@Kavanagh:2012; @Kavanagh:2013a], taking an agnostic approach to the DM speed distribution means that we do not know what fraction of WIMPs lie above the energy thresholds of the experiments. While this does not adversely impact the reconstruction of the WIMP mass, it does result in a strong degeneracy, such that only lower limits can be placed on the cross section using such methods. In any case, the cross section appears in the event rate (Eq. \[eq:Rate\]) in the degenerate combination $\rho_0 \sigma_p$. Uncertainties on the local DM density $\rho_0$ are at present on the order of a factor of 2 (see e.g. [@Iocco:2011; @Bovy:2012; @Zhang:2013; @Nesti:2013]) and thus any reconstruction of the cross section would be subject to the same systematic uncertainty. In this work, we focus instead on reconstructing the WIMP mass and the shape of the speed distribution. For concreteness, we use the values $\sigma_p = 10^{-45} \textrm{ cm}^2$ and $\rho_0 = 0.3 \textrm{ GeV cm}^{-3}$ throughout this work.

Experimental benchmarks {#sec:experiments}
-----------------------

In order to generate mock data sets, we consider three idealized mock experiments, loosely based on detectors which are currently in development. As previous work has shown [@Kavanagh:2012; @Peter:2013a], the WIMP mass and speed distribution are degenerate when data from only a single experiment is considered. However, this degeneracy can be broken by including data from additional experiments with different nuclear target masses. The three target materials we consider here are Xenon, Argon and Germanium. We describe each experiment in terms of its nucleon number $A$, fiducial detector mass $m_\textrm{det}$, efficiency $\epsilon$ and energy sensitivity window $\left[E_\textrm{min}, E_\textrm{max}\right]$. We incorporate the effects of detector sensitivity, analysis cuts and detector down-time into the value of the efficiency $\epsilon$, which we take to be energy independent for simplicity. We consider a total exposure time for all experiments of $t_\textrm{exp} = \textrm{ 2 years}$. The experimental parameter values used in this work are summarized in Tab. \[tab:experiments\].

[c|m[1.2cm]{}m[1.7cm]{}m[1.5cm]{}m[1.7cm]{}]{} Experiment & Target Mass, $A$ & Detector Mass (fid.), $m_\textrm{det}$/kg & Efficiency, $\epsilon$ & Energy Range/keV\
Xenon & 131 & 1100 [@Aprile:2012a] & 0.7 [@Aprile:2012b] & 7-45 [@Aprile:2010]\
Argon & 40 & 1000 & 0.9 [@Benetti:2007] & 30-100 [@Grandi:2005]\
Germanium & 73 & 150 [@Bauer:2013b] & 0.6 [@Bauer:2013a] & 8-100 [@Bauer:2013a]\

The exact parameter values we used in this work do not strongly impact the results we present. However, it is important to note that the total mass and exposure of the experiments will affect the total number of events observed. This in turn will affect the precision of the reconstructions. For example, we have chosen a total Argon mass of 1000 kg. This is the stated target for Argon-based experiments which are in development (e.g. Ref. [@Badertscher:2013]), though at present typical fiducial masses for Argon prototypes are of the order of 100 kg [@Grandi:2005]. The data we have generated does not represent the ‘high-statistics’ regime: across all three experiments the total number of events observed is roughly 200-300 with as few as 10 events in the Germanium detector for some scenarios. Using a smaller exposure (or equivalently a smaller interaction cross section) will reduce the precision of the results, but should not introduce any additional bias. We also briefly consider the impact of a *larger* number of events in Sec. \[sec:Recon\].

Parameter sampling
------------------

We make parameter inferences using a combination of Bayesian and frequentist statistics. Bayes theorem for the probability of a particular set of theoretical parameters $\boldsymbol{\Theta}$ given the observed data $\textbf{D}$ is:

$$P(\boldsymbol{\Theta}|\textbf{D}) = \frac{P(\boldsymbol{\Theta}) P(\textbf{D}|\boldsymbol{\Theta})}{P(\textbf{D})}\,,$$

where $P(\boldsymbol{\Theta})$ is the prior on the parameters and $P(\textbf{D})$ is the Bayesian evidence, which acts as a normalizing factor and has no impact on parameter inference. We summarize the priors used in this work in Tab. \[tab:priors\]. We also summarize in Tab. \[tab:MultiNest\] the MultiNest sampling parameters used.

[m[1in]{}|cc]{} Parameter & Prior type & Prior range\
$m_\chi / \textrm{ GeV}$ & log-flat & $\left[10^{0}, 10^{3}\right]$\
$\sigma_p / \textrm{ cm}^2$ & log-flat & $\left[10^{-46}, 10^{-42}\right]$\
$\left\{a_k\right\}$ & linear-flat & $\left[-50, 50\right]$\
$R_{BG} / \textrm{dru}$ & log-flat & $\left[10^{-12}, 10^{-5}\right]$\

       Parameter         Value
  ------------------- -----------
   $N_\textrm{live}$     10000
      efficiency         0.25
       tolerance       $10^{-4}$

  : Summary of the MultiNest sampling parameters used in this work.[]{data-label="tab:MultiNest"}

The factor $P(\textbf{D}|\boldsymbol{\Theta})$ is simply the likelihood of the data given the parameters $\boldsymbol{\Theta}$. In Sec. \[sec:Parametrization\] and Sec. \[sec:mass\], we consider the effects of varying the form of the parametrization and of varying the input WIMP mass. In order to eliminate the effects of Poisson noise, we use Asimov data [@Cowan:2013] for these sections. This means that we divide the energy window of each experiment into bins of width 1 keV. We then set the observed number of events $N_{o,i}$ in bin $i$ equal to the expected number of events $N_{e,i}$. In this case, we use the binned likelihood, calculated for $N_b$ energy bins: $$\label{eq:binnedL}
\mathcal{L}_b = \prod_{i = 1}^{N_b} \frac{N_{e,i}^{N_{o,i}} \textrm{e}^{-N_{e,i}}}{N_{o,i}!}\,.$$

In Sec. \[sec:stats\] and Sec. \[sec:Recon\], we consider many realisations of data, including the effects of Poisson noise. We therefore use the extended likelihood which has previously been used by both the Xenon [@Aprile:2011] and CDMS [@Ahmed:2009] collaborations, which for a single experiment is given by: $$\label{eq:unbinnedL}
\mathcal{L} = \frac{N_e^{N_o} \textrm{e}^{-N_e}}{N_o!} \prod_{i = 1}^{N_o} P(E_i)\,,$$ where the expected number of events is given by: $$\label{eq:N_expected}
N_e = \epsilon m_\textrm{det}t_\textrm{exp}\int_{E_\textrm{min}}^{E_\textrm{max}} \frac{\textrm{d}R}{\textrm{d}E}\, \textrm{d}E\,,$$ and the normalised recoil spectrum is given by: $$\label{eq:eventdistribution}
P(E) = \frac{ \epsilon m_\textrm{det} t_\textrm{exp}}{N_e} \frac{\textrm{d}R}{\textrm{d}E}\,.$$ The total likelihood is then the product over all experiments considered.

Using nested sampling, we can extract the full posterior probability distribution of the parameters $P(\boldsymbol{\Theta}|\textbf{D})$, as well as the likelihood $\mathcal{L}(\Theta)$. However, we often want to make inferences not jointly for all parameters but for only a subset (treating the remaining as nuisance parameters). If we conceptually partition the parameter space into the parameters of interest $\boldsymbol{\psi}$ and the remaining nuisance parameters $\boldsymbol{\phi}$, we would like to make inferences about the values of $\boldsymbol{\psi}$, without reference to the values of $\boldsymbol{\phi}$. One option for doing this is to calculate the marginalized posterior distribution, obtained by integrating the posterior probability over the parameters we are not interested in:

$$P_m(\boldsymbol{\psi}) = \int P(\boldsymbol{\psi}, \boldsymbol{\phi}) \, \textrm{d}\boldsymbol{\phi}\,.$$

This method performs well for small numbers of observations (compared to the number of free parameters in the fit). We take the mode of the distribution to be the reconstructed parameter value and construct p% *minimal* credible intervals, which include those parameter values with $P_m(\boldsymbol{\psi}) \geq h$, where $h$ is chosen such that p% of the probability distribution lies within the interval. The marginalized posterior method is used in Sec. \[sec:stats\] and Sec. \[sec:Recon\], where in some cases the number of events observed in an experiment is less than 10.

An alternative method is to calculate the profile likelihood. This is obtained by maximizing the full likelihood function over the nuisance parameters: $$\label{eq:profilelikelihood}
\mathcal{L}_p(\boldsymbol{\psi}) = \max_{\boldsymbol{\phi}} \mathcal{L}(\boldsymbol{\psi},\boldsymbol{\phi})\,.$$ For a large number of observations, we can take the value which maximizes $\mathcal{L}_p$ as the reconstructed value and construct confidence intervals using the asymptotic properties of the profile likelihood. We use the profile likelihood for parameter inferences in Sec. \[sec:Parametrization\] and Sec. \[sec:mass\], as the Asimov data sets provide a large number of measurements of $N_{e,i}$ over a large number of bins. The profile likelihood can also lead to less noisy reconstructions than the marginalized posterior, especially when the dimensionality of the parameter space becomes high, as in Sec. \[sec:Parametrization\] and Sec. \[sec:mass\].

Testing the parametrization {#sec:Parametrization}
===========================

We now consider the two questions:

how many basis functions are required and

which polynomial basis should be used?

In order to answer these questions, we use the two benchmark distribution functions illustrated in Fig. \[fig:VaryingN\_distributions\]. We have chosen these benchmarks not because they are necessarily realistic distribution functions but because they should be difficult to fit using standard techniques and fitting functions (e.g. [@Lisanti:2010]). The first distribution (referred to as ‘bump’) is a SHM distribution with the addition of a small bump, which contributes just 3% of the total WIMP population and could correspond to a small sub-halo or stream [@Vogelsberger:2009]. This should be difficult to fit because it represents only a very small deviation from the standard scenario. The second distribution (referred to as ‘double-peak’) has a sharp and rapidly varying structure, which we anticipate should be difficult to capture using a small number of basis functions.

![Benchmark speed distributions used in Sec. \[sec:Parametrization\] to test the performance of the parametrization as a function of the number and type of basis functions.[]{data-label="fig:VaryingN_distributions"}](SpeedDistributions-VaryingN.pdf){width="49.00000%"}

Varying the number of basis functions
-------------------------------------

We first investigate how the reconstructed WIMP mass $m_\textrm{rec}$ and uncertainty varies with the number of basis functions $N$. For now, we fix our choice of basis to shifted Legendre polynomials, as used in Paper 1:

$$P_k(v) = L_k\left(2\frac{v}{v_\textrm{max}} - 1\right)\,,$$

where $L_k$ is the Legendre polynomial of order $k$, and $v_\textrm{max}$ is a cut off for the parametrization. We should choose $v_\textrm{max}$ to ensure that $f_1(v)$ is negligible above the cut off. However, too high a choice of $v_\textrm{max}$ will result in $f_1(v)$ being close to zero over a large range of the parametrization, making fitting more difficult. We use the value $v_\textrm{max} = 1000 {\,\textrm{km s}^{-1}}$, while lies significantly above the Galactic escape speed.

The lower panel of Fig. \[fig:BUMP\_LEG\] shows the best fit mass and 68% confidence intervals as a function of $N$, using as input a WIMP of mass 50 GeV and the ‘bump’ distribution function. The reconstructed mass very rapidly settles close to the true value, using as few as three basis functions. This is because adding the bump near $v \sim 500 {\,\textrm{km s}^{-1}}$ still leaves the mean inverse speed relatively smooth, so a large number of basis function are not required. The correct mass is reconstructed and we emphasize in the lower panel of Fig. \[fig:BUMP\_LEG\] that the reconstruction is stable with the addition of more basis functions.

We should also consider how the quality of the fit changes as a function of $N$. We would expect that adding fit parameters should always lead to a better fit. Eventually, the fit should be good enough that adding additional basis functions will no longer improve it significantly. We can then be confident that our reconstruction is accurate and not an artifact of using too few basis functions. In order to investigate this, we utilise the Bayesian Information Criterion (BIC) [@Schwarz:1978], which is given by:

$$BIC = 2N_p\textrm{ln}(N_m) - \textrm{ln}(\mathcal{L}_\textrm{max}) \, ,$$

where $N_p$ is the number of free parameters, $N_m$ is the number of measurements or observations and $\mathcal{L}_\textrm{max}$ is the maximum likelihood value obtained in the reconstruction. For the case of binned data, $N_m$ corresponds simply to the total number of energy bins across all experiments. This criterion penalises the inclusion of additional free parameters and in comparing several models, we should prefer the one which minimises the BIC.

![Bayesian information criterion (BIC) as a function of the number of basis functions for an underlying ‘bump’ distribution function, 50 GeV WIMP and using Legendre polynomial basis functions (upper panel). Also shown (lower panel) are the reconstructed WIMP mass (dashed blue line), 68% confidence interval (shaded blue region) and underlying WIMP mass (solid horizontal black line).[]{data-label="fig:BUMP_LEG"}](VaryingN_BUMP_LEG.pdf){width="49.00000%"}

The upper panel of Fig. \[fig:BUMP\_LEG\] shows the BIC (in arbitrary units) as a function of the number of basis functions for the ‘bump’ distribution function. The BIC is comparable for the cases of $N=2$ and $N=3$, indicating that the quality of the fit is improved slightly by the addition of another basis function. However, adding further basis functions does not have a significant impact on the maximum likelihood, leading to an increase in the BIC. This coincides with the stabilization of the reconstructed mass around the true value and we conclude that only two or three basis functions are required to provide a good fit to the data.

Figure \[fig:DP\_LEG\] shows the corresponding results for the ‘double-peak’ distribution function. Here, we note that the bias induced by using too small a number of basis functions is larger than for the case of the ‘bump’ distribution, due to the more complicated structure in this case. The BIC is minimized for $N=7$, indicating that additional basis functions do not significantly improve the quality of the fit to data. This suggests that the shape of the speed distribution can be well fit by $N\geq7$ basis functions. As shown in the lower panel of Fig. \[fig:DP\_LEG\], the reconstruction of the WIMP mass is stable around the true mass for these values of $N$.

![As Fig. \[fig:BUMP\_LEG\] but for an underlying ‘double-peak’ distribution function.[]{data-label="fig:DP_LEG"}](VaryingN_DP_LEG.pdf){width="49.00000%"}

We propose that such a procedure should be used in the case of real data should a dark matter signal be observed at multiple detectors. We have shown that by analyzing the reconstructed mass as a function of $N$ we can recover the true mass and that by using the BIC we can be confident that we have obtained an adequate fit to data.

Choice of basis functions
-------------------------

We now consider the second question posed at the start of Sec. \[sec:Parametrization\]: which polynomial basis should be used? We see immediately that a naive power series of the form

$$\textrm{ln}f(v) \approx a_0 + a_1 v + a_2 v^2 + a_3 v^3 + ...\,,$$

is not practical for the purposes of parameter estimation. Higher powers of $v$ will have rapidly growing contributions to $\textrm{ln} f$, meaning that the associated coefficients must be rapidly decreasing in order to suppress these contributions. Fitting to the SHM using just 5 terms, the range of values for the $a_k$ in the case of a simple power series would span around 13 orders of magnitude. Ideally, we would like to specify an identical prior on each of the coefficients. However, in this scenario this would result in a highly inefficient exploration of the parameter space when some of the terms are so small.

This problem can be significantly improved by rescaling $v$. We choose to rescale by a factor of $v_\textrm{max} = 1000 {\,\textrm{km s}^{-1}}$, and cut off the distribution function at $v_\textrm{max}$. The basis functions $(v/v_\textrm{max})^k$ are now less than unity by construction and the coefficients $a_k$ are now dimensionless:

$$\textrm{ln}f(v) \approx a_0 + a_1 (v/v_\textrm{max}) + a_2 (v/v_\textrm{max})^2 + a_3 (v/v_\textrm{max})^3 + ...\,.$$

We now address the problem of *conditioning* of the polynomial basis (see e.g. Refs. [@Gautschi:1978; @Wilkinson:1984]). Conditioning is a measure of how much the value of a polynomial changes, given a small change in the coefficients. For a well-conditioned polynomial, small changes in the coefficient are expected to lead to small changes in the value of the polynomial. This is ideal for parameter estimation as it leads to a more efficient exploration of the parameter space. Orthogonal polynomial basis functions typically have improved conditioning [@Gautschi:1978] and we consider two specific choices: the Legendre polynomials which have already been considered and the Chebyshev polynomials. The Chebyshev polynomials are used extensively in polynomial approximation theory [@Mason:2002] and are expected to be well conditioned [@Gautschi:1978].

We have checked that the reconstruction results using Chebyshev polynomials are largely indistinguishable from the case of Legendre polynomials for both the ‘bump’ and ‘double-peak’ distributions and as a function of $N$. This leads us to conclude that the accuracy of the reconstruction is independent of the specific choice of basis. However, the reconstruction was much faster in the case of the Chebyshev basis. This is illustrated in Fig. \[fig:times\], which shows the time taken for reconstruction of the ‘bump’ benchmark as a function of $N$. The time taken grows much more slowly for the Chebyshev basis (roughly as $N^2$) than for the Legendre basis (roughly as $N^3$). We have also checked that this difference is not an artifact of how we calculate the basis functions. These results indicate that this choice of basis provides both reliable and efficient reconstruction for the WIMP mass and we therefore use the Chebyshev basis in the remainder of this work.

![Time taken (using 4 processors in parallel) for the reconstruction of the ‘bump’ benchmark, as a function of number of basis functions. The time taken using the Chebyshev basis (blue squares) grows more slowly with $N$ than for the Legendre basis (red triangles).[]{data-label="fig:times"}](RunTimes.pdf){width="49.00000%"}

Varying $m_\chi$ {#sec:mass}
================

In previous work [@Kavanagh:2013a], this parametrization method was only tested for a single WIMP mass of $50 \textrm{ GeV}$. Here, we extend this analysis to a wider range of WIMP masses. We generate Asimov data for WIMP masses of 10, 20, 30, 40, 50, 75, 100, 200 and 500 GeV and reconstruct the best fit WIMP mass $m_\textrm{rec}$ and 68% and 95% confidence intervals from the profile likelihood. We use the SHM as a benchmark distribution function and use a fixed number of $N=5$ basis functions. The results are shown in Fig. \[fig:VaryingM\], along with the line $m_\textrm{rec} = m_\chi$ for reference.

For large values of $m_\chi$, the shape of the event spectrum becomes independent of $m_\chi$ [@Green:2008], which results in a widening of the confidence intervals as the WIMP mass increases. For low mass WIMPs, fewer events are observed in each bin, again resulting in wider confidence intervals. It should be noted that for this analysis we have used Asimov data, in which the exact (non-integer) number of events is recorded in each bin. For low mass WIMPs, this means that the spectrum (and therefore the correct WIMP mass) is still well reconstructed using Asimov data, in spite of the small number of events. The tightest constraints are obtained when the input WIMP mass is close to the masses of several of the detector nuclei (in the range 30-80 GeV). There also appears to be no bias in the WIMP mass: the reconstruction matches the true mass across all values considered.

![Reconstructed WIMP mass $m_\textrm{rec}$ (central dashed blue line) as a function of input WIMP mass $m_\chi$ as well as 68% and 95% intervals (inner and outer blue dashed lines respectively). The line $m_\textrm{rec} = m_\chi$ (solid red line) is also plotted for reference.[]{data-label="fig:VaryingM"}](VaryingM.pdf){width="49.00000%"}

An alternative parametrization method was proposed in Ref. [@Kavanagh:2012], in which the *momentum* distribution of halo WIMPs was parametrized. For a given speed distribution, the corresponding momentum distribution may be broad and easily reconstructed for high mass WIMPs. However, for low mass WIMPs the momentum distribution would be much narrower, owing to their lower momenta. The momentum parametrization method therefore performs poorly for low mass WIMPs. The parametrization presented in this paper does not suffer from similar problems.

So far, we have only considered idealized direct detection experiments. We now apply the method to more realistic mock detectors, taking into account the effects of finite energy resolution, as well as unrejected background events. We assume here that each experiment has a gaussian energy resolution with fixed width $\sigma_E = 1 \textrm{ keV}$, such that the observed event rate for recoils of energy $E$ is given by:

$$\frac{\textrm{d}R}{\textrm{d}E} = \int_{0}^{\infty} \frac{1}{\sqrt{2 \pi} \sigma_E}\exp\left\{-\frac{(E-E')^2}{2\sigma_E^2}\right\} \frac{\textrm{d}{R'}}{\textrm{d}E'} \, \textrm{d}E'\,,$$

where the primed event rate is the underlying (perfect resolution) rate. We also assume a constant flat background rate for each experiment $R_\textrm{BG} = 10^{-6}$ events/kg/keV/day (which has been suggested as a possible background rate for Xenon1T [@Aprile:2010] and WArP-100L [@Grandi:2005]) when generating mock data sets. However, we allow the flat background rate in each experiment to vary as free parameters during the fit.

We have chosen relatively generic resolution and background parameters in this work, because the precise details of energy resolution and background shape and rate will depend on the specific experiment under consideration. Instead, we hope to show that the inclusion of more realistic experimental setups does not introduce an additional bias or otherwise spoil the good properties of the method presented here.

Figure \[fig:VaryingM\_real\] shows the reconstructed mass as a function of input mass in this more realistic scenario. The 68% and 95% confidence intervals are now wider and the reconstructed mass does not appear to be as accurate. For input masses above $\sim$100 GeV, the uncertainties become very wide, with only a lower limit of $m_\textrm{rec} > 20 \textrm{ GeV}$ being placed on the WIMP mass. Due to the poorer energy resolution the shape of the energy spectrum is less well-determined. In addition, a flat background contribution can mimic a higher mass WIMP, as it leads to a flatter spectrum. This leads to a strong degeneracy, as a wide range of mass values can provide a good fit to the data. For high input masses, the profile likelihood is approximately constant above $m_\textrm{rec} \sim 20 \textrm{ GeV}$, indicating that there is no sensitivity to the underlying WIMP mass.

In spite of this, the true mass values still lie within the 68% and 95% confidence intervals. In addition, the poor values for the reconstructed mass for heavy WIMPs are a side effect of the loss of sensitivity. Because the profile likelihood is approximately flat, the maximum likelihood point is equally likely to be anywhere within the 68% interval. These effects would be present even if we had considered a fixed form for the speed distribution. However, when we allow for a range of possible speed distributions, the effects become more pronounced. These results show that for more realistic experimental scenarios, the method presented in this paper remains reliable over a range of masses, though its precision may be significantly reduced.

![As fig. \[fig:VaryingM\] but including the effects of finite energy resolution and non-zero backgrounds, as described in the text.[]{data-label="fig:VaryingM_real"}](VaryingM_real.pdf){width="49.00000%"}

Statistical properties {#sec:stats}
======================

We now consider the impact of statistical fluctuations on the reconstruction of the WIMP mass. In reality, the number of events observed $N_o$ at a given experiment will be Poisson distributed about the expected value $N_e$, while the observed distribution of recoil energies will not exactly match that expected from the calculated event rate. The fundamental statistical limitations of future direct detection experiments have been studied in detail in Ref. [@Strege:2012]. In this work, we generate 250 realisations of data from the mock experiments described in Tab. \[tab:experiments\]. Each realisation of the mock data is generated as follows:

1.  Calculate the number of expected events $N_e$, given $\left\{m_\chi, \sigma_p, f(v)\right\}$, using Eq. \[eq:N\_expected\],

2.  Pick the number of observed events $N_o$ from a Poisson distribution with mean $N_e$,

3.  Pick recoil energies $\left\{E_1, E_2, ..., E_{N_o}\right\}$, from the distribution $P(E)$ in Eq. \[eq:eventdistribution\],

4.  Repeat for all three experiments.

For each realisation, we then use the method described in Sec. \[sec:ParameterRecon\] (using $N = 5$ basis functions) to reconstruct the WIMP mass and 68% and 95% credible intervals. Figure \[fig:Realisations\] shows the distribution of reconstructed masses for an input mass of 50 GeV for three benchmark speed distributions: SHM, SHM+DD and Lisanti et al. as described in Sec. \[sec:DDRate\]. In all three cases, the reconstructions are peaked close to the true value, regardless of the underlying distribution. For the SHM+DD distribution, the spread of reconstructions is slightly wider (with more reconstructions extending up to higher masses). This is due to the smaller number of events for this benchmark, making the data sets more susceptible to Poisson fluctuations.

In order to assess the accuracy of the reconstructed value of the mass $m_\textrm{rec}$, we also calculate the bias $b$ for each realisation:

$$\label{eq:bias}
b = \textrm{ln}(m_\textrm{rec} / \textrm{GeV}) - \textrm{ln}(m_\textrm{true} / \textrm{GeV})\,.$$

We compare the logarithms of the mass values because we have used logarithmically-flat priors on the WIMP mass. In Tab. \[tab:bias\] we show the average bias across all 250 realisations for each of the three benchmark distributions. In all three cases, the average bias is consistent with zero. Even in the SHM+DD case, which shows larger fluctuations away from the true value, there is no statistical bias.

![Distribution of the reconstructed mass $m_\textrm{rec}$ for 250 mock data sets generated using several benchmark speed distributions, defined in Sec. \[sec:DDRate\]. These are the SHM (top), SHM+DD (middle) and Lisanti et al. (bottom) distributions. The input WIMP mass of $m_\chi = 50 \textrm{ GeV}$ is shown as a vertical dashed red line.[]{data-label="fig:Realisations"}](SHM_ensemble.pdf "fig:"){width="49.00000%"} ![Distribution of the reconstructed mass $m_\textrm{rec}$ for 250 mock data sets generated using several benchmark speed distributions, defined in Sec. \[sec:DDRate\]. These are the SHM (top), SHM+DD (middle) and Lisanti et al. (bottom) distributions. The input WIMP mass of $m_\chi = 50 \textrm{ GeV}$ is shown as a vertical dashed red line.[]{data-label="fig:Realisations"}](DD_ensemble.pdf "fig:"){width="49.00000%"} ![Distribution of the reconstructed mass $m_\textrm{rec}$ for 250 mock data sets generated using several benchmark speed distributions, defined in Sec. \[sec:DDRate\]. These are the SHM (top), SHM+DD (middle) and Lisanti et al. (bottom) distributions. The input WIMP mass of $m_\chi = 50 \textrm{ GeV}$ is shown as a vertical dashed red line.[]{data-label="fig:Realisations"}](LIS_ensemble.pdf "fig:"){width="49.00000%"}

[m[1in]{}|c]{} Benchmark speed distribution & Mean bias $\langle b \rangle$\
SHM & 0.002 $\pm$ 0.008\
SHM+DD & 0.005 $\pm$ 0.007\
Lisanti et al. & 0.01 $\pm$ 0.01\

We also test the *coverage* of the credible intervals which have been constructed. For a $p\%$ credible interval, we expect that the true parameter value of the WIMP mass will lie within the interval in $p\%$ of realisations. In this case, we say that the method provides *exact coverage*. However, if the true parameter lies within the interval in fewer than $p\%$ of realisations, our reconstructed credible intervals are too narrow and provide *undercoverage*. Alternatively, we obtain *overcoverage* when the true parameter lies within the interval more often that $p\%$ of the time. Table \[tab:coverage\] shows the coverage values for the $68\%$ and $95\%$ intervals obtained in this section. In each case, there is very close to exact coverage. We have also checked that these intervals only provide exact coverage for the true WIMP mass of 50 GeV. Other values of $m_\textrm{rec}$ are contained within the intervals less frequently than the true value, again indicating that this parametrization allows for unbiased and statistically robust reconstructions of the WIMP mass.

[m[1in]{}|cc]{} Benchmark speed distribution & 68% coverage & 95% coverage\
SHM & 71 $\pm$ 3 % & 94 $\pm$ 3 %\
SHM+DD & 68 $\pm$ 3 % & 91 $\pm$ 4 %\
Lisanti et al. & 70 $\pm$ 3 % & 95 $\pm$ 3 %\

Reconstructing $f_1(v)$ {#sec:Recon}
=======================

Using the method described in this paper, we can obtain the posterior probability distribution for the coefficients $\left\{ a_1, ..., a_{N-1}\right\}$ given the data, which we refer to as $P(\textbf{a})$. We would like to be able to present this information in terms of the distribution function $f_1(v)$ in order to compare with some known distribution or look for particular features in the distribution. However, due to the fact that the distribution function is normalized, the values of $f_1$ at different speeds will be strongly correlated. We illustrate here how robust comparisons with benchmark distributions can be made.

As a first step, we can attempt to sample from the $P(\textbf{a})$, in order to obtain $P(f_1(v))$. This is the probability distribution for the value of $f_1$ at a particular speed $v$, marginalizing over the values of $f_1$ at all other speeds. We can repeat for a range of speeds to obtain 68% and 95% credible intervals for the whole of $f_1(v)$. The result of this procedure is presented in Fig. \[fig:f\], for a randomly selected realisation from the SHM ensemble of Sec. \[sec:stats\]. The underlying SHM distribution is shown as a solid line, while the 68% and 95% marginalized intervals are shown as dark and light shaded regions respectively. In this naive approach, we see that there is little shape information which can be recovered from the reconstruction, with only upper limits being placed on the speed distribution.

![Reconstructed speed distribution for a single realisation of data, generated for a 50 GeV WIMP. 68% and 95% credible intervals are shown as dark and light shaded regions respectively, while the underlying SHM distribution function is shown as a solid blue line.[]{data-label="fig:f"}](f_SHM.pdf){width="49.00000%"}

This method performs poorly because, as initially mentioned in Sec. \[sec:ParameterRecon\], we have no information about the fraction of dark matter particles below the energy threshold of our experiments. If this fraction is large, the event rate for a given cross-section is suppressed. However, increasing the cross-section will increase the total event rate. There is thus a degeneracy between the shape of the speed distribution and the cross-section, meaning that we can only probe the shape of $f_1(v)$, rather than its overall normalization. This degeneracy has not been accounted for in Fig. \[fig:f\]. We can attempt to correct for this by adjusting the normalization of $f_1(v)$. If we fix $f_1(v)$ to be normalized to unity above $v_a$ (where $v_a \approx 171 {\,\textrm{km s}^{-1}}$ is the lowest speed probed by the experiments for a WIMP of mass 50 GeV), we can compare the shapes of the underlying and reconstructed distribution functions. This is illustrated in Fig. \[fig:f\_scaled\], which shows that we now broadly reconstruct the correct shape of $f_1(v)$. Below $v_a$, the value of $f_1(v)$ is poorly constrained, because the experiments provide no information about the shape of the distribution below theshold.

There remain several issues with this approach. In order to utilize this method, we must know the approximate value of the lowest speed probed by the experiments. However, this value is set by the WIMP mass. We could determine $v_a$ using the reconstructed WIMP mass, but this would be subject to significant uncertainty. In addition, direct reconstructions of the speed distribution are easily biased. The upper limit of the energy windows of the experiments corresponds to a particular WIMP speed (for a given WIMP mass). WIMPs above this speed still contribute to the total event rate, but contribute no spectral information. The reconstructed shape of the high speed tail of the distribution is therefore not constrained by the data, but may affect the reconstructed value of $f_1$ at lower speeds.

![Reconstructed speed distribution for the same realisation of data as Fig. \[fig:f\]. In this case, we have also normalized $f_1(v)$ to unity above $v_a \approx 171 {\,\textrm{km s}^{-1}}$ (vertical dashed line). This is the lowest speed accessible to the experiments for a WIMP of mass 50 GeV. 68% and 95% credible intervals are shown as dark and light shaded regions respectively, while the underlying SHM distribution function is shown as a solid blue line.[]{data-label="fig:f_scaled"}](f_SHM_scaled_line.pdf){width="49.00000%"}

An alternative approach is to reconstruct the mean inverse speed $\eta(v)$ (defined in Eq. \[eq:eta\]) at some speed $v$. Because $\eta(v)$ is an integral function of $f_1$, it is less prone to bias as it takes into account the full shape of the distribution at speeds greater than $v$. However, we do not know the normalization of $f_1$ and so we must normalize $\eta$ appropriately. For each point sampled from $P(\textbf{a})$, we calculate $\eta$. We then divide by $\alpha(v)$, the fraction of WIMPs above speed $v$, calculated using the same parameter point: $$\label{eq:alpha}
\alpha(v) = \int_{v}^{\infty} f_1(v') \, \textrm{d}v'\,.$$

We will write this rescaled mean inverse speed as $\eta^*(v) = \eta(v)/\alpha(v)$. The value of $\eta^*(v)$ is a measure of the shape of the distribution function above $v$. However, information about the normalization of the distribution has been factored out by dividing by $\alpha(v)$. We no longer need to know the value of $v_a$ in order to obtain information about the shape of the distribution at higher speeds. We may still need to decide the speed down to which we trust our reconstruction, but this no longer relies on an arbitrary choice of $v_a$ to normalize the reconstructions at all speeds.

In Fig. \[fig:eta\_stats\], we plot the mean reconstructed value of $\eta^*$ at several values of $v$, using 250 realisations of the 50 GeV SHM benchmark. We also show the mean upper and lower limits of the 68% credible intervals as errorbars. The form of $\eta^*$ for the SHM is shown as a solid blue line. In all cases except for $v=100 {\,\textrm{km s}^{-1}}$, the mean reconstructed value is close to the true value, indicating that $\eta^*$ can be reconstructed without bias using this method. At low speeds, the reconstructed value deviates from the true value. In addition, the credible intervals lead to *under*coverage in the $v=100 {\,\textrm{km s}^{-1}}$ case. However, this point lies below the lowest speed to which the experiments are sensitive and therefore we cannot trust the reconstruction at this low speed. We have checked that for the remaining values of $v$ the method provides exact or overcoverage, indicating that at higher speeds we can use $\eta^*$ as a reliable and statistically robust measure of the shape of the distribution.

![Mean reconstructed values of the rescaled mean inverse speed $\eta(v)/\alpha(v)$ at several values of $v$, calculated over 250 realisations of data using a 50 GeV WIMP and underlying SHM distribution function. Errorbars indicate the mean upper and lower limits of the 68% credible intervals. The underlying form of $\eta(v)/\alpha(v)$ obtained from the SHM is shown as a solid blue line.[]{data-label="fig:eta_stats"}](Eta.pdf){width="49.00000%"}

In the case of a single realisation of data, we would like to compare the probability distribution for $\eta^*(v)$ (obtained from $P(\textbf{a})$) to the value calculated from some test distribution. We note that several distributions may produce the same value of $\eta^*(v)$ at a given value of $v$. Thus, we may fail to reject a distribution function which is not the true distribution. However, if the calculated value of $\eta^*(v)$ does lie outside the $p\%$ interval, we can reject it at the $p\%$ level.

We can increase the discriminating power of this method by repeating this reconstruction over all speeds and checking to see if the benchmark value of $\eta^*$ is rejected at any value of $v$. The result of this procedure is shown in Fig. \[fig:eta\] for a single realisation of data generated using an SHM distribution (the same data as in Figs. \[fig:f\] and \[fig:f\_scaled\]). We plot the 68%, 95% and 99% credible intervals as shaded regions, as well as the values of $\eta^*(v)$ calculated from several benchmark speed distribution. We will focus on the intermediate speed range ($v \gtrsim 200 {\,\textrm{km s}^{-1}}$), as we do not know *a priori* the lowest speed to which the experiments are sensitive.

![Rescaled mean inverse speed $\eta(v)/\alpha(v)$, reconstructed from a single realisation of data using a 50 GeV WIMP and underlying SHM distribution function. At each value of $v$ we calculate 68%, 95% and 99% credible intervals (shown as shaded intervals). We also show the calculated values of $\eta(v)/\alpha(v)$ for several possible benchmark speed distributions: SHM (solid blue), SHM+DD (dashed green), Lisanti et al. (dot-dashed red) and stream (dotted magenta). The benchmark curves are truncated when the underlying distribution function goes to zero.[]{data-label="fig:eta"}](SHM_lores.pdf){width="49.00000%"}

The reconstructed intervals are consistent with a range of possible distribution functions. The SHM and SHM+DD distributions are identical over a wide range of speeds. This is because above $\sim 200 {\,\textrm{km s}^{-1}}$, the two distributions differ in normalization but not in shape. Differences appear between the two at low speeds where their shapes diverge. The Lisanti et al. distribution results in a larger deviation from the SHM, but not sufficiently large to differentiate between the two distributions given the size of the uncertainties. Finally, the stream distribution results in a significantly different form for $\eta^*(v)$. At approximately $400 {\,\textrm{km s}^{-1}}$, the curve for the stream distribution lies outside the reconstructed 99% credible interval. We can therefore use this method to reject the stream distribution at the 99% confidence level.

Figure \[fig:eta\_hires\] shows the results of a reconstruction using a larger exposure. In this case, we generate data using the Lisanti et al. distribution and an exposure increased by a factor of $2.5$, resulting in approximately 1000 events across the three detectors. As expected, the resulting credible intervals are now substantially narrower. The stream distribution now lies significantly outside the 99% interval. In Fig. \[fig:eta\_hires\_zoom\], we show the same results, but focusing in on the region around $v \sim 400 {\,\textrm{km s}^{-1}}$. At certain points, the SHM and SHM+DD distributions now lie outside the 95% credible interval, suggesting that with a number of events of the order of 1000, we may be able to reject these benchmarks.

![As Fig. \[fig:eta\], but using as input a Lisanti et al. speed distribution and an exposure time which is 2.5 times longer.[]{data-label="fig:eta_hires"}](LIS_hires.pdf){width="49.00000%"}

![As Fig. \[fig:eta\_hires\], but focusing on the region around $v \sim 400 {\,\textrm{km s}^{-1}}$. Notice that in the range $400-550 {\,\textrm{km s}^{-1}}$, both the SHM and SHM+DD curves lie at or below the lower limit of the 95% credible interval.[]{data-label="fig:eta_hires_zoom"}](LIS_hires_zoom.pdf){width="49.00000%"}

While the method displayed in Fig. \[fig:f\_scaled\] allows the approximate shape of the speed distribution to be reconstructed, reconstructions of $\eta^*(v)$ allow more statistically robust statements to be made about the underlying speed distribution. In particular, Fig. \[fig:eta\_hires\_zoom\] illustrates that with larger exposures deviations from Maxwellian speed distributions can be detected in a model-independent fashion.

Conclusions {#sec:Conclusions}
===========

We have studied in detail the parametrization for the local dark matter speed distribution introduced in Paper 1. This method involves writing the logarithm of the speed distribution as a polynomial in speed $v$ and fitting the polynomial coefficients (along with the WIMP mass and cross section) to the data. We have attempted to disentangle in this paper the influence of different benchmark speed distributions, different benchmark WIMP masses and different forms for the parametrization. We summarize our conclusions as follows:

-   We have shown that the reconstruction of the WIMP mass is robust under changes in the number of basis functions $N$. We have used the Bayesian Information Criterion (BIC) to compare models with different values of $N$ and have shown that minimizing the BIC allows us to determine how many basis functions are required for a reliable reconstruction. We have also demonstrated that the results of the method do not depend strongly on the choice of basis functions, but that the speed of reconstructions may improved by using the Chebyshev polynomial basis.

-   We have shown that the method leads to unbiased reconstructions of the WIMP mass for masses in the range 10-500 GeV. Including realistic experimental parameters, including non-zero backgrounds and finite energy resolution, reduces the precision of these reconstructions. In particular, for large values of the input mass, we can only place a lower limit of approximately 20 GeV on the reconstructed mass. This is significantly lower than in the idealized case, where we can typically constrain the WIMP mass to be heavier than around 50 GeV.

-   We have used several ensembles of data realisations to demonstrate the statistical properties of the method, including unbiased reconstructions and exact coverage of the WIMP mass.

-   We have presented several ways of displaying the reconstructed WIMP speed distribution using this method. In order to make robust statistical inferences about the speed distribution, we calculate the probability distribution of $\eta(v)/\alpha(v)$. This is the mean inverse speed $\eta(v)$, which appears in the direct detection event rate (eq. \[eq:Rate\]), rescaled by the fraction of WIMPs $\alpha(v)$ above speed $v$. This can be used as a measure of the *shape* of the distribution function, from which the unknown normalization has been factored out. We can then compare to the expected value of $\eta(v)/\alpha(v)$ from a given benchmark speed distribution, allowing us to distinguish between different underlying models.

We have shown that this parametrization method is statistically robust and works well over a large range of input parameters, both in terms of particle physics and astrophysics. The inclusion of more realistic experimental parameters does not introduce any additional bias, but does reduce the precision of reconstructions. We obtain unbiased estimates of the WIMP mass over large numbers of data sets. Finally, we have shown that we can distinguish different forms of the speed distribution. With around 1000 events, it may be possible to detect minor deviations from the Standard Halo Model and begin to search for more interesting structure in the speed distribution of the Milky Way.

The author thanks Anne M. Green and Mattia Fornasa for helpful comments. BJK is supported by STFC. Access to the University of Nottingham High Performance Computing Facility is also gratefully acknowledged.
