A new halo-independent approach to dark matter direct detection analysis

Uncertainty in the local dark matter velocity distribution is a key difficulty in the analysis of data from direct detection experiments. Here we propose a new approach for dealing with this uncertainty, which does not involve any assumptions about the structure of the dark matter halo. Given a dark matter model, our method yields the velocity distribution which best describes a set of direct detection data as a finite sum of streams with optimised speeds and densities. The method is conceptually simple and numerically very efficient. We give an explicit example in which the method is applied to determining the ratio of proton to neutron couplings of dark matter from a hypothetical set of future data.

Uncertainty in the local dark matter velocity distribution is a key difficulty in the analysis of data from direct detection experiments. Here we propose a completely new approach for dealing with this uncertainty, which does not involve any assumptions about the structure of the dark matter halo. Given a dark matter model, our method yields the velocity distribution which best describes a set of direct detection data as a finite sum of streams with optimized speeds and densities. The method is conceptually simple and numerically very efficient. We give an explicit example in which the method is applied to determining the ratio of proton to neutron couplings of dark matter from a hypothetical set of future data. Introduction. The direct search for dark matter (DM) in shielded underground detectors is a promising strategy not only for confirming the particle nature of DM, but also for measuring its key properties, such as mass and couplings to nucleons. A central problem in the analysis of such experiments, however, is the uncertainty in the DM velocity distribution f (v) [1][2][3]. Numerical simulations indicate that assuming an isotropic and isothermal halo may not be a good approximation [4] and that in addition there may be localized streams of DM [5] as well as a DM disk co-rotating with the stars [6].
Such uncertainties are particularly important for light DM, which only probes the tail of f (v) [7,8], but they also significantly reduce the amount of information that can be inferred about general DM candidates. In fact, since changes in the halo structure can mimic changes in the DM parameters, a single direct detection experiment is insufficient to determine even the DM mass. To make progress, one has to find a way to combine information from different target materials and quantify the impact of astrophysical uncertainties [9][10][11].
The standard approach to this problem is to parametrize the uncertainties in f (v) and scan (or marginalize) over the associated parameters [12][13][14][15][16][17] (for alternative strategies see [18][19][20]). This approach suffers from the problem that it is unclear whether the chosen parametrization of the halo is sufficiently general. Moreover, direct detection experiments do not probe f (v) directly, but instead probe the velocity integral g(v min ) = v>vmin f (v)/v d 3 v. Therefore, it is typically necessary to perform large numbers of numerical integrations over the DM velocity distribution, making this method numerically slow. Efficient scans then often require a Bayesian approach, with the need to motivate prior distributions for DM parameters.
In this letter we propose a new method for dealing with uncertainties in the DM velocity distribution. Instead of parametrizing f (v), we directly parametrize g(v min ), so that predicted event rates depend on the parameters in a very simple way. This approach allows us to have a very large number of free parameters and removes the need to make any assumptions about the form of f (v). Consequently, any conclusions drawn from it will be completely robust in the face of astrophysical uncertainties. Moreover, our method involves a frequentist rather than a Bayesian approach, so that no prior distributions for any of the parameters need to be proposed.
General framework for direct detection. The differential event rate in a given experiment with respect to recoil energy E R is given by: where F (E R ) is the appropriate nuclear form factor (taken from [21]) and µ nχ is the reduced DM-nucleon mass. Furthermore, with A and Z being the mass and charge number of the target nucleus and f p /f n denoting the ratio of DMproton to DM-neutron couplings. The minimum velocity required for an energy transfer of size E R is v min (E R ) = m N E R /(2µ 2 ), where m N is the mass of the target nucleus and µ is the reduced mass of the DM-nucleus system. The rescaled velocity integral [20,22] is given by: where σ n is the DM-neutron scattering cross section, m χ is the DM mass, f (v) is the local DM velocity distribution in the lab frame and we take ρ = 0. [E i , E i+1 ] is given by [40]: where eff (E R ) is the effective exposure after including detector acceptance and is the detector response function [23].
Finding the optimum velocity integral. We observe from Eq. 1 that differential event rates depend on f (v) in the same way for all experiments, namely through the monotonically decreasing functioñ g(v min ) [41]. When discussing astrophysical uncertainties it is therefore more efficient to work directly with the velocity integral [20,22,24].
Our goal is to find the functiong(v min ) that best describes a given set of data from direct detection experiments for a proposed DM model. We want to make no assumptions on the functional form ofg(v min ) apart from it being monotonically decreasing [25]. We therefore proceed as follows: • We determine the region of v min -space probed by the experiments under consideration by converting recoil energies into v min -values. We then divide this region into N s evenly spaced intervals of the form [v j , v j+1 ].
• We take the ansatz thatg(v min ) is a monotonically decreasing sum of N s steps, i.e. we introduce N s constants Fig. 1). The monotonicity requirement implies that 0 ≤ g j ≤ g j−1 for all j [42].
• With this ansatz, and defining E j implicitly by v j = v min (E j ), Eq. 3 can be written in the simple form are calculable constants that depend only on the experimental details and the assumed DM properties, but are independent of astrophysics.
• We now define the usual χ 2 test statistic where N i is the experimentally observed number of events in the ith bin.
• Now, we numerically find the step heights g j that minimize χ 2 (g), maintaining the requirement of monotonicity 0 ≤ g j ≤ g j−1 . Although N s can be rather large (see below), it is nevertheless easily possible to find the global minimum of χ 2 (g). The reason is that the allowed region of g is convex and the Hessian of χ 2 (g) is positive semi-definite everywhere within. Consequently, any local minimum of χ 2 (g) is automatically a global minimum.
• Finally, we take N s → ∞. In practice, we take N s so large that further increases yield negligible improvements to χ 2 .
In typical examples of interest, we find that N s 30 is sufficient for this purpose. For such values of N s the minimization takes roughly a few seconds on a standard desktop computer.
Effectively, we decompose the DM velocity distribution into a large number of streams with different densities and speeds. Since any continuous function may be approximated arbitrarily well by a sum of step functions, this method effectively finds the best possible form for g(v min ) for a given set of data and model parameters.
Interestingly, we have found that, even in the limit that N s → ∞, the number of non-zero steps in the optimized halos (i.e. the number of steps with g j = g j+1 ) always remains smaller than the number of bins of data. This follows from the fact that, if we divide an optimized g(v min ) into "flat" sections with differing heights h i , then these must all satisfy ∂χ 2 /∂h i = 0. It may be checked that these equations can only all be satisfied if either the number of flat sections is smaller than the number of bins, or if the predictions P i can be made to perfectly match the observations N i . The latter possibility, however, is extremely unlikely for realistic cases including Poisson fluctuations. As a result, taking the limit N s → ∞ does not actually lead to more and more steps being added to the optimized halos, but rather allows for finer adjustments of the endpoints of the flat sections. Taking large N s thus turns out to be a useful trick for determining the optimal endpoints for the flat sections in a numerically efficient way.
We can now repeat our procedure for different sets of DM parameters and thereby find the minimum of χ 2 (g), calledχ 2 , for every point in parameter space. The bestfit values for the DM parameters are then determined by finding the global minimum ofχ 2 . We can then define ∆χ 2 =χ 2 −χ 2 min , which we have confirmed to follow a χ 2 -distribution with the degrees of freedom given by the number of fitted DM parameters. Consequently, ∆χ 2 can be used to define confidence intervals as usual.
One particular example is shown in Fig. 2, where we have generated mock data for a set of future experiments (taken from [9], see below) using a non-standard DM velocity distribution that contains a 25% contribution from a dark disk with velocity dispersion and lag relative to the baryonic disk both equal to 50 km s −1 [6]. For this example, we have assumed m χ = 800 GeV, f p /f n = 1 and σ n = 7 × 10 −45 cm 2 , and neglected Poisson fluctuations.
We observe that the best-fit values for the DM parameters obtained from our method (purple cross) perfectly coincide with the true values (white star) and the bestfitg(v min ) matches perfectly the assumed velocity integral (see inset). For comparison, we also show the bestfit values and 90% confidence regions obtained from two alternative methods, namely simply assuming the SHM with fixed parameters and optimizing the fit over σ n only (green) or assuming the SHM, but optimizing the fits over both σ n and the velocity dispersion σ dis (orange). It may be seen that both of these methods incorrectly exclude the true DM parameters at 90% confidence level.
Example: Hypercharged Dark Matter. We will now focus on a concrete model to show how our method is used in practice. If DM carries Standard Model hy-percharge (and is in an appropriate representation of SU (2) L so that it has an electrically neutral component), then it will generically interact with nuclei via tree-level Z-boson exchange, which results in a coupling ratio of f p /f n ∼ −0.04. As discussed in [26], current direct detection constraints on such DM particles require masses greater than about 10 8 GeV, and future experiments will be able to see a signal for m χ up to 10 10 − 10 11 GeV. Such heavy DM particles most simply obtain an appropriate relic abundance by having masses close to the reheating temperature of the universe. Since the coupling strength of hypercharged DM is fixed, detection of a signal would immediately reveal the DM mass through the scattering rate, and this would then give otherwise unobtainable information about the thermal history of the universe. In this scenario it would be of crucial importance to confirm that DM-nucleon scattering is mediated by Z-bosons. We will therefore use hypercharged DM as an example to show how our method may be used to determine f p /f n in a halo-independent way.
The most likely target materials for ton-scale future direct detection experiments are Xenon, Germanium and Argon [9]. Since their neutron to proton number ratios differ by less than 15%, determining f p /f n by looking at the relative rates of a signal on these elements will not be easy. Moreover, planned experiments using these elements will probe different regions of v min -space, not only because the energy thresholds may differ, but also because heavier target nuclei are sensitive to smaller velocities.
Consequently, astrophysical uncertainties severely affect the determination of f p /f n . If, for example, future data shows an excess of events in Xe-based experiments compared to Ar-based experiments, this observation could either be due to preferential coupling of DM to neutrons, as in the hypercharged scenario, or due to the DM velocity distribution decreasing more rapidly than in the SHM. Nevertheless, as long as there is non-negligible overlap between the regions of v min -space probed by different experiments, it will be possible to determine f p /f n in a halo-independent way, given sufficient exposure.
This situation is illustrated in the top panel of Fig. 3. The data points correspond to hypothetical measurements from future experiments employing Xe (blue), Ge (green) and Ar (purple) targets, assuming fermionic DM with hypercharge 1/2 and m χ = 7 × 10 7 GeV, compatible with the constraint from LUX [27], as well as a DM velocity distribution given by the SHM. For this plot, we have chosen a bin width of 10 keV and taken all experimental details from [9], assuming an energy-independent effective exposure. For Xe we only include the lowest 6 bins, because its rapidly decreasing form factor cuts off the event rate at higher energies. We map the experimental data ontog(v min )-values bin by bin using Eq. 1 [20] and test the incorrect model hypothesis that f p /f n = 1.
Because of this incorrect hypothesis the data points from Xe predict values ofg(v min ) that are large compared to the SHM prediction with best-fit normalization (orange dashed line), while the Ar predictions are relatively small. The optimum velocity integralg(v min ) as obtained from our method (red solid line), therefore clearly deviates from the SHM prediction and falls more steeply to reproduce this trend. In the absence of Poisson fluctuations, the value ofχ 2 associated to the best-fit velocity integral in this example isχ 2 = 1.05. Had we made the "correct" hypothesis that f p /f n = −0.04, the best-fit velocity integral would be identical to the SHM prediction used to generate the data, givingχ 2 min ≈ 0 (see Fig. 1). The difference ∆χ 2 = 1.05 describes the extent to which f p /f n = 1 is disfavored by the data. In this particular example, the hypothesis is excluded at the 69% confidence level. Of course, Poisson fluctuations are expected to modify our conclusions. One possible example for the effect of such fluctuations is shown in the bottom panel of Fig. 3, together with the corresponding best-fit velocity integral. If we include Poisson fluctuations, we find a median exclusion of the hypothesis that f p /f n = 1 at the 66% confidence level. With a probability of about 22%, the fluctuations are such that we can exclude f p /f n = 1 with at least 90% confidence. In Fig. 4 we show the 90% confidence limits on f p /f n which may be obtained from various exposures at Xe (rescaling the exposures at Ge and Ar from [9] correspondingly [43]). The data is generated assuming the same parameters for hypercharged DM as above, and we have minimized χ 2 over both g and the hypothesized DM mass. In order to show what may be accomplished in a typical case, we have again ignored Poisson fluctuations for this plot.
To conclude this section, we note that the power to exclude f p /f n = 1 comes largely from the fact that Gebased experiments have overlap in v min -space with both Xe and Ar targets. In the absence of a Ge target, we find (neglecting Poisson fluctuations) ∆χ 2 < 0.001. We come to the important conclusion that essentially no information can be inferred about f p /f n in a halo-independent way when using only Xe and Ar targets.
Applications and future directions. The method we have presented is extremely general and given a data set can be used to determine the best-fit velocity distribution for a wide variety of possible model hypotheses. There is no obstacle to performing analyses of many more complicated particle physics models, such as inelastic [28] or exothermic DM [29,30] or DM with long-range [31] or momentum-suppressed [32,33] interactions. Similarly, our method could be applied to comparing annually modulating signals of DM, such as observed by DAMA [34] and CoGeNT [35], and possibly also to constraining the modulation fraction.
Another exciting prospect is to apply our method to experimental results that give contradictory information when interpreted in terms of the SHM (such as LUX, CoGeNT, SuperCDMS [36] and CDMS-Si [37]) to understand if a different DM velocity distribution can bring these experiments into agreement. Here, however, there are two important complications. The first issue is that, for experiments with unclear compatibility, one would be interested in determining the goodness of fit at the actual best-fit point in addition to the determination of confidence intervals in DM parameter space. Typically, the value of the χ 2 statistic at a best fit point follows a χ 2 distribution with number of degrees of freedom given by the number of observations (i.e. bins) minus the number of free parameters. Our requirement of monotonicity, however, makes the notion of the number of free parameters in our halo fits somewhat unclear. In particular, we know that it generally remains impossible for us to fit data arbitrarily well even in the limit of an infinite number ofg(v min ) steps. The second issue is that, for experiments observing small event rates, binning the data and using χ 2 -methods becomes unreliable. In principle, there is no obvious obstacle to using an unbinned extended maximum likelihood method [38] to determine the optimumg(v min ), but doing so would appear to make future progress on the goodness-of-fit question difficult. We leave these problems to future work.