Taking Halo-Independent Dark Matter Methods Out of the Bin

We develop a new halo-independent strategy for analyzing emerging DM hints, utilizing the method of extended maximum likelihood. This approach does not require the binning of events, making it uniquely suited to the analysis of emerging DM direct detection hints. It determines a preferred envelope, at a given confidence level, for the DM velocity integral which best fits the data using all available information and can be used even in the case of a single anomalous scattering event. All of the halo-independent information from a direct detection result may then be presented in a single plot, allowing simple comparisons between multiple experiments. This results in the halo-independent analogue of the usual mass and cross-section plots found in typical direct detection analyses, where limit curves may be compared with best-fit regions in halo-space. The method is straightforward to implement, using already-established techniques, and its utility is demonstrated through the first unbinned halo-independent comparison of the three anomalous events observed in the CDMS-Si detector with recent limits from the LUX experiment.


Introduction
Despite the extraordinary progress in understanding the fundamental forces and building blocks of the Universe, the particle nature of dark matter (DM) remains unknown. Unravelling this puzzle remains one of the key tasks facing particle physics. Among the many facets of this wide-ranging endeavor, interactions between familiar matter and DM particles are searched for in underground direct detection experiments. Specifically these experiments search for the energy which would be deposited when a DM particle strikes a nucleus. Progress in this direction has been rapid and experiments have been extremely successful in pushing to ever lower interaction strengths, providing much needed knowledge about this aspect of DM phenomenology.
A few experiments have observed some potential signals of DM scattering, such as the long-standing DAMA annual modulation [1], the CoGeNT excess and modulation [2,3], CRESST-II excess [4], and most recently the CDMS-Si excess [5]. However, results from other experiments which have not observed any excess over expected backgrounds have increasingly put tension on DM interpretations of the positive hints, with the recent result from the LUX experiment excluding the simplest possibility, spin-independent elastically scattering DM where the DM couples equally to protons and neutrons [6][7][8][9]. 1 A lesson learned from studying these past DM hints is that the interplay between signals and constraints at different detectors may depend heavily on the local velocity distribution of DM, making this unknown a particularly troubling (or in some cases useful) nuisance parameter [15][16][17]. To mitigate this uncertainty, methods which allow the comparison of scattering rates at different detectors irrespective of the DM velocity distribution were developed [18,19] 1 Further analysis of the LUX results reveal that DM interpretations of the CDMS-Si excess with unequal DM couplings to protons and neutrons [7][8][9][10] now face increased tension with the LUX results. Models with exothermic scattering [9,[11][12][13] are now also in considerable tension with the LUX results, however there is no tension between LUX and an interpretation of the CDMS-Si excess in terms of a DM sub-component such as exothermic double-disk dark matter [14]. and subsequently extended to treat detectors with multiple target nuclei [20], detector energy resolution effects [21], annual modulation signals [22], and inelastic DM scattering [23]. While these methods have been of great utility in scrutinizing past DM hints, they will also be extremely important if a true signal of DM scattering begins to emerge, since the initial stages of discovery would begin with a small statistical excess within a particular detector. The ability to compare this to limits from other detectors in a fashion independent of the DM velocity distribution will be critical in assessing the validity of a signal.
While the halo-independent methods are very effective in interpreting null results from DM searches in order to place unambiguous limits on the allowed scattering rates at other detectors, the interpretation of an emerging DM signal using current halo-independent methods is open to some ambiguities. The current methods require that candidate DM scattering events be grouped into bins of recoil energy. The total rate in each bin is then mapped into a halo-independent rate to be compared with the limits from other detectors. For many applications this method is appropriate, however for an emerging DM signal it is not ideal for the following reasons.
• State-of-the-art detectors achieve expected backgrounds which are very low, typically expecting O( 1) background events in the DM acceptance region. As each new experimental run often leads to less than an order of magnitude increase in sensitivity, an emerging DM signal will likely come in the form of a small number of events. Many more events may follow with further experimental runs, but it is unlikely that the discovery of DM will begin with a large number of events. The binning of a small number of events is undesirable, since it is ambiguous and introduces sensitivity to the choice of bins. Hence, methods which rely on binning will not be optimal in the early stages of DM discovery.
• Current and future DM direct detection technology typically achieves excellent energy resolution. As the uncertainty in the energy of each candidate DM scattering event is likely to be small, bins wider than the energy resolution can only lead to the loss of important information about each event, effectively reducing the interpreted resolution and efficacy of the detector. Ideally, as much information as possible about each event should be retained in any comparison between candidate DM events and constraints from other detectors. For an emerging discovery, halo-independent methods which do not rely on binning are desirable.
In this work a new halo-independent method for analyzing candidate DM events is proposed which, by the above arguments, would be useful in the early stages of a DM discovery, and beyond. This builds on previous methods and relies on well-known properties of the integral over the velocity distribution of DM. The method allows for candidate DM events to be interpreted as best-fit points, with associated confidence intervals, for the DM velocity integral. These best-fit points and confidence intervals are shown to hold over all possible DM halos, and are in this sense halo-independent. Once determined, the implied values of the DM velocity integral can then be compared to limits from other detectors, allowing a halo-independent comparison between candidate DM signals and null DM experiments, free from the need to bin events and the ambiguities this introduces.
In Sec. 2 the halo-independent methods are reviewed. The calculation of constraints from null experiments is reviewed in Sec. 2.1 and the new method for an unbinned haloindependent interpretation of candidate DM events is developed and clarified in Sec. 2.2.
Comparisons between positive signals and null results are discussed in Sec. 2.3. In Sec. 2.4, we describe how the halo-independent information for one specific DM mass may be simply and unambiguously mapped to other DM masses, avoiding the proliferation of limit plots and calculations. The reader only interested in a short explanation of how to apply the methods can proceed directly to Sec. 2.5 where all necessary calculation steps for setting limits and for interpreting signals are briefly set out. In Sec. 3 the new unbinned halo-independent methods are applied to the three anomalous events observed in the CDMS-Si detector and compared to the current constraints from XENON10 and LUX. Finally, in Sec. 4 conclusions and suggestions for areas of future development are presented. App. A contains a proof that our method works equally well for both the idealized case of perfect energy resolution and the more realistic case of finite experimental energy resolution.

Halo-Independent Analysis Methods
The differential event rate 2 at a direct detection experiment is where m χ is the DM mass, m n the nucleon mass, µ nχ the nucleon-DM reduced mass, σ n the DM-nucleon scattering cross-section, ρ χ the local density, N A is Avogadro's number, F (E R ) is the nuclear form factor which accounts for loss of coherence as the DM resolves sub-nuclear distance scales, C T (A, Z) = (f p /f n Z + (A − Z)) is the usual coherent DM-nucleus coupling factor, (E R ) is the detector efficiency, and G(E R , E R ) is the detector resolution function. The velocity integral is where f (v) is the DM velocity distribution, and v E is the Earth's velocity, both in the galactic frame. We ignore the small time-dependence introduced by the Earth's motion around the Sun. For elastically scattering DM the minimum DM velocity required to produce a nuclear recoil energy where µ N χ is the nucleus-DM reduced mass. As is now standard, the constant factors which are common to all DM detectors are absorbed into a rescaled velocity integral An observation critical to the halo-independent methods, first noted in [18,19], is that because the velocity integrand is positive definite,g(v min ) is a monotonically decreasing function of v min for any DM halo. This observation becomes very powerful in developing halo-independent methods for the comparison of multiple experiments, as now described.

Constrainingg(v min ) -null results
Before considering the possibility of positive DM search results it is worthwhile to first consider the case of null experiments which can be used to constrain the velocity integralg(v min ).
We follow the discussion of [18]. Once a specific value of the DM mass m χ is chosen it is possible to place limits on the velocity integralg(v min ). If, at some reference minimum velocity v ref , the velocity integral is non-zerog(v ref ) = 0 then, since the velocity integral is monotonically decreasing, the unique form for the velocity integral which minimizes the total number of events for a giveng(v ref ) = 0 is

Discoveringg(v min ) -positive results
The most sensitive, and arguably least ambiguous DM detectors, strive to keep backgrounds low enough that O(1) background events are expected in a given run. They also typically have excellent energy resolution, such that ∆E R /E R 1. These factors combined suggest that the initial emergence of a DM discovery will likely be in the form of a relatively small number, N O , of events observed at discrete energies E i . To establish the consistency of such a scenario it will be important to then compare this potential DM discovery with limits from other experiments, ideally in a context free of uncertainties in the DM halo. Clearly an optimal route is to compare constraints ong(v min ) from the null experiments (described in Sec. 2.1) with the non-zero values ofg(v min ) hinted at by the emerging DM signal.
For positive signals, all current methods require the ad-hoc choice of a set of energy bins and then the calculation of upper and lower limits on the signal within these bins using the observed events and estimated backgrounds. These energy bins and preferred rates in each bin are then converted into v min -space bins and preferred values ofg(v min ) in each bin, subject to the constraint that the velocity integral is monotonically decreasing. The problems with such a method are immediately apparent. For emerging hints the number of events in the energy range of the detector will be small and binning a small number of events is a statistically questionable exercise from the outset, open to ambiguities and introducing issues with bin choice. Also, if a detector has good energy resolution, then valuable information is lost by binning the data in bins much larger than the experimental resolution, reducing the efficacy of any interpretation of the DM hint. Most crucially, binning data in bins of width much greater than the experimental resolution may lead to misinterpretation of the halo-independent constraints on this DM hint. Conversely, choosing bins of width much smaller than the energy resolution would, in the limit of a small number of events, smear single events across bins.
Ideally, it would be possible to map an emerging DM hint tog(v min ) − v min space in a way which preserves as much information as possible. In the case of detectors with excellent energy resolution this is of the utmost importance. But even for detectors with poor energy resolution there is information in the positions of the events and maintaining that information means employing methods which avoid binning the data.

The Method
A method commonly used in fitting a model with free parameters to unbinned data is the extended maximum likelihood method [24] which is desirable over the standard likelihood method as the normalization of a given rate is taken into account. When applied to a DM direct detection experiment which has observed N O events, in the energy range [E min , E max ], the extended likelihood is where dR T /dE R contains signal and background components and is the total number of events expected for a given set of parameters. We may compare different parameter choices by considering the log-likelihood, L = −2 log(L) which is minimized for a good fit and grows with decreasing quality of fit. Discarding constants irrelevant to the fitting procedure we have Using the DM rate, in terms ofg(v min ), as presented in Eq. (2.1), and including a background component where the first term accounts for the (small) estimated backgrounds and the last term the DM signal. There now appears to be a barrier to calculating L since there are an infinite set of possible DM halos to consider as one must also make a choice of the form ofg(v min (E R )) not only at each event, but over the whole range of measurable energies since the total number of events is calculated as the integral over this energy range. For simplicity let us first consider the case with perfect energy resolution G(E R , E R ) = δ(E R − E R ). A given set of events corresponds to a set of N O hypothetical values of g i ≡g(v min (E i )) as well as the form ofg(v min (E R )) interpolating between theg i . However, Eq. (2.8) penalizes against the total number of events predicted, since L increases as N E increases. Thus, sinceg(v min (E R )) is monotonically decreasing, the best fit out of all possible DM halos is the one which minimizes the total number of events predicted in any interval E i−1 < E R < E i between events. This is accomplished by choosing a constant valuẽ Figure 1: A schematic representation of all halo possibilities forg(v min ). If an experiment observes a number of events consistent with DM scattering, in this case three events of energy E i , then hypothetical values ofg(ṽ i−1 < v min ≤ṽ i ) =g i may be chosen where the positions of the stepsṽ i are given by v min (E i ) in the case of perfect energy resolution, and are allowed to float as free parameters if the energy resolution is non-zero. The solid blue curve will always minimize the extended log-likelihood, both in the case of perfect energy resolution and also with resolution effects included as demonstrated in App. A. Conversely the dashed red curve corresponds to the worst possible fit out of all halos, which is infinitely bad if the velocity integral between v low and v 1 is taken to infinity. Here, v low (v high ) is the velocity that corresponds to the low (high) energy threshold of the experiment. To determine the range of halos implied by the DM candidate events the parametersg i andṽ i may be varied, consistently choosing the solid blue curve in the likelihood, in order to determine the best-fit values and confidence intervals forg i . This form ofg(v min ) is quite robust. Indeed, in App. A we prove using variational techniques that the best-fitg(v min ) is still a sum of N O step functions even in the case of a very general resolution function; the only difference is that the positionsṽ i of the steps may now shift to the right of their position in the scenario with perfect energy resolution, v i ≥ v min (E i ). Thus, in all cases of interest, the form of the velocity integral which minimizes the extended likelihood for N O observed events is a sum of at most N O step functions, 4 whose 2N O free parameters (heights and positions) may be determined numerically in a straightforward manner, or analytically in the case of perfect energy resolution.
To calculate the log-likelihood it helps to define the differential background rate evaluated at the energy of each event E i , and which is the differential scattering rate at each event E i . Here E i are the positions of the steps in the halo velocity integralg (written as a function of recoil energy E R ) satisfying E i = E i in the case of perfect energy resolution. Another useful quantity is which is simply the total number of DM events expected. 5 In terms of these quantities (which depend on theg i and theẼ j ) the extended log-likelihood now decomposes as,  15) contains all of the information required to find the best-fit values and confidence intervals for the DM halo integral. To find the best-fit valuesg i,min and the best-fit positions of the steps E i,min , the likelihood may be numerically minimized to find L min , subject to the monotonicity constraint which must be imposed for any DM interpretation i.e.g i,min ≥g i+1,min . The confidence intervals in eachg i , denoted ∆g ± i , may be found by determining the extremum values satisfying L(g i ± ∆g ± i ) = L min + ∆L, for some ∆L which is determined from the statistical confidence desired. The other values ofg j =i and the positions of the steps E j should also be allowed to vary when determining the extremum values. It should be noted that in determining the confidence intervals, the monotonicity constraint must still be imposed, thus in determining ∆g ± i the otherg i =i may not always take their best-fit values.
As an infinite number of possible halos have been discarded, one may wonder whether this method actually captures the full ranges forg i at the desired confidence level. For a giveng i , a non-minimal halo not saturating the monotonicity constraint, i.e. one for which g (v min (E i−1 < E R < E i )) >g i , would only increase the value of N T and therefore the loglikelihood, meaning that a smaller range ofg i would be allowed with respect to the global minimum of the likelihood. Thus, rather than testing all possible halos to determine the bestfit values of, and allowed range of, theg i one can instead make the the minimal (saturating) choice,g(v min (E i−1 < E R < E i )) =g i . The best-fit points found this way will be the best fit out of all possible halos and the confidence intervals ∆g ± i necessarily encompass the maximally allowed ranges. This means that the envelope of allowedg i captures all halos for which the extended likelihood is within ∆L of the minimum.

Comparing with null results
Although a DM hint may suggest non-zero values ofg i for each anomalous event, it is desirable to compare these values in a halo-independent way with constraints from detectors which do not observe a signal. As described in [18], and Sec. 2.1, the most conservative limit on the velocity integral at a specific value of v min = v ref , denotedg(v ref ), may be determined by considering limits on the functiong( Calculating limits ong(v min ) in this way, and the best-fit values and confidence intervals forg i suggested by a DM hint using the method above, leads to plots such as Fig. 2, showing experimental limits and the best-fit values and confidence envelope for the velocity integral. It should be emphasized that the envelope ofg(v min ) does not imply that any curve passing through the envelope will have a log-likelihood value of L ≤ L min + ∆L, but it does imply that there exists a curve which passes through any single point in the envelope within a confidence interval satisfying L ≤ L min + ∆L. Furthermore, no curve with L ≤ L min + ∆L lies outside the envelope.
The most important information in any such plot is the interplay between the limits curve and the preferred envelope in the velocity integral. Consider a point on the lowest boundary of the envelope ing(v min ) at a point v min , denotedg − (v min ). The halo which leads to this value ofg(v min ) at v min , but would predict the smallest possible number of events in any detector, corresponds tog(v min ) =g − (v min )Θ(v min − v min ). However, it is precisely this halo shape which has been constrained by the null experiment. Hence, if a single point along the lowest boundary of the preferred envelope forg(v min ) is excluded by a null experiment then there is no halo within ∆L of the minimum of the likelihood for which the hint could be consistent with the null experiment. In other words, it is excluded by the null experiment independent of any uncertainties in the DM halo.

Varying m χ
The halo-independent methods are clearly of great value in comparing experimental results whilst avoiding the significant uncertainties in the velocity distribution of the DM. One perceived weakness of this approach is that it appears the calculations must be performed under the hypothesis of a single DM mass, m χ , and to consider a different DM mass m χ the entire calculation must be repeated again, leading to a proliferation of plots when presenting the results. However, assuming the detector is built from a single material, once limits and best-fit velocity integrals have been calculated for a single DM mass m χ , it is simple to map them to the analogous quantities for a different mass m χ .
Let us first consider the energy of a scattering event. The minimum DM velocity required is given by Eq. (2.3) which, for a specific scattering energy, immediately gives the relationship between v min (E R ) for a DM mass m χ and v min (E R ) for a DM mass m χ , mapping a point on the v min axes for m χ to a point on the v min axes for m χ while preserving the ordering of the scattering events. It should be noted that this mapping is nucleusdependent, and shifts limits and hints from different detectors by differing amounts. Furthermore, as halo-independent limits and best-fit points are calculated assuming a flat velocity integral between any neighboring events, the total number of events predicted between any two events only changes by a global normalization factor. This normalization can be found from Eqs. This has important implications for the presentation of DM direct detection results: if new DM limits, or hints, are presented in plot form in the halo-independent framework for a specific DM mass, then a single plot alone contains all of the information required for all DM masses. Thus, if an experimental collaboration released such a plot it would be possible to study halo-independent limits for any DM mass. Even if many details of the experimental analysis are not publicly available, this would enable the robust application of the DM results to different halo-independent scenarios by external groups.
As the shift in the normalization affects allg equally, the same minimum value of the log-likelihood Eq. (2.15) will be found for any DM mass, and the halo-independent method contains no information on the preference of data from a single experiment for a specific DM mass. A preferred mass may only be determined by appealing to a specific halo, or through requirements on the upper limits on v min due to the galactic escape velocity, or by combining date from multiple experiments.

In a Nutshell
It is useful at this point to summarize the steps required to perform an unbinned haloindependent analysis with real data. To calculate exclusion contours ing(v min ) from a null experiment, it is only necessary to calculate limits in the usual way, with the exception that at a point v ref the usual velocity integral is replaced withg and an upper bound is calculated for the constantg(v ref ). This process, which was described in [18] is repeated for different values of v ref to build up an exclusion contour.
In the case of an experiment with good energy resolution which observes an excess of events over an expected O( 1) background events, it is only necessary to assume the velocity integralg(v min ) takes the form of at most N O step functions with undetermined heights and positions. It is then necessary to calculate µ i ,μ i , for each event and the total number of events predicted N T , where these quantities are defined in Eq. (2.11), Eq. (2.12), and Eq. (2.13). With these quantities in hand one simply varies the heights and positions of the steps to find the minimum of the sum L/2 = N T − N O i=1 log (μ i + µ i ) with the additional constraints that g i,min ≥g i+1,min . The uncertainty on these determinations, at a given confidence level, is given by finding the variations, ±∆g v min , which saturate L = L min + ∆L (also allowing the positions of the steps to vary) to construct an envelope of preferred values forg(v min ). The ranges ±∆g i encapsulate the full envelope of possibilities of all DM distributions which are monotonic and satisfy L = L min + ∆L.
Once these limits and hints have been calculated and compared for a particular DM mass they have effectively been compared for all masses, assuming the DM scatters elastically and the detector consists of a single target. A section of the lower boundary of the preferred halo envelope for CDMS-Si is incompatible with the null LUX results, meaning that there is no DM halo compatible with the LUX results for which the extended likelihood is within ∆L of the best-fit halo. If the DM-nucleon couplings are tuned to maximally suppress scattering on xenon, a DM interpretation is still inconsistent with the LUX results. The curve for the SHM is also shown, giving a good fit to the CDMS-Si data as well as a curve for the best-fit halo which minimizes the extended likelihood.
halo-independent unbinned comparison between the CDMS-Si excess and the recent LUX results.
The S2-only XENON10 analysis [25] is used, with the ionization yield Q y also taken from [25]. We take the detector resolution function G(E R , E R ) to be a Gaussian with energydependent width ∆E R = E R / E R Q y (E R ). The acceptance is 95%, and the exposure is 15 kg days. Yellin's 'Pmax' method [26] is used to set limits.
The LUX collaboration have recently announced results from the first run [6]. The estimated LUX background distributions are not yet publicly available, making a profile likelihood ratio (PLR) test statistic analysis impossible. In [27] it was shown that for light DM the vast majority of nuclear recoil events would actually lie below the mean of the AmBe and Cf-252 nuclear recoil calibration band. The reason for this is that for a given low S2 signal the S1 signal is likely to have appeared above threshold due to a Poisson fluctuation. As there are no events in the region expected for light DM scattering (or equivalently low energy events) the DM event detection efficiency provided in [6] can be used to calculate the total number of expected events for a light DM candidate and then a Poisson upper limit can be set for zero observed events. We find excellent agreement with the estimated limits from [8] and good agreement with the official LUX results for the light DM region.
For CDMS-Si three events were found in 140.2 kg days of data [5]. We take the detector resolution function G(E R , E R ) to be a Gaussian and assume a conservative detector resolution of 0.5 keV. The acceptance is taken from [5]. The background contributions are taken from [28] with normalization such that surface events, neutrons, and 206 Pb, give 0.41, 0.13, and 500 600 700 800 0.08 events respectively. The best-fit points are calculated following the method described in Sec. 2.2 and confidence intervals are calculated for a variation ∆L = 1, where L is the total log-likelihood.
In Fig. 2 we show the halo-independent constraints on an elastically scattering spinindependent DM scattering interpretation of the CDMS-Si events. There is tension between the CDMS-Si excess and the LUX results independent of the DM halo if the DM couples equally to protons and neutrons. The lower energy events may still be consistent with LUX, however with a reduced number of events the significance of the excess is reduced. Even when couplings are tuned to maximally suppress scattering on xenon [10], an elastically scattering DM interpretation of the highest energy event is in tension with the LUX results independent of any uncertainties in the DM halo velocity distribution. The best-fit halo is also shown and interestingly this takes the form of two step functions. Although there are three events, the Gaussian smearing leads to a best-fit halo which only has two steps, whereas in the case of perfect energy resolution there would be three steps.
Thus, independent of uncertainties in the DM halo, and free from uncertainties introduced by binning the three anomalous CDMS-Si events, a DM interpretation of this excess faces considerable tension with the LUX results. This tension is reduced if the DM-proton and neutron couplings are tuned to maximally suppress scattering on xenon, however even when exploiting this freedom there is still tension with the LUX results. Ag(v min ) curve is also shown for the SHM to demonstrate that the CDMS-Si events give a good fit to the SHM.
In Fig. 3 we show the result of using the mapping, Eq. (2.16) and Eq. (2.17), from the points in Fig. 2 for m χ = 9 GeV to curves for other masses, demonstrating that an exclusion curve, or best-fit points, for a single DM mass contains all of the information necessary to translate the curve of best-fit points in a halo-independent way for different masses. This confirms that the presentation of new experimental results in a halo-independent manner for a single DM is a very efficient way to communicate the halo-independent information.

Conclusions
The DM direct detection field continues to evolve rapidly. The richness and effectiveness with which this dark frontier is explored relies on multiple experiments and detection strategies being employed. If an experiment begins to observe events consistent with DM scattering it will be crucial to confirm or refute this possibility with a separate independent experiment which uses different techniques and a different target nucleus. Previously developed haloindependent methods significantly reduce the systematic errors in such a comparison by eliminating the uncertainties due to the unknown DM velocity distribution. In this work these methods have been extended to enable a halo-independent analysis of candidate DM events without having to resort to event binning, which is inappropriate for a small number events and for detectors with good energy resolution, as would be expected in the circumstances of an emerging DM discovery. This method was developed for the simplest scenario of elastically scattering DM, however it would be interesting to extend it to include non-minimal scenarios such as inelastic or exothermic DM, or non-isotropic scattering.
The method we have described uses the standard approach of minimizing the extended likelihood, which has the advantage of being a well known technique in the field and thus straightforward for experimental collaborations to implement. Furthermore, it has the feature that results from multiple experiments can be straightforwardly added to the likelihood to carry out a combined analysis, although we have not studied such combinations here. This is true for both excesses and limits. It would be interesting to see if other statistical techniques, which do not require binning, give similar results. In addition to being straightforward for the experimental collaborations to implement, and reducing one of the systematic uncertainties that plague the interpretation of their results, we reemphasize that this does not come at the expense of complicating the presentation of their results. For DM scattering elastically in a detector with a single target, the results need only be presented for a single DM mass as this contains all necessary information; the extension to other masses is straightforward to calculate from the results for a single mass. In addition this method provides a haloindependent analogue of the usual comparison between limits and preferred regions.
Finally, as a test example we applied our technique to the recent results from CDMS and LUX. In accordance with expectations an unbinned halo-independent analysis of the three anomalous CDMS-Si events shows that for elastic, spin-independent scattering the CMDS-Si events are in tension with the null results from the LUX detector. If a DM interpretation of the CDMS-Si excess is to be found with no tension from the LUX results, this analysis suggests it will require non-standard DM scenarios.

Note added
While this work was in final preparation Ref. [29] was made public. While also concerned with halo-independent analyses Ref. [29] is based on event binning and, as emphasized in [29], is only applicable in scenarios with a large number of events. The method presented here is motivated by studying emerging DM hints in the limit of a few events. Although of most utility in the case of few events, it is applicable for an arbitrary number of events.

A Optimal halos and finite energy resolution
For an experiment with finite energy resolution G(E R , E R ), one may worry that due to smearing effects, the halo integral which minimizes the log-likelihood is no longer a sum of step functions, but perhaps a more complicated function whose many free parameters preclude a simple numerical minimization of the kind described in Sec. 2.5. Here we present a proof to the contrary -for any physically reasonable resolution function, the only effects of smearing are to shift slightly the positions of the steps ofg(E R ) away from the measured energies E i , and possibly to merge some of the steps. In particular, the optimal halo integral is still a sum of at most N O step functions.
Although we have in mind Gaussian smearing, this analysis holds for any reasonable form of the resolution function. We define a physically reasonable resolution function G(E R , E R ) to have the following properties: (ii) As a function of E R for fixed E R , G(E R , E R ) has a single local maximum at E R = E R and no other local extrema.
Property (i) simply states that the resolution function is normalized and doesn't change the total number of events. Property (ii) states that G has a single peak where the detected energy equals the true energy, and no other structure. Property (iii) is a technical assumption which will be used in the arguments below, and states that if G is flat on some interval, it must vanish. A normalized Gaussian resolution function G(E R , E R ) ∝ e −(E R −E R ) 2 /2σ 2 clearly satisfies all three properties, 6 as does a delta function G(E R , E R ) = δ(E R − E R ). Certain models for energy resolution may violate property (iii), for example a "top hat" shape where G(E R , E R ) is constant in some interval about E R and zero everywhere else, but one may always assume some infinitesimal deviations from flatness which otherwise has no measurable effect.
To simplify the notation, we write the differential scattering rate (2.1) as where we have absorbed the form factor, efficiency, and all prefactors into K(E R ). For reasonable choices of the form factor and efficiency functions, K(E R ) > 0 for all E R within 6 A Gaussian with an energy-dependent width σ(ER) also satisfies these properties as long as the form of σ(ER) is physically reasonable. For example σ(ER) ∼ 1/ √ ER in the XENON experiment, and G(ER, E R ) satisfies property (ii) as long as the region of ER close to zero is avoided. the experimental sensitivity, and in addition dK/dE R is small. We have also writteng(E R ) as a function of E R directly rather than v min . Note however that since v min is a monotonic function of E R ,g(E R ) is also monotonic function of E R . Consider now the expression for the log-likelihood (2.8), written in the suggestive form Property (i) above ensures the resolution function G does not appear in the first integral.
We can now view the log-likelihood minimization as a variational problem: minimize the functional L[g] with respect to the functiong(E R ), subject to the monotonicity constraint dg/dE R ≤ 0. The subject of variational problems with inequality constraints may be somewhat unfamiliar to physicists, but is well-known in economics and related fields; the solution is given by the Karush-Kuhn-Tucker conditions [30,31], which generalize the concept of Lagrange multipliers. In a similar fashion to imposing an equality constraint with a Lagrange multiplier, we can impose the inequality constraint by introducing an auxiliary function q(E R ) and modifying the log-likelihood, The solution that minimizes the log-likelihood while satisfying the monotonicity constraint will satisfy where is the total differential event rate at E i . But the left-hand side of Eq. (A.7) depends on E R through G(E i , E R ) while the right-hand side is constant; by property (iii), this is impossible. 8 This proves thatg(E R ) must be flat except at isolated points E j ; in other words, it is a sum of step functions.
To determine the number and position of the points E j , we can read Eq. (A.3) as a differential equation for q: The solution to this equation depends on the γ i , which in turn depend on the full solution functiong(E R ), so we cannot integrate this equation directly. In fact this turns out not to be necessary. By the complementarity condition, we must have q( E j ) = 0, but to preserve the positivity condition (A.5), we must also have at the roots E j of q. Taking the derivative of Eq. (A.9), and using the assumption dK/dE R ≈ 0, the condition on the second derivative becomes By property (ii), G(E i , E R ) will be peaked at E R = E i , so for E R close but not equal to E i , the ith term in the sum should dominate. The derivative of a peaked function is negative to the right of the peak, so to satisfy the inequality, we must have E j > E i . In other words, the positions of the steps shift to the right slightly.
The positions E j are given by solving dq/dE R | E j = 0, which we already used as Eq. (A.7) to derive the shape of g. Reading (A.7) now as an algebraic equation for E R , the left-hand side is a sum of N O peaked functions, weighted by various factors γ i . A function of this form will cross a horizontal line at most 2N O times, but only N O of these solutions will satisfy the convexity condition (A.11) and could qualify as roots of q. For sufficiently large γ i , the peak of G(E i , E R ) may dip below the line of height 1, so there may be fewer than N O solutions; the same is true if two peaks are close enough to one another to merge together. Furthermore, it may be the case that both conditions (A.10) are satisfied at E j , but q( E j ) = 0, in which case there would be no step at E j . We conclude that in the case of physically reasonable G(E R , E R ), the optimal halo integralg(E R ) is given by a sum of at most N O step functions.
Rather than integrate Eq. (A.9), one may simply use this knowledge to perform at most a 2N O -parameter numerical minimization of the log-likelihood subject to the monotonicity constraint ong.