Assessing Compatibility of Direct Detection Data: Halo-Independent Global Likelihood Analyses

We present two different halo-independent methods to assess the compatibility of several direct dark matter detection data sets for a given dark matter model using a global likelihood consisting of at least one extended likelihood and an arbitrary number of Gaussian or Poisson likelihoods. In the first method we find the global best fit halo function (we prove that it is a unique piecewise constant function with a number of down steps smaller than or equal to a maximum number that we compute) and construct a two-sided pointwise confidence band at any desired confidence level, which can then be compared with those derived from the extended likelihood alone to assess the joint compatibility of the data. In the second method we define a"constrained parameter goodness-of-fit"test statistic, whose $p$-value we then use to define a"plausibility region"(e.g. where $p \geq 10\%$). For any halo function not entirely contained within the plausibility region, the level of compatibility of the data is very low (e.g. $p<10 \%$). We illustrate these methods by applying them to CDMS-II-Si and SuperCDMS data, assuming dark matter particles with elastic spin-independent isospin-conserving interactions or exothermic spin-independent isospin-violating interactions.


Introduction
Astrophysical and cosmological evidence indicate that roughly 85% of the matter in the Universe is in the form of dark matter (DM) most likely composed of yet unknown elementary particles. Arguably the most extensively studied DM particle candidate is a weakly interacting massive particle (WIMP), which offers both theoretical appeal and hope for near-future detection. Most of the matter in our own galaxy resides in a spheroidal dark halo that extends much beyond the visible disk. Direct DM detection experiments represent one of the primary WIMP search methods currently employed. These experiments attempt to measure the recoil energy of nuclei after they collide with DM particles bound to the galactic dark halo passing through Earth. The current status of DM direct detection experiments remain ambiguous, with three experiments observing a potential DM signal and all others reporting upper bounds, some of which appear to be in irreconcilable conflict with the putative detection claims for most particle candidates [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18].
Interpreting the results of DM direct detection experiments typically requires assumptions on the local DM density, the DM velocity distribution, the DM-nuclei interaction, and the scattering kinematics. The uncertainties associated with these inputs can significantly -1 -

JCAP10(2016)029
affect the expected recoil spectrum (both in shape and magnitude) for a particular experiment, as well as the observed compatibility between experimental data. Attempts have been made to remove the astrophysical uncertainty from direct DM detection calculations, and compare data in a "halo-independent" manner, by translating measurements and bounds on the scattering rate into measurements and bounds on a function we will refer to asη(v min , t) common to all experiments, which contains all of the information on the local DM density and velocity distribution (see e.g. ).
The functionη(v min , t) depends on the time t and a particular speed v min . The physical interpretation of v min depends on the type of analysis being used. If the nuclear recoil E R is considered an independent variable, then v min is understood to be the minimum speed necessary for the incoming DM particle to impart a nuclear recoil E R to the target nucleus (and thus it depends on the target nuclide T through its mass m T , v T min = v min (E R , m T )). This has been the more common approach [19,21,24]. Alternatively, one can choose v min as the independent variable, in which case E T R is understood to be the extremum recoil energy (maximum for elastic scattering, and either maximum or minimum for inelastic scattering) that can be imparted by an incoming WIMP traveling with speed v = v min to a target nuclide T . Note that for elastic scattering off a single nuclide target the two approaches are just related by a simple change of variables. We will choose to treat v min as an independent variable for the remainder of this paper, as this choice allows us to account for any isotopic target composition by summing terms dependent on E T R (v min ) over target nuclides T , for any fixed detected energy E .
Early halo-independent analyses were limited in the way they handled putative signals. Only weighted averages on v min intervals of the unmodulated component ofη(v min , t), η 0 (v min ), and of the amplitude of the annually modulated component,η 1 (v min ), (see eq. (2.15) below) were plotted against upper bounds in the v min −η plane (see e.g. [19,21,22,27]). This type of analysis leads to a poor understanding of the compatibility of various data sets.
Recently, attempts have been made to move beyond this limited approach of taking averages over v min intervals by finding a best fitη 0 function and constructing confidence bands in the v min −η plane [31,46], from unbinned data with an extended likelihood [47]. One can then compare upper bounds at a particular confidence level (CL) with a confidence band at a particular CL to assess if they are compatible (see [46] for a discussion). From now on, when an upper index 0 or 1 is not written,η(v min ) is understood to beη 0 (v min ).
An alternative approach to analyzing the compatibility of data has been studied in [36] using the "parameter goodness-of-fit" test statistic [48,49] derived from a global likelihood (an alternative approach is taken in [37]). In [36], the compatibility of various experiments within a particular theoretical framework was determined by obtaining a p-value from Monte Carlo (MC) simulated data, generated under the assumption that the true halo function is the global best fit halo function. This approach has an advantage in that one can make quantitative statements about the compatibility between the observed data given a dark matter candidate model. However, this procedure assigns only a single number to the whole haloindependent parameter space, and we would like to have the ability to assess compatibility of the data with less restrictive assumptions on the underlying halo function.
In this paper we extend the approaches of [36] and [46] by using the global likelihood function to assess the compatibility of multiple data sets within a particular theoretical model across the halo-independent v min −η parameter space. This is done with two distinct approaches. First, we extend the construction of the halo-independent pointwise confidence band presented in [46] to the case of a global likelihood function, consisting of one (or more) The resultant global confidence band can be compared directly with the confidence band constructed from the extended likelihood alone, to assess the joint compatibility of the data for any choice of DM-nuclei interaction and scattering kinematics. The drawback of this method is that it cannot quantitatively address the level of compatibility of the data sets. To address this concern we also propose an extension of the parameter goodness-of-fit test, which we will refer to as the "constrained parameter goodness-of-fit" test, that quantifies the compatibility of various data sets for a given DM particle candidate assuming the halo functionη(v min ) passes through a particular point (v * ,η * ). By calculating the p-values for each (v * ,η * ) throughout the v min −η plane, one can construct plausibility regions, such that for any halo function not entirely contained within the plausibility region the data are incompatible at the chosen level, e.g. p < 10%.
In section 2 we review the procedure for constructing the best fit halo functionη BF and confidence band from an extended likelihood. Readers familiar with [46] may wish to skip this section and go directly to section 3, which discusses how the construction of the best fit halo function and confidence band is altered when dealing with a global likelihood function that is the product of one (or more) extended likelihoods and an arbitrary number of Poisson or Gaussian likelihoods. In section 4, we use the methods discussed in section 3 to construct the best fit halo and global pointwise confidence band, for the combined analysis of CDMS-II-Si and SuperCDMS data assuming elastic isospin-conserving [50][51][52] and exothermic isospinviolating spin-independent (SI) interactions [32,35]. Section 5 introduces the "constrained parameter goodness-of-fit" test statistic and the construction of the plausibility regions. This method is illustrated using CDMS-II-Si and SuperCDMS data, assuming elastic isospinconserving spin-independent interactions. We conclude in section 6.
2 Review of the Extended Maximum-Likelihood Halo-independent (EHI) analysis method

Generalized halo-independent analysis
The differential rate per unit of detector mass as a function of nuclear recoil energy E R for dark matter particles of mass m scattering off a target nuclide T with mass m T is given by where ρ is the local dark matter density, C T is the mass fraction of the nuclide T in the detector, f (v, t) is the dark matter velocity distribution in Earth's frame, and dσ T /dE R is the WIMP-nuclide differential cross section in the lab frame. When multiple target elements are present in the detector, the differential rate is To allow for the possibility of inelastic DM-nuclei scattering, we consider a DM particle scattering to a new state of mass m = m + δ, where |δ| m, and δ > 0 (< 0) describes endothermic (exothermic) scattering. In the limit µ T |δ|/m 2 1, v min (E R ) is given by , that can be imparted to a target nucleus by a DM particle traveling at speed v in Earth's frame, given by Eq.
(2.4) shows that for endothermic scattering there exists a nontrivial kinematic endpoint, given by the DM speed v T δ = 2δ/µ T , below which incoming DM particles cannot induce nuclear recoils. When multiple targets are present in a detector, we use v δ to denote the minimum of all v T δ . For exothermic and elastic scattering v δ = 0. Experiments do not actually measure the recoil energy of a target nucleus, but rather a proxy for recoil energy (e.g. the number of photoelectrons detected in a photomultiplier tube) denoted E . The differential rate as a function of the detected energy E is given by where the differential rate in eq. (2.1) has been convolved with the efficiency function (E R , E ) and the energy resolution function G T (E R , E ), which together give the probability that a detected recoil energy E resulted from a true recoil energy E R .
Upon changing the order of integration, one can express the differential rate in detected energy as where dH T /dE is given by (2.7) and we define dH Here, we only consider differential cross sections that depend on the speed of the WIMP v = |v|. The cross section depends only on the speed v if the incoming WIMPs and the target nuclei are unpolarized and the detector response is isotropic, as is most common. In eqs. (2.6) and (2.7), we have incorporated the parameter σ ref which denotes the overall strength of the interaction. For example in the case of the SI interaction, with differential cross section given by where A T and Z T are the atomic and charge numbers of nuclide T , f n and f p are the neutron and proton couplings, and F T (E R ) is the form factor normalized to F T (0) = 1 (taken here to be Helm form factor), we will choose σ ref = σ p , the WIMP-proton cross section.

JCAP10(2016)029
A halo-independent analysis relies on the separation of the astrophysical parameters from the particle physics and detector-dependent quantities. Here we follow [27]. Let us definẽ Differentiating both sides of eq. (2.10) gives which upon insertion into eq. (2.6) leads to Using the fact thatη(∞, t) = 0 and dH/dE (E , v δ ) = 0, integration by parts of eq. (2.12) results in the following expression for the differential rate where we have now defined the differential response function dR/dE as (2.14) η(v min , t) is a function of time due to the annual rotation of the Earth around the Sun. If one now makes the approximatioñ η(v min , t) η 0 (v min ) +η 1 (v min ) cos(2π(t − t 0 )/year) (2.15) and integrates the differential rate over the energy range of interest, the unmodulated component R 0 and annual modulation amplitude R 1 of the rate are given by where α = 0 or 1, and the energy-integrated response function R is given by In the event that R [E 1 ,E 2 ] (v min ) is a well-localized function in v min , measurements on unmodulated and modulated rate can be used to infer the average values ofη 0 andη 1 over a v min interval. This is the case for DM-nuclei differential cross sections proportional to 1/v 2 (e.g. the typical SI and SD contact interactions). Should the differential cross section not be of this form, one may need to regularize the energy-integrated response function as described in [27].

Extended maximum likelihood analysis
It was initially proven in [31], that if there is no uncertainty in the measurement of recoil energies in a single nuclide target, then the extended likelihood, given by is maximized by a non-increasing piecewise constantη 0 (v min ) function (which we call simplỹ η(v min )) with at most N O (the number of observed events) steps. N E [η] in eq. (2.19) is the total number of expected events, and E a is the observed energy of event a. This proof was generalized to the case of realistic energy resolution and arbitrary target composition in [46]. The generalized proof presented in [46] applies the Karush-Kuhn-Tucker (KKT) conditions, which are only valid for systems with an objective function of finite number of variables subject to a finite number of constraints, to the likelihood functional in eq. (2.19) by discretizing the variable v min , applying the KKT conditions, and then taking the continuum limit.
Here, we will briefly review the conclusions presented in [46]. If one defines the quantity then instead of maximizing the likelihood, one can equivalently minimize L[η]. The KKT conditions, applied to eq. (2.20) and taken in the continuum limit, lead to the following: A direct consequence of eq. (2.24) is thatη(v min ) is a piecewise constant function with the locations of the steps given by the v min values which satisfy q(v min ) = 0. For this reason, we need to analyze the behavior of q(v min ). Eq. (2.19) and (2.21) can be used to show that where we have defined the following quantities:

JCAP10(2016)029
For the extended likelihood function in eq. (2.19), the behavior of the terms in eq. (2.25) were studied in [46] to determine how many steps can appear in the best fitη function. We briefly review their behavior here (see [46] for additional details).
Consider first the v min -dependence of dH/dE (see eq. (2.7)), which appears in both the integrand of ξ(v min ) and in H a (v min ). If the differential cross section is proportional to v −2 , as is the case for the standard SI and SD contact interactions, the only velocity dependence of dH/dE is in the integration range . For these interactions, as v min increases, the integration covers a larger portion of the parameter space where the integrand is non-zero. At large values of v min , the entire region where the integrand is non-zero is included in the integration and dH/dE becomes constant. For a fixed value of E , one would expect the integrand of dH/dE to be a well-localized function of E R (i.e. an observed recoil E can only result from a narrow range of E R values). For this reason, the terms H a (v min ) appear as step-like functions in v min .
The term ξ(v min ) contains an additional integration of dH/dE over E . The only dependence on E appears in the factor (E , E R )G T (E , E R ), which describes the probability a detected recoil energy E is the result of some true recoil energy E R . For small values of v min , only a narrow range of recoil energies are integrated over and thus ξ will be quite small (i.e. , for v min values such that E T,+ R (v min ) is below threshold). As v min increases, the integration range widens and ξ(v min ) steadily increases. Eventually, the entire region where the integrand is nonzero is included in the integration, and ξ(v min ) becomes constant.
The only term dependent on the halo function is γ a [η], which only alters the relative contribution of each step-like function to q(v min ). The functionη BF (v min ) can only be discontinuous when q(v min ) = 0, which is equivalent to saying the steps ofη BF occur where the step-like functions H a (v min )/γ a [η] touch ξ(v min ) from below. Since there is a single term of the form H a (v min )/γ a [η] for each observed event, the number of steps appearing inη BF must be less than or equal to the number of observed events, N O .

Construction of the best fit halo function and confidence band from an extended likelihood
In this section we briefly review the construction of the best fit functionη BF (v min ) and the confidence band for an extended likelihood [46]. Let us define the function where v = (v 1 , . . . , v N O ) and η = (η 1 , . . . ,η N O ), and the various v a andη a specify the location and height of each step. Here, we have defined the piecewise constant functionη N O as Using the result of the previous section, minimizing the functional L[η], and thus finding the best fitη(v min ), is now reduced to minimizing We can define the confidence band as the region filled by all possibleη functions satisifying where L min is the minimum of L[η], and ∆L * corresponds to the desired confidence level. However, in practice, finding allη functions satisfying eq. (2.32) is not possible. Instead, let us consider the possible subset ofη functions which minimize L[η] subject to the constraint Now let us define L c min (v * ,η * ) to be the minimum of L[η] subject to the constraint in eq. (2.33), and If the point (v * ,η * ) lies within the confidence band, then there should exist at least oneη function passing through this point which satisfies ∆L[η] ≤ ∆L * . Should this be the case, it follows that ∆L c min (v * ,η * ) ≤ ∆L * . Alternatively, if ∆L c min (v * ,η * ) ≥ ∆L * , one can state that there does not exist a singleη which satisfies ∆L[η] ≤ ∆L * . Thus the confidence band can be constructed by finding the values of (v * ,η * ) which satisfy ∆L c min (v * ,η * ) ≤ ∆L * . This condition defines a two-sided interval aroundη BF for each v min value (with v min = v * ), and the collection of those intervals forms a pointwise confidence band in v min -η space, which we are simply calling the confidence band.
To understand the meaning of ∆L c min , let us first discretize the continuous variable v min into a collection of K discrete values v min = (v 0 min , . . . , v K−1 min ). The likelihood functional in eq. (2.20) then becomes a function of the K−dimensional vector η = (η 0 ,η 1 , . . . ,η K−1 ) which defines the piecewise constant functionη(v min ; η) given bỹ With this discretization, the constraint on (v * ,η * ) in eq.
whereη i are theη i values which maximize the likelihood function L(η 0 , . . .η K−1 ) ≡ L[η(v min ; η)] subject to the constraintη k =η * , andη i maximize L without the constraint. ∆L k,c min (η * ) now defines the −2 ln of the profile likelihood ratio with one parameter (η k ), and thus by Wilks' theorem the distribution of ∆L k,c min (η * ) approaches the chi-square distribution with one degree of freedom in the limit where the data sample is very large. If we now recover the continuum limit by taking K → ∞, we see that ∆L k,c min (η * ) approaches ∆L c min (v * ,η * ). Thus the construction of the confidence band is equivalent to finding the collection of confidence intervals inη * for each v * at a given CL for which ∆L c min < ∆L * . Assuming that ∆L c min is chi-square distributed, the choices ∆L * = 1.0 and ∆L * = 2.7 correspond to the confidence intervals ofη at 68% and 90% CL, respectively, for each v min value. In [46] it was shown that the constrained best fit halo functionη c BF defining L c min (v * ,η * ) is a piecewise constant function with at most N O + 1 steps, with the additional step potentially appearing at (v * ,η * ). An in-depth discussion of the interpretation of the confidence band constructed from the profile likelihood ratio is provided in [46].

JCAP10(2016)029 3 Extension of EHI analysis to a global maximum likelihood
In this paper we extend the analysis presented in [46] to make statistically meaningful statements about the data of multiple experiments in a halo-independent manner. Specifically, we (i) extend the formalism of constructing a pointwise confidence band from a profile likelihood in halo-independent parameter space to a global likelihood function (this section), and (ii) propose a method for creating plausibility regions, constructed from a new family of test statistics which can assess the compatibility of multiple data sets under the assumption that the halo functionη(v min ) passes through each (v * ,η * ) point (see section 5). To accomplish these tasks one must first understand how to find the best fit halo function and constrained best fit halo function from a global likelihood.
In this section we extend the procedure of [46] to the global likelihood function, defined by the product of some number N exp of individual likelihood functions, α = 1, 2, . . . N exp , The procedure of [46] relies on the fact that an extended likelihood function is maximized by a non-increasing piecewise constantη BF (v min ) function with a finite number of points of discontinuity. As discussed below, the methods and reasoning of [46] extend to a global likelihood, if it includes at least one extended likelihood. Thus, the global likelihood function we will work with for the remainder of the paper is where L EHI is an extended likelihood (EHI stands for "extended halo-independent" [46]) as in eq. (2.19) and, for each α, L α represents Poisson likelihoods, or Gaussian likelihoods j , and n We now prove that global likelihoods of the form eq. (3.2) are maximized by nonincreasing piecewise constantη functions with at most N steps, The KKT conditions in eq. (2.21)-(2.24) apply equally to any likelihood function L. The KKT condition in eq. (2.24) implies thatη BF is constant in an open interval where q(v min ) = 0. Thus if the q(v min ) function given by eq. (2.21) has only a finite number of isolated zeros within a range, the best fitη in this range should be a piecewise constant function with steps located at the zeros of q(v min ). Therefore, the problem of determining the potential number of steps ofη BF is equivalent to counting the maximum possible number of isolated zeros of the q(v min ) function.
For the global likelihood in eq. (3.2), q(v min ) is given by for Poisson likelihoods of the form in eq. (3.3), and for Gaussian likelihoods in eq. (3.4). Changing the functionη(v min ) only alters the sign and magnitude of the prefactor of ξ The v min dependence of Q (α) [η, v min ] exclusively appears in the functions ξ In appendix A we prove that above a certain value of v min , given by the minimum v µ low (see appendix A.1 for definition), the zeros of q(v min ) in eq. (3.6) are isolated, and the maximum number of isolated zeros is given by eq. (3.5). However, in practice the number of steps is smaller than N and can be determined by studying the functional form of the functions ξ EHI (v min ), H EHI a (v min ), and ξ (α) j (v min ) (which are independent ofη). In appendix B we prove the uniqueness of the best fit halo function,η BF .
An explicit example of the q(v min ) function and its components is shown in figure 1 for the case of CDMS-II-Si combined with SuperCDMS data. For SuperCDMS we have taken a one-bin Poisson likelihood, summing over all detectors in table 1 of [18], the contribution from which to q(v min ) is shown in green. Also included in figure 1 are the contributions to q(v min ) arising from ξ EHI (v min ) (red) and the summation over the H a (v min )/γ a [η] (blue). Figure 1 shows that q(v min ) goes to 0 at v min 510 km/s and 580 km/s, denoting the locations of the steps ofη BF (shown later in figure 2).
We would like to emphasize that all of the aforementioned arguments have relied on having a global likelihood that contains at least one extended likelihood. This likelihood has the essential feature of contributing anη-dependent term and anη-independent term to q(v min ), with different functional dependences on v min .  [46]. The red crosses represent the 68% CL intervals of the averagedη arising from binning the CDMS-II-Si events into 2 keVnr bins between 7 and 13 keV (see e.g. [25,28,34], in which we take the horizontal bars to be the v min range where 90% of the area under is contained). The results shown assume a 9 GeV DM particle scattering elastically through a SI isospin-conserving contact interaction (f n /f p = 1).

JCAP10(2016)029
In order to construct the two sided confidence band, we compute at each value of v min = v * the two sided interval defined by whereL G (v * ,η * ) is the maximum of the global likelihood subject to constraint eq. (2.33), andL G is the maximum of the global likelihood. Using the same arguments of section 2.3 and assuming that ∆L c G,min is chi-square distributed, the distribution of ∆L c G,min has one degree of freedom and ∆L * = 1.0 and ∆L * = 2.7 for the 68% and 90% CL intervals, respectively. In section 4.2 of [46] it was shown that if L is maximized by anη BF function with a maximum of N steps, then L(v * ,η * ) (i.e. L subject to the constraint thatη(v min ) passes through the point (v * ,η * )) is maximized by a halo function, which we call the constrained best fitη c BF , with a maximum of (N + 1) steps, one of which could occur at v min = v * . This proof applies to L G and L G (v * ,η * ) as well.

Global likelihood analysis of CDMS-II-Si and SuperCDMS data
Here we apply the formalism described in section 3 using the global likelihood function in eq. (3.2) with an extended likelihood [47] for the three events observed by CDMS-II-Si [7], and a 1-bin Poisson likelihood for SuperCDMS [18]. To obtain background estimates for CDMS-II-Si, we take the normalized background distribution functions from [53] and rescale them such that 0.41, 0.13, and 0.08 events are expected from surface events, neutrons, and 208 Pb respectively (see [7]). Since the resolution function for silicon in CDMS-II has not been measured, we take the energy resolution function for germanium from eq.1 of [54].
In addition to implementing the full SuperCDMS data in table 1 of [18] (11 events observed, 6.56 expected background events, 577 kg-days of exposure), we also use a subset of the SuperCDMS data which neglects the observed events (and the exposure) from tower 5 (4 events observed, 5.33 expected background events, 412 kg-days of exposure). The SuperCDMS collaboration acknowledges that tower 5 had a malfunctioning guard electrode which resulted in a poor understanding of the background in this tower. We will use the label "SuperCDMSLT5" for this analysis (where LT5 stands for "Less Tower 5").
The data analysis used throughout this paper is included in the CoddsDM software [55], an open-source Python program for the analysis of dark matter direct detection data.
In the left panel of figure 2 we show the 68% (dark red) and 90% (light red) CL confidence bands, calculated assuming ∆L c min (v * ,η * ) is χ 2 distributed with one degree of freedom, for the combined analysis of CDMS-II-Si and SuperCDMS, assuming a 9 GeV DM particle scattering elastically off nuclei with a SI isospin-conserving contact interaction. Also shown in figure 2 is the globalη BF function (dark red line), theη BF function for CDMS-II-Si data alone (blue line), the SuperCDMS 90% upper limit (dark yellow), and the upper and lower boundaries of the 68% (black dashed) and 90% (black solid) CL confidence bands obtained using CDMS-II-Si data alone (these coincide with those presented in figure 3 of [46]). Notice that the confidence bands are unbounded from above for v min 275 km/s and v min 400 km/s, for the global analyses and CDMS-II-Si analyses respectively (the lower boundaries of the confidence bands are, however, well defined asη(v min ) is a non-increasing function). This is because q(v min ) = 0 in these intervals (i.e. the experiment/experiments are not sensitive to recoils imparted from DM traveling at these speeds), and thus theη BF is actually undetermined.

JCAP10(2016)029
Since the purpose of plotting these functions is to compare the compatibility of putative and null signals, we extendη BF in our plots to this region, in the most conservative way (i.e. constant). The red crosses in figure 2 represent the 68% CL intervals (vertical bars) of averagedη over corresponding v min intervals (indicated by horizontal bars) arising from binning the CDMS-II-Si events into 2 keVnr bins between 7 and 13 keV (see e.g. [25,28,34], except we take the horizontal bars to be defined by the v min range where 90% of the area To determine the all upper bounds onη 0 arising throughout this paper from the Super-CDMS data, we follow the procedure first outlined in [19,21]. Using the fact thatη 0 (v min ) is a non-increasing function, this procedure argues the smallest possible function passing through a point (v 0 ,η 0 ) is the downward step-functionη 0 Θ(v 0 − v min ). With this in mind, eq. (2.16) can be rewritten such that an upper bound on the observed rate in the energy range [E 1 , E 2 ] can be translated into an upper boundη lim (v min ) onη 0 , using . (4.1) This limit is conservative in that everyη 0 function lying above the bound is excluded by the data, but not allη 0 functions lying below the bound are allowed by the data. The values of R lim used in this paper are determined using the Feldman-Cousins approach [56]. Assuming a Poisson distribution for both SuperCDMS (n = 11, b = 6.56) and SuperCDMSLT5 (n = 4, b = 5.33) and an energy range [E 1 , E 2 ] corresponding to the quoted experimental range (i.e. E 1 = 1.6 keV and E 2 = 10.0 keV), this leads to 90% CL upper limits on the number of DM events µ lim of 11.25 and 3.33 events respectively. The value of R lim can then be obtained by dividing µ lim by the exposure of the relevant experiment. The globalη BF function is shifted to lower values ofη by over an order of magnitude relative to theη BF found using CDMS-II-Si data alone, and is outside the 68% and 90% CL confidence bands of CDMS-II-Si alone. Similarly, theη BF for CDMS-II-Si alone (in blue) is incompatible with the 68% and 90% Cl global confidence bands. Furthermore, in the range 360 km/s v min 480 km/s the 68% CL global confidence band has no overlap with the 68% CL confidence band of CDMS-II-Si.
The right panel of figure 2 is the same as the left panel but using SuperCDMSLT5 instead of SuperCDMS. The globalη BF function has shifted to slightly lower values ofη (relative to the SuperCDMS analysis), as have both confidence bands, but the general conclusions are the same -namely, there appears to be a strong level of incompatibility between the results arising from the global likelihood and those found using only CDMS-II-Si data. We also note that the increased conflict between CDMS-II-Si and SuperCDMSLT5 has resulted in the 90% CL confidence band extending down toη 0 (i.e. no DM) at low values of v min , as opposed to having a well defined non-zero lower boundary for the case of SuperCDMS.
We present one final illustration of this method in figure 3 for a 3.5 GeV DM particle with exothermic scattering (δ = −50 keV) and a Ge-phobic SI interaction (f n /f p = −0.8) [32]. This example has been chosen to illustrate how the globalη BF and confidence bands behave in the case of non-conflicting data sets. As expected, the results from the global likelihood analysis of CDMS-II-Si and SuperCDMSLT5 are nearly identical to the results obtained from CDMS-II-Si alone, with the only significant change occurring at low values of v min , where the upper bound of SuperCDMSLT5 is in conflict with the confidence bands of CDMS-II-Si alone.

Constrained goodness-of-fit analysis
The global likelihood analysis presented in the previous section always produces a best fit halo function and confidence band, even when considering conflicting data sets. A particular goodness-of-fit test has been proposed in [48,49] to assess the compatibility of different data sets in the framework of a given theoretical model. This so called "parameter goodness-of-fit" (PG) test was used in [36] to gauge the compatibility of CDMS-II-Si, SuperCDMS, and LUX data, in a halo-independent way. It is defined as whereL G is the maximum of the global likelihood andL α is the maximum of the likelihood of experiment α. If theη BF of all individual experiments would coincide, then q P G = 0. On the other hand a strong disagreement between theη BF of individual experiments would lead to a large value of q P G . Thus q P G quantifies the degree of compatibility of all data sets under the assumption of a particular DM particle model. To provide a quantitative statement about the compatibility, the p-value of the observed data was obtained from a MC simulation, assuming the globalη BF is the true halo model [36]. This procedure assigns a single number, a single p-value, to the whole halo-independent parameter space, and we would like to identify regions of this space whereη(v min ) functions may lead to better or worse compatibility among data sets. With this purpose in mind, we define a family of test statistics similar to q P G , one for each point in parameter space, using the profile likelihood, defined as the likelihood maximized subject to the constraint in eq. (2.33), i.e.η(v * ) =η * (it is the continuum limit of the numerator inside the square bracket in eq. (2.36)). We will then define a p-value for every point in the halo independent parameter space. We define the -14 -  "constrained parameter goodness-of-fit" test statistic as whereL c G (v * ,η * ) is the global profile likelihood andL c α (v * ,η * ) is the profile likelihood of experiment α. q c P G tests the compatibility of the different data sets under the assumption that η(v min ) passes through (v * ,η * ). To infer the probability distribution for q c P G (v * ,η * ) we use a Monte Carlo simulation, assuming the true halo model is given by the global best fit halo function that maximizes L G under the constraintη(v * ) =η * . We call "constrained best fit halo function"η c BF the function that maximizes a likelihood subjected to this constraint. There is a differentη c BF (v min ) function for each (v * ,η * ) point (which certainly fulfills the conditioñ η c BF (v * ) =η * ), one for the global likelihood and one for each single experiment extended likelihood. The p-value for a given (v * ,η * ) is then obtained by comparing the observed value of q c P G to the distribution constructed from O(10 3 ) simulated data sets (for each choice of (v * ,η * )).
We have only developed a method for maximizing the Poisson and Gaussian likelihoods subject to the constraintη(v * ) =η * for a single bin Poisson/Gaussian likelihood. In this case, the likelihood is maximized by an expected number of dark matter eventsν 1 . In order to maximize the constrained likelihood, one needs to consider whether v * lies above or below the experimental threshold. If v * is below threshold, a halo function passing through (v * ,η * ) produces a minimum number of 0 observed events (withη =η * Θ(v * −v min )), and a maximum number ν max of events given by the flat halo functionη(v min ) =η * . If v * is above threshold, a halo function passing through (v * ,η * ) produces a minimum number ν min of observed events whenη =η * Θ(v * − v min ), and there is no limit on the maximum number of observed events becauseη can be unbounded from above for v min < v * . Ifν (α) j lies between the minimum and maximum number of predicted events forη(v min ) passing through (v * ,η * ) in each case, then the maximum of the constrained likelihood is the maximum of the likelihood. Otherwise, the maximum of the constrained likelihood is calculated using ν max or ν min , depending on the respective case above.
The probability distributions of q c P G are shown in figure 4 for the combination of CDMS-II-Si and SuperCDMSLT5, for a SI contact interaction with (m, δ, f n /f p ) = (9 GeV, 0 keV, 1), . Plausibility region (light purple) generated from the constrained parameter goodnessof-fit test statistic for CDMS-II-Si and SuperCDMSLT5 (p-value larger than 10%), compared with the confidence bands (red shaded) generated for CDMS-II-Si data alone (left) [46] and the global confidence bands (red shaded) constructed in section 4 (right). The plausibility regions are crossed over because halo functions entirely contained within these regions are not necessary allowed by our test, i.e. do not necessarily lead to a compatibility of the data sets at the level of p >10%. However, for any halo function not entirely contained within the plausibility region the data sets are incompatible at the chosen level (p <10%). Also shown areη BF for CDMS-II-Si alone (blue), theη BF resulting from the global likelihood analysis (dark red), and the v min -averaged CDMS-II-Si data (crosses) as described in section 4.
for v * = 400 km/s (left) and 500 km/s (right) withη * chosen on the globalη BF curve. The observed value of q c P G in figure 4 are indicated by the dashed red line. The p-values roughly correspond to 2.8% for v * = 400 km/s, and 0.5% for v * = 500 km/s. While the probability distributions shown in figure 4 do not appear to approach 1 in the limit x → 0, there are in fact a large number of simulations which yield extremely small values of q c P G that are not depicted (the probabilities do in fact equal 1 at x = 0). This happens because the global best fit halo function predicts less than one observed event in CDMS-II-Si, which leads to many simulations in which 0 events are observed by CDMS-II-Si. In turn, this implies the global constrained best fit halo function and the constrained best fit halo function for CDMS-II-Si are the same, as they can only have a single step at the location of (v * ,η * ). For SuperCDMSLT5, the expected background is larger than the number of observed events, and thus the profile likelihood of SuperCDMSLT5 is relatively insensitive to halo functions that predict small numbers of DM events. Consequently, it is not uncommon to find lnL c G (v * ,η * ) α lnL c α (v * ,η * ). Figure 4 already demonstrates a high level of incompatibility between the CDMS-II-Si and SuperCDMSLT5 data sets for the assumed WIMP candidate, because the global η BF (v min ) cannot produce a large p-value, say larger than 10%. We can construct intervals at each v min = v * in which the probability of obtaining a q c P G value larger than the one observed is ≥ 10%. By joining these intervals we build regions in (v min ,η) which are referred to as "plausibility" regions.
Let us now clarify the meaning of the plausibility regions. A halo functionη(v min ) is a non-increasing continuous function which must be defined for any value of v min . Consequently,  Figure 6. Same as the right panel of figure 5 but for a 3.5 GeV DM particle with exothermic scattering (δ = −50 keV) and a Ge-phobic SI interaction (f n /f p = −0.8) [32] (as in figure 3). Halo functions ηv min entirely contained within the plausibility region (light purple) lead to a compatibility of the data sets at the chosen level (p >10%). For those not entirely contained within the plausibility region the data sets are incompatible at the chosen level (p <10%).
any halo function not entirely contained within the plausibility region passes though points with p < 10 %, and thus for these functions the data sets are incompatible at the chosen level (p < 10%). However, halo functions that are entirely contained within a plausibility region are not necessarily allowed by our test, i.e. do not necessarily lead to compatibility of all data at the chosen level. The issue is that the true halo model adopted at each point within a plausibility region, namely theη c BF of the profile global likelihood at each point, may also pass through points outside the plausibility region and be rejected by our test. If so, the p-value evaluation at the particular point in the plausibility region is inconsistent. This is the case for all points in the plausibility regions (light purple) shown in figure 5. The regions are crossed by thin black lines to indicate that halo functions entirely contained within them are not guaranteed to lead to compatibility of the data sets. However, if the true halo model adopted at all points within a plausibility region are entirely contained within it, the p-value calculation is reliable and halo functions entirely contained with this region are allowed by our test. This is the case of the plausibility region in figure 6 (shown in light purple).
The plausibility region for p ≥ 10% arising from the constrained parameter goodness-offit test for the combination of CDMS-II-Si and SuperCDMSLT5 (light purple region) is compared in figure 5 to the confidence bands (red shaded regions) generated from the global likelihood described in section 4 (right panel) and the confidence bands generated with CDMS-II-Si data alone (left panel). Also shown are the global best fit halo functionη BF (dark red) and the best fit halo function for CDMS-II-Si alone (blue). The left panel of figure 5 shows that there does not exist a halo function in the CDMS-II-Si confidence bands that can describe the compatibility of the observed data sets. Not all halo functions contained within the 90% global confidence band in the right panel of figure 5 are excluded by the plausibility region, but the 68% region is entirely excluded, as is the global best fit halo function. As explained above the plausibility regions are crossed over because the functions entirely included within them are -17 -

JCAP10(2016)029
not allowed by our test (while those passing through points outside them are rejected). By contrast, we can see in figure 6 how the plausibility region includes the entire global confidence bands (as well as the bands for CDMS-II-Si alone, which in this case are nearly identical, see figure 3) in the case of non-conflicting data sets. This is the example of a 3.5 GeV DM particle with exothermic scattering (δ = −50 keV) and a Ge-phobic SI interaction (f n /f p = −0.8) [32]. The plausibility region provides in this case a further indication of compatibility of the CDMS-II-Si and SuperCDMSLT5 data sets for this particular DM particle model, besides the near complete overlap of the global and single experiment confidence bands.
A comment is in order regarding figure 6. While generating the probability distributions at large values ofη and v min (above the 90% CL band), we found the predicted number of events in both experiments was too large for our computational methods to work. We resorted to using a nearest-neighbor extrapolation at fixed v min to generate the probability distributions in this region. We found that in this region of the v min −η plane, the probability distribution changes slowly with respect to the observed value of q c P G , and thus we believe we obtained a good estimate of the upper boundary of the plausibility region. This extrapolation was only used above the 90% CL confidence band boundary.

Conclusions
In this paper we have presented two distinct methods to assess the joint compatibility of data sets for a given DM particle model across halo-independent parameter space, using a global likelihood consisting of at least one extended likelihood and an arbitrary number of Gaussian or Poisson likelihoods. We have illustrated these methods by applying them to CDMS-II-Si and SuperCDMS data, assuming WIMP candidates with SI contact interactions.
The first method is a natural extension of the procedure presented in [46], in which a best fit halo function and pointwise confidence band are constructed from the profile likelihood ratio. Here we have proven that the best fit halo functionη BF for the global likelihood we studied is a piecewise constant function with the number of steps at most equal to the number of unbinned data points plus the number of data bins in all the single likelihoods, and argued why in practice the number of steps is smaller than this maximum number (see section 3 and appendix A). A best fit piecewise constant halo function had already been found in the literature (see [36]) for a global likelihood of the type we use, but as a curiosity without any explanation (or proof of uniqueness). In addition to showing how to find the best fit halo functionη BF and that this function is unique (see appendix B), here we have shown for the first time how to construct two-sided confidence bands at any CL for the type of global likelihood we studied. As an illustration of the method we have found the best fit halo function and the 68% and 90% CL confidence bands assuming two different choices for the DM particle model parameters m, δ, and f n /f p . The choice of a 9 GeV DM particle scattering elastically (δ = 0) with an isospin-conserving coupling (f n /f p = 1) leads to an apparent incompatibility between the observed CDMS-II-Si events and the SuperCDMS upper limit, in agreement with previous published results (see e.g. [34,46]). This incompatibility can be assessed by comparing the overlap or lack thereof of the global confidence bands with those of CDMS-II-Si alone. As shown in figure 2, at the 68% CL, it is not possible to find a halo function passing through both confidence bands. The situation is very different for a 3.5 GeV DM particle with exothermic scattering (δ = −50 keV) and a Ge-phobic SI interaction (f n /f p = −0.8) [32], for which the data sets are compatible. As shown in figure 3 the global and CDMS-II-Si alone confidence bands practically coincide.

JCAP10(2016)029
The drawback of this method is that it cannot provide a quantitative measurement of the level of incompatibility of the various data sets that comprise the global likelihood. To address this concern, we have proposed in section 5 a second method in which we construct a "plausibility region" arising from the global likelihood, using an extension of the parameter goodness-of-fit test [36,48,49], that we refer to as the "constrained parameter goodness-of-fit" test. By evaluating the ratio of the global profile likelihood and the product of the individual profile likelihoods (assumingη(v * ) =η * ), a plausibility region can be constructed by grouping together regions of parameter space for which, at each point (v * ,η * ), our observed test statistic has a p-value e.g. ≥ 10%. This p-value was determined using a probability distribution constructed with Monte Carlo generated data assuming the true halo function is the constrained best fitη c BF of the profile global likelihood, i.e. the halo function that maximizes the global likelihood subject to the constraintη(v * ) =η * . For any halo function not entirely contained within this plausibility region the data are incompatible for the assumed DM particle model at the assumed level (e.g. p < 10%). For halo functions entirely contained within the plausibility region the data sets are compatible at the chosen level only if the contained best fit at each point within the region are also entirely contained within the region. We have demonstrated this method for a 9 GeV DM particle scattering elastically with an isospin conserving coupling and for the aforementioned Ge-phobic particle candidate. The results are shown in figures 5 and 6 respectively. In the first case the confidence bands are largely outside the plausibility region, while in the second case the confidence bands are entirely included in the plausibility region and any halo function entirely contained within the plausibility region lead to a compatibility of the data sets at the chosen level (p > 10%).
Together these two methods provide complementary assessments of the compatibility of the data given a particular dark matter model, across the v min −η halo-independent parameter space. We expect these tools to prove useful for future direct dark matter searches both to test compatibility of different data sets as to provide a guidance of which type of halo functions provide a better or worse compatibility of all the data.

A The zeros of the q(v min ) function
We are first going to argue that the zeros of q(v min ) are only isolated above a certain v min range where all terms in the sum defining q(v min ) are zero. We will then find the maximum possible number of isolated zeros, although the actual number of zeros can be much smaller than the maximum.

A.1 The zeros are isolated above a certain v min value
The terms defining q(v min ) in eq. (3.6) are either positive semidefinite, e.g. ξ EHI (v min ) and some of the terms proportional to ξ smooth discussion of the behavior of these functions, let us introduce the label µ, which will be used to denote either a quantity associated with the EHI experiment or an experiment-bin pair (α, j). This way, quantities like ξ µ (v min ) can either represent ξ EHI (v min ) or ξ (α) j (v min ). Each term in the sum in eq. (3.6) has a different v min -dependence. In particular the ξ µ (v min ) functions (note the general behavior of ξ (α) j (v min ) is identical to that of ξ(v min ) discussed in section 2.2) are zero below certain values of v min , which we will refer to as v µ low , strictly increase with v min (although the second derivative may exhibit sign changes), until at some value of v min , call it v µ high , they plateau and become constant. The v µ low and v µ high of each ξ µ (v min ) function, as well as the height of the plateau, depend on theoretical framework and the specifics of the experiments (e.g. the scattering kinematics, the differential cross section, the energy resolution functions, etc.). The H EHI a (v min ), also described in section 2.2, are instead upward step-like functions, starting from zero at low v min , with the steps appearing roughly at the v min values corresponding to the detected energy of the events observed in the EHI experiment.
In addition to having unique v min -dependencies, each of the terms in eq. (3.6) has uniquely definedη-dependent coefficient. Thus the terms are all independent of each other and have very different functional forms.
For values of v min below the minimum v µ low , i.e. where all the terms in eq. (3.6) are zero, q(v min ) is zero, which impliesη BF (v min ) is undetermined. This is not detrimental to the arguments we have made as it reflects the fact that experiments under consideration do not probe the halo function at these values of v min . Notice that in order to have non-negative q(v min ) values, the v low of some of the positive terms must be smaller than the smallest v low of all negative terms.
For values of v min larger than the minimum v µ low , zeros of q(v min ) can appear where the modulus of the sum of all negative terms in eq. (3.6) touches from below the sum of all positive terms in eq. (3.6) (recall that q(v min ) is a non-negative function). The positive terms consist of different ξ µ (v min ) (most of them multiplied byη-dependent coefficients). Thus, in general, the sum of all positive terms behaves as a monotonically increasing function starting from zero at the lowest v µ low (lowest of all positive terms) and plateauing to a constant value at the largest v µ high (again considering only positive terms). The negative terms in eq. (3.6) include the step-like H a (v min ) (multiplied byη dependent coefficients), which each add a "step-like" feature to the modulus of the sum of negative terms, and some of the ξ (α) j (v min ) dependent terms (multiplied byη dependent negative coefficients). Depending on the nature of these negative ξ (α) j (v min ) terms, they could add "shoulder-like" features, arising from changes in the sign of the second derivative, to the modulus of the sum of negative terms. The modulus of the sum of negative terms also plateaus above the largest v µ high (largest of all negative terms). The plateau of the sum of positive terms and the plateau of the modulus of the sum of all negative terms are entirely independent of each other, and thus the possibility that the two plateaux would coincide to produce q(v min ) = 0 is completely unrealistic since they both depend on entirely different experimental features. Furthermore, for most realistic cases, the maximum value of v µ high is larger than the galactic escape velocity, and thusη(v min ) should be zero in this region. Since these plateaus cannot feasibly coincide, q(v min ) cannot equal 0 above the largest v µ high . Typically isolated zeros of q(v min ) would happen when some of the "step-like" or "shoulder-like" features of the modulus of the sum of negative terms in eq. (3.6) touch from below the monotonically increasing sum of all positive terms in eq. (3.6). Alternatively, if the -20 -

JCAP10(2016)029
sum of the positive terms has a region of negative curvature, it may be possible for that this sum could reach towards and touch the modulus of the sum of negative terms from above.
A practically impossible conspiracy between terms dependent on different experiments, energy intervals, andη functions would be required for q(v min ) to be zero in an extended v min interval above the minimum v µ low , a conspiracy which would not survive infinitesimal changes in any of the elements defining each term in eq. (3.6). We include in appendix A.3 a more mathematically rigorous proof illustrating why extended zeros of q(v min ) cannot exist above the minimum v µ low . In the following we only consider the possibility that q(v min ) has a finite number of isolated zeros.

A.2 Maximum number of isolated zeros of the function q(v min ) for a global likelihood
Before counting the number of isolated zeros of q(v min ), let us introduce the notion of a "generic" solution. We say that a solution is generic if small changes in the quantities that define it do not affect the existence of the solution. In our context, the quantities defining the solutions are the input parameters and functions given to fully specify ξ EHI , ξ (α) j , and H EHI a , e.g. the efficiency function (E , E R ), the energy resolution function G T (E , E R ), the differential cross section dσ T /dE R , and the exposure M T for each experiment and bin.
Let us briefly demonstrate the importance of the concept of generic solutions by considering the number of isolated zeros that can arise in the linear combination of two functions f (x) and g(x) which do not have the same functional form, since they are assumed to be derived from two independent experimental setups (i.e. changes in the experimental quantities of one experiment may affect e.g. f (x), but do not affect g(x) in the same manner). For an adjustable parameter λ, it is possible for f (x) and λg(x) to have a generic point of osculation, i.e. a point where f (x) = λg(x) and f (x) = λg (x), at which the Wronskian W [f, g] vanishes In fact, W [f, g] could vanish in more than one point, say x 1 , x 2 , . . . x n , or in various intervals. In this case the value of λ can be chosen so that f (x 1 ) = λg(x 1 ) at one of those discrete points, say x 1 . This point of osculation defines an isolated zero of the function [f (x) − λg(x)], with zero slope. Having two points of osculation, say x 1 and x 2 , would require  except here we will treat the λ m as free parameters. The argument above ensures that there could be at most one generic point of osculation between ξ EHI (v min ) and λ m X m (v min ), or between λ m X m (v min ) and λ k X k (v min ) for k = m. Here, the coefficients λ m are adjustable parameters, equivalent to a multidimensional generalization of the parameter λ in the above example. In the context of eq. (3.6), one can identify the λ m with the halo-dependent quantities, e.g. 1/γ[η] and the factors in the square bracket of eq. (3.7) and eq. (3.8).

JCAP10(2016)029
and v (k ) (λ 1 , . . . , λ k−1 , λ k+1 , . . . , λ k −1 , λ k +1 , . . . , λ n ) ≡ v (k ) (λ 1 , . . . , λ k −1 , λ k +1 , . . . , λ n )| λ k =λ k , (A. 12) which are respectively induced from the functions defined on M In a similar way, if all coefficients λ k could be freely adjusted the intersection of all (n − 1)-dimensional sub-manifolds, is generically a zero-dimensional sub-manifold of M n , i.e. a set of discrete points. For one of these points, which we call (λ 1 , . . . ,λ n ), we can define the function which has n isolated zeros, with zero slope. Here, n = N ≡ N EHI + α N (α) bin (see eq. (3.5)), i.e. the number of events observed by the EHI experiment plus the total number of bins employed by all Poisson and Gaussian experiments. This is what we wanted to prove. However we have so far assumed the coefficients λ m could all be freely adjusted. This is not true, however, and the actual number of isolated zeros of q(v min ) (with q (v min ) = 0) will be in most circumstances much smaller than the maximum N .
In fact, the coefficients λ m are quantities derived from a halo functionη. All points in M n that can be actually realized from halo functionsη form a continuous subset S of the manifold M n . The maximum number of the individual sub-manifolds M k n−1 passing through a point in S gives the maximum number of actual possible steps in the best fitη function. This number can be determined by carefully considering the functional form of ξ EHI , ξ (α) j , and H EHI a , and is in general smaller than N .
A.3 Argument against non-isolated zeros of q(v min ) Here, we provide a more mathematically rigorous proof for why q(v min ) cannot have nonisolated zeros above the minimum v µ low . Using eq. (2.21) we can equate the functional derivative of L[η] to the derivative of q(v min ) as δL δη(v min ) = ∂ ∂v min q(v min ) . (A.15) We will begin by assuming that there exists some interval [v 1 , v 2 ] above the minimum v µ low in which q(v min ) = 0, and prove by contradiction that this cannot be the case. Let us introduce an infinitesimal perturbation δF (v) in the speed distribution F (v), that is only non-zero in the interval [v 1 , v 2 ], and define the quantity

JCAP10(2016)029
Since q(v) = 0 in the interval [v 1 , v 2 ] and δF (v) = 0 outside of this interval, ∆ is trivially zero. Since the halo functionη linearly depends on F (v), the induced change in the halo functionη is and its derivative is given by Performing integration-by-parts on the integral in eq. (A. 16) gives where the last line is obtained from the definition of the functional derivative. This expression for ∆ implies that any change of F (v) introduced above the minimum v µ low in an interval where q(v min ) = 0 would not change the value of the likelihood. However, by definition, this perturbation necessarily introduces a constant shift inη at all v min values below v 1 . It is inconceivable for a likelihood function to be invariant under such a rigid shift ofη. Therefore, our original assumption must have been false, and there cannot exist intervals above the minimum v µ low in which q(v min ) = 0. The interested reader might wonder why the proof presented above does not preclude the existence of isolated zeros. In this case, the modification to the speed distribution δF (v) must take the form of a constant times a delta function (i.e. nonzero only at the location of the isolated zero). The resultant speed distribution F (v) + δF (v) would no longer be a smooth function, and thus cannot possibly be representative of a true speed distribution (no physical process would produce a component with zero dispersion).

B The uniqueness of the best-fit halo function
Here, we show that the halo functionη BF (v min ) maximizing a global likelihood functional having at least one extended likelihood as a factor is unique in the v min range wherein the experiments in consideration can probe the value of the halo function. The proof consists in showing that the second directional derivatives of the functional L ≡ −2 ln L with respect to variations of theη function are all positive.

B.1 Statement of the proof
We start by stating two properties of the global likelihood (and the individual likelihoods) considered in this paper. First, the likelihood depends on the halo function only through physically observable quantities, which are either the scattering rate in a bin

JCAP10(2016)029
Expanding the functional L around the best fit halo up to the second order we get (B.11) Here we defined the index A to run over all data points, namely all bins or single events of all experiments considered, so A runs from 1 to N . Specifically, the summation over A runs over the observed events in the extended likelihoods and the bins j of all experiments α with Poisson and Gaussian likelihoods. The index i indicates instead each constant portion of the best fit halo function, between the steps at v min values v i−1 and v i . The maximum number of steps was found in appendix A to be N , so we can take take the number of steps to be equal to N and consider some of the step heights to be zero. In this way, i also runs from 1 to N . The coefficients K Ai are given by for extended likelihoods, for Poisson likelihoods, and (B.14) for Gaussian likelihoods. In the last two equations the index A accounts for the experimentbin pairs indexes (α, j). Notice that the quantities K Ai can be interpreted as the components of N non-zero vectors K A in a vector space with dimension N with components denoted by i. We can also consider ∆η i to be the components of a vector ∆ η with the same number of dimensions of the K A vectors. Each vector ∆ η is a possible infinitesimal variation of the heights of the stepsη i around a best fit halo function. Eq. (B.9) is then a sum of the squares of the inner products of two vectors K A and ∆ η.
Notice also that the vectors K A are generically linearly independent, because there is no reason that the experiment-specific quantities K Ai should be dependent upon information contained in a different bin or experiment.
Since the vectors K A are generically linearly independent there is no non-zero vector ∆ η orthogonal to all of them. This implies that there is no infinitesimal variation of the heights of the stepsη i around a best fit halo function for which the second order variation of L vanishes. This proves that the likelihood functional L is not invariant under any infinitesimal variation around the best fit halo function which would lead to another best fit halo function. Since all the second directional derivatives of L around a best fit halo function are positive the best fit halo function must be unique.