Studying generalised dark matter interactions with extended halo-independent methods

The interpretation of dark matter direct detection experiments is complicated by the fact that neither the astrophysical distribution of dark matter nor the properties of its particle physics interactions with nuclei are known in detail. To address both of these issues in a very general way we develop a new framework that combines the full formalism of non-relativistic effective interactions with state-of-the-art halo-independent methods. This approach makes it possible to analyse direct detection experiments for arbitrary dark matter interactions and quantify the goodness-of-fit independent of astrophysical uncertainties. We employ this method in order to demonstrate that the degeneracy between astrophysical uncertainties and particle physics unknowns is not complete. Certain models can be distinguished in a halo-independent way using a single ton-scale experiment based on liquid xenon, while other models are indistinguishable with a single experiment but can be separated using combined information from several target elements.


Introduction
Over the past decade dark matter (DM) direct detection experiments have improved their sensitivity by an order of magnitude every two years and this trend is expected to continue for the near future [1][2][3][4][5][6][7][8][9]. While there is at present no clear evidence for the interactions of DM particles in nuclear recoil detectors [10][11][12][13][14], it is perfectly conceivable (and in fact predicted by many models for DM) that hundreds of events will be observed by 2020. Once a signal is seen in one or several direct detection experiments, the challenge will be to identify those models of DM that allow for a good fit of the experimental data and to determine the preferred values of the underlying parameters, such as the mass of the DM particle.
At the same time, however, it has become clear that the interactions of DM particles with nuclei can be significantly more complicated than suggested by the simple division into spin-independent and spin-dependent interactions. In general the scattering cross section can also depend on the momentum transfer and the relative velocity between the DM particle and the nucleus. A full classification of all possible DM interactions in the non-relativistic limit requires no less than 28 different scattering operators as well as a large number of nuclear response functions [49,50], which can significantly affect the expected signals in direct detection experiments [46,[51][52][53][54][55][56][57][58][59][60]. The possible importance of the additional DM-nucleon scattering operators makes the issue of astrophysical uncertainties even more pressing, because a nonstandard velocity distribution can potentially mimic non-standard DM interactions.
Fortunately, the degeneracy between astrophysical uncertainties and particle physics unknowns is not complete. For example, standard spin-independent interactions induce a differential event rate that decreases monotonically with increasing recoil energy for any DM velocity distribution [33]. The observation of a non-monotonic differential event rate could hence not be attributed to astrophysical uncertainties and would instead point strongly towards non-standard interactions. It should therefore always be possible to obtain at least some basic information on the nature of DM interactions even when accounting for astrophysical uncertainties.
In the present paper we develop a framework to combine the full formalism of nonrelativistic effective interactions with state-of-the-art halo-independent methods. Our approach is based on the idea of decomposing the velocity distribution into a finite sum of streams with velocity v j and then varying the normalisation of each stream [30,31]. We will show that, even for non-standard interactions, it is always possible to calculate a matrix D ij such that the number of events N i in the ith bin of a given experiment can be calculated via simple matrix multiplication N i = D ij g j with the g j being determined by the normalizations of the streams. This simple relation allows to determine the velocity distribution that best describes a given set of data for assumed particle physics properties.
We can repeat this procedure for different particle physics assumptions in order to study whether changes in the velocity distribution can compensate for changes in the particle properties of DM and thus reduce our ability to determine these properties unambiguously. The aim is to quantify what information can be inferred on the coupling structure of DM in a halo-independent way. For example, if no good fit to a given set of data can be found for any DM velocity distribution, the corresponding particle physics assumptions can be disfavoured without the need to make assumptions on the astrophysical distribution of DM. This approach makes it possible to analyse direct detection experiments for arbitrary DM interactions, in particular DM interactions with non-standard momentum and velocity dependence, independent of astrophysical uncertainties.
To illustrate the general formalism we study a representative set of DM-nucleon interactions. These consist of the standard spin-independent and spin-dependent scattering scenarios, interactions induced by an anapole or a magnetic dipole moment of DM, as well as a dipole interaction involving a new heavy mediator. This set of interactions, though certainly not exhaustive, covers all the central aspects relevant for a halo-independent investigation of non-standard interactions between DM and nucleons. We discuss the scattering rates induced by these models in the context of three future direct detection experiments based respectively on xenon, germanium and iodine targets. We first study whether a single (xenon-based) experiment can distinguish the different models of DM in a halo-independent way and then focus on the question whether the complementarity of several different target materials can improve the distinction.
A similar analysis of the possibility to distinguish different DM models using future direct detection experiments has been performed in [58]. While considering a larger set of DM interactions and more different experiments, the analysis makes much more specific assumptions on the velocity distribution of DM. Our results extend the findings from [58] to arbitrary DM velocity distributions (and furthermore take into account effects from finite energy resolution).
This paper is structured as follows. In section 2 we review the basic formalism for direct detection, including the central ideas underlying the halo-independent methods for DM-nucleon interactions with standard velocity dependence. We then explain how to generalise such a framework to more complicated DM interactions and derive the central formulas used to calculate the matrix D ij . In section 3, after defining the set of interaction scenarios discussed in this work, we provide a qualitative illustration of the halo-independent interpretation of direct detection data in the context of these models. We then introduce methods for quantifying the goodness of halo-independent fits to experimental data. In sections 4 and 5 we apply the method to various combinations of direct detection experiments, and discuss in particular which of the above-mentioned interaction scenarios can be distinguished in a halo-independent way. Section 6 provides a summary of our findings. Additional details on the parametrization of non-relativistic effective interactions and the calculation of event rates are given in the appendices.

Extended halo-independent methods
This section describes how the integral of the DM velocity distribution can be used to study astrophysical uncertainties in direct detection experiments. After briefly recalling the formalism of writing the event rates in terms of velocity integrals, we review the approach for analysing direct detection data in a halo-independent way for the case where the differential cross section has the standard velocity dependence dσ/dE R ∝ v −2 . Crucially for the remainder of this work, we then generalise the method to more complicated cross sections. Additional technical details related to this section can be found in the appendices.

Calculating event rates from velocity integrals
The differential event rate with respect to recoil energy in a direct detection experiment is given by where ρ is the local DM density, m χ and m T are the DM and target nucleus mass, f (v) is the local DM velocity distribution evaluated in the Galactic rest frame, v E (t) is the velocity of the Earth relative to the Galactic rest frame and v = |v|. For elastic scattering the minimum velocity required for a DM particle to transfer the energy E R to a given nucleus T is is the reduced mass of the DM-nucleus system. The differential cross section dσ/dE R encodes the details of the interactions between DM particles and nuclei. As discussed in detail in Appendix A, for a very general set of effective operators, the differential cross section can always be decomposed into a contribution proportional to v −2 and a velocity-independent contribution: Defining the velocity integrals 1 we can then simply write the differential event rate as If the second term of this equation is absent, all information on the DM velocity distribution is encoded in the velocity integral g(v min ). As a result, this quantity is particularly convenient for studying the impact of astrophysical uncertainties on observables in direct detection experiments, such as the number of events R i in a bin of the form [E i , E i+1 ]. In order to study how R i depends on the velocity integral, it was proposed in [30] to parametrize g(v min ) as a piecewise constant function with N s steps. 2 By choosing N s large enough, this allows for an arbitrarily good approximation to any possible velocity integral. Concretely, as suggested in [30], we divide the region of v min -space probed by the experiments under consideration into N s intervals of the form [v j , v j+1 ], where the number of steps N s can be as large as 50. In each interval, the velocity integral is taken to be constant: (2.6) The fact that the velocity integral is both positive and monotonically decreasing leads to the requirement 0 ≤ g j ≤ g j−1 for all j. It was then shown in [30] that it is possible to encode all experimental details and particle physics properties of DM into a matrix G ij such that R i = j G ij g j . Given the very simple relation between the parameters g j and the observables R i it is then straight-forward to find the velocity integral that leads to the best agreement with data, effectively profiling out all halo uncertainties. 1 Here we make two approximations. First, we neglect the small time dependence of vE and second we assume that dσ/dER only depends on v 2 and not on the direction of v (e.g. because the target is unpolarised and the direction of recoil tracks cannot be distinguished).
2 Note that this is equivalent to writing the speed distribution f (v) as a sum of δ-functions. In other words, we decompose the velocity distribution into a sum of streams with fixed velocities and negligible velocity dispersion.

Halo-independent methods for generalised dark matter interactions
At first sight, the situation becomes more complicated if both g(v min ) and h(v min ) contribute to the differential event rate. The crucial observation is however that g(v min ) and h(v min ) are not independent. In fact, we can make use of the fact that g [39,47] In other words, it is possible to calculate h(v min ) directly from the usual velocity integral g(v min ). This observation is key to generalising existing halo-independent methods to models with a more complicated velocity dependence. To this end, we write the velocity integral as where Θ(x) is the Heaviside step function and l j ≡ g j − g j+1 , with the g j defined via equation (2.6), setting in addition g Ns+1 ≡ 0. We can then use equation (2.7) to calculate h(v min ): The parameters g j and h j are then related by a simple linear transformation: with the matrix F given by F jj = v 2 j+1 , F jj = 0 for j > j and F jj = v 2 j +1 − v 2 j for j > j (see also appendix B.1).
In order to calculate predicted event rates for a given direct detection experiment we need to know the probability p(E D ; E R ) that a DM interaction with true (i.e. physical) recoil energy E R leads to a signal with detected energy E D . In terms of this function, and using our ansatz for g(v min ), we find where E j is defined implicitly by v j = v min (E j ). If we are interested in the number of events R i predicted in a given bin, we simply need to multiply dR/dE D with the total exposure κ and integrate over E D : We now make the simplifying assumption that the detection efficiency depends only on the physical recoil energy E R (and not on E D ) and that fluctuations can be approximated by a Gaussian distribution with standard deviation given by ∆E(E R ). In this case we can write The integration over E D can then be performed analytically, and one obtains where with erf(x) = (2/ √ π) x 0 exp(−t 2 ) dt. We can further simplify equation (2.14) by making use of the relation h k = j F kj g j derived above. It is then possible to combine the two matrices G ij and H ij into one matrix The final expression is thus formally identical to the one used in [30]. In other words, it is possible to obtain the same concise expression relating the binned event rates and the discretized velocity integral even if the differential cross section has a more complicated velocity dependence. The matrix D ij encodes all relevant experimental details and the assumed DM properties, but it is independent of the DM velocity distribution. Once the matrix D ij has been calculated, it is possible to calculate the predicted number of events for each bin for any given velocity integral parametrized by the coefficients g j . These predictions can be used to calculate the likelihood for a specific choice of the g j . Denoting the expected number of background events in each bin by B i and the total number of observed events by N i , this likelihood is given by Since all entries of D ij are positive, one can show that this likelihood is a convex function of the parameters g j defined on a convex domain and therefore any local minimum is guaranteed to be a global minimum. In other words, in spite of the large number of steps N s it is possible to unambiguously minimise the likelihood with respect to the DM velocity distribution and thereby find the astrophysical parameters that give the best fit to a given set of observations. 3 The equations above generalise directly to the case of several experiments (labelled by α = 1, . . . , N exp ). In this case, the matrix D (α) ij must be calculated separately for each experiment. 4 One can then calculate the likelihood L (α) for each experiment given the expected number of background events B i . The total likelihood for given g j is then calculated by multiplying together the likelihoods for the individual experiments.
Additional technical details on how to calculate the matrices G ij and H ij are provided in the appendices. Appendix A.4 presents the method we use in order to implement equations (2.15)-(2.17) in a numerically efficient way. In appendix B.2 we discuss the problem of determining the appropriate ranges of v min -space for DM interactions with non-standard velocity dependence.

Studying generalised DM interactions
In the remainder of this paper we illustrate the possible applications of the framework developed above by addressing the following question: Given a set of experimental data, what information can be inferred on the coupling structure of DM without making any assumptions on its astrophysical distribution? To answer this question, we will generate sets of mock data for various DM models and then determine the goodness-of-fit for a number of different coupling structures. The resulting information can be used both to understand qualitatively whether certain coupling scenarios are distinguishable in principle as well as to quantify how strongly a specific coupling scenario is favoured or disfavoured by a given set of data.
In this section, we present the models that we consider for this study and illustrate the importance of halo-independent methods. Furthermore, we explain the quantitative method used to compare the different coupling structures. In section 4 we then focus on the simplest case, where only a single experiment is sensitive to the interactions of DM. More complicated scenarios with more than one experiment are discussed in section 5.

Models
We apply the halo-independent approach introduced in section 2 to five different benchmark models of DM interactions: standard spin-independent and spin-dependent scattering, DM featuring an anapole or a magnetic dipole moment, as well as DM interacting via a "dark magnetic dipole moment". Our main reasoning for choosing these scenarios, which will be discussed in more detail below, is that they encompass the most relevant nuclear response functions and shapes of recoil spectra, and thus should constitute a representative set for describing the phenomenology of possible future direct detection data. Moreover, each of our benchmark scenarios can arise in well-motivated UV complete models of DM, making our choice also motivated from the model building perspective. We remark that our set of scenarios is similar to the one discussed in [53,58].
In the following, we define each of the benchmark models in terms of an effective DMnucleon Lagrangian, expressed via the non-relativistic operators O τ λ [49,57] (see also equation (A.1)). For the generation of mock data for future experiments, we assume in all scenarios a true DM mass m DM = 50 GeV, as well as a standard Maxwell-Boltzmann velocity distribution, with the most probable speed v 0 = 220 km/s, the Galactic escape velocity v 0 = 544 km/s and the mean velocity of the Earth v E = 230 km/s. Furthermore, we fix the normalization of the interaction strength (e.g. the total scattering cross section in the standard SI/SD case) such that it complies with all current direct detection bounds. Concretely, we consider the experiments LUX [10], SuperCDMS [13], SIMPLE [61], COUPP [62], PICASSO [63] and PICO-2L [12], all of them implemented as in [60], and choose the normalization of the interaction rate in each benchmark model to be a factor of two below the most sensitive exclusion bound.
Standard spin-independent and spin-dependent scattering (SI/SD) A large variety of DM models lead to the standard spin-independent and/or spin-dependent interaction, which are described by the (relativistic) effective DM-nucleon Lagrangians f (SI) Nχ χN N and f (SD) Nχ γ µ γ 5 χN γ µ γ 5 N , respectively. In the notation of [57], the resulting non-relativistic Lagrangians read respectively. 5 Concretely, for the spin-independent case, we will consider a benchmark scenario with fixed neutron-to-proton coupling ratio f n /f p = 1, corresponding e.g. to Higgs exchange, with an absolute strength given by σ Note that while we restrict ourselves to specific values of f n /f p for the purpose of generating the mock data, we will in some cases allow arbitrary neutron-to-proton coupling ratios for the halo-independent fits that we perform.

Anapole moment (AM)
Next, we consider DM interacting dominantly via an anapole moment, corresponding to the DM-photon interaction Aχγ µ γ 5 χ ∂ ν F µν , which constitutes the only possible electromagnetic moment of a Majorana fermion; see e.g. [64] for a concrete realization. The coupling of the anapole moment to the nucleon electromagnetic current then leads to [53] where e is the electric charge, Q p = 1 and Q n = 0 is the charge of the proton and neutron in units of |e|, respectively, and µ N is the magnetic moment of the nucleon N in units of the nuclear magneton e/(2m p ). This interaction structure gives rise to a scattering cross section with a non-standard dependence on v 2 and q 2 : it contains a spin-independent part with an additional factor of v 2 ⊥ ≡ v 2 − q 2 /(4µ 2 T ), as well as a contribution rising as q 2 , which is enhanced for nuclei with a significant magnetic dipole moment. For the generation of the mock data, we again assume m DM = 50 GeV, as well as an anapole moment A = 3.5 · 10 −5 GeV −2 , compatible with all current direct detection bounds.

Magnetic dipole moment (MDM)
Dirac DM can interact with photons via a magnetic dipole moment given by D magnχ σ µν χF µν (there is also an equivalent operator for complex scalar DM). Scattering induced by the magnetic dipole moment of DM is known to be the leading contribution to the total interaction rate e.g. in scenarios where a Dirac DM particle couples to a Standard model lepton and a heavy scalar mediator via a Yukawa interaction [64,65]. The effective non-relativistic DM-nucleon Lagrangian in this scenario is given by [53] Compared to the standard spin-independent or spin-dependent scattering, this interaction structure leads to a differential scattering cross section which is enhanced by 1/E R at small recoil energies, implying a steeper recoil spectrum. For our benchmark scenario, we choose D magn = 9.3 · 10 −7 e fm. Similar to the magnetic dipole moment, DM could also interact with nucleons via an electric dipole moment (which however would require significant CP violation in the dark sector) [66,67]. We have found that for the purpose of this work, magnetic and electric dipole DM are essentially equivalent from a phenomenological perspective, and hence for simplicity we only include the former scenario in our list of benchmark models.

Dark magnetic dipole moment (DMDM)
Lastly, we consider the case of a "dark magnetic dipole moment", which is obtained by replacing 1/q 2 → 1/M 2 in the MDM scenario given by equation (3.3) [68]. Physically, this corresponds to a situation in which DM interacts via a dipole operator with a new heavy gauge boson of mass M , as opposed to the MDM scenario involving the exchange of a photon. Redefining the overall normalization of the interaction strength as D DMDM (which has units GeV −3 ), we obtain Phenomenologically, this model has the interesting feature that it leads to a recoil spectrum rising like q 2 at small momentum transfers, strikingly different from the standard spinindependent or spin-dependent scenarios. For the generation of the mock data corresponding to this benchmark model, we assume D DMDM = 1.6 · 10 −4 GeV −3 .

Mapping to v min -space
For a given velocity integral the models introduced above make very different predictions for the number of expected events in different experiments as well as for the shape of the individual recoil spectra. Thus, one would expect that they can be easily distinguished once a sizeable number of DM scattering events have been observed. Such a reconstruction of the DM properties, however, is often only possible if a specific form of the velocity integral is assumed. Unfortunately, incorrect assumptions on the velocity integral may strongly bias the reconstruction and may even lead to the exclusion of the correct model of DM [22,23,30]. The halo-independent method introduced in section 2 does not require such assumptions and therefore avoids false exclusions. A direct consequence is that it becomes harder in such an approach to distinguish different models of DM. Whenever changes in the velocity distribution can compensate for changes in the particle physics properties of DM, it is conceivable that a good fit to the data can be obtained even for an incorrect model of DM. It is therefore of great importance to understand what experimental set-up is necessary to ensure that different models of DM can be distinguished even when allowing for arbitrary velocity distributions.
A convenient way to illustrate these issues is to map experimental data into v minspace [34,36,37]. Here we briefly review the general idea of this mapping and then use it to discuss a few illustrative examples for the importance of halo-independent methods. For this purpose let us consider an experiment that has observed a number of DM scattering events, allowing it to infer the differential event rate dR/dE R at a certain recoil energy E 0 . Under the assumption of a specific DM model, this measurement can be used to infer the value of the velocity integral g where A is the mass number of the nucleus, F (E R ) is the standard SI form factor and σ n is the DM-neutron scattering cross section at zero momentum transfer. Substituting this expression into equation (2.5), one finds To simplify this expression, we absorb all quantities that do not depend on the target material into the rescaled velocity integral The rescaled velocity integral can now be calculated from the measured differential event rate by inverting eq. (3.6):g This translation can be performed for further measurements of the differential event rate at different recoil energies and for measurements made by additional experiments. 6 Since all  Figure 1: Reconstruction of the velocity integral in comparison with the underlying mock data mapped to v min -space. In the top row the same model (SI interactions) has been used to generate the data and to reconstruct the velocity integral, while in the bottom row an incorrect model (MDM interactions) has been assumed for the reconstruction. In the right column Poisson fluctuations have been applied to the mock data to illustrate their effect on the best-fit velocity integral. Note that the error bars on the data points are shown for illustration only and are not used directly for the calculation of the best-fit velocity integral, which is determined from the binned event rates as described in section 2.
experiments probe the same velocity integral, it should be possible to fit the inferred values forg(v min ) with a single monotonically decreasing function. If no good fit can be found, the assumed model of DM can be ruled out in a halo-independent way. In figure 1 we illustrate this mapping using mock data from a single experiment based on liquid xenon (Xe), see section 4 for details. The top-left panel shows the simplest case, where we use the same model, namely SI interactions with m χ = 50 GeV, to generate the mock data and to perform the mapping to v min -space. In the absence of Poisson fluctuations, one can see that the inferred values ofg(v min ) agree perfectly with the SHM used to generate the mock data (indicated by the red dashed line). 7 The best-fit velocity integral for N s = 30 is indicated by the purple solid line. As expected, it agrees well with the SHM. In the top-right panel, we have included Poisson fluctuations in the data. As a result, the best-fit velocity integral now differs from the SHM, because it attempts to follow the random fluctuations in the data. The important observation is that even the best-fit velocity integral cannot necessarily achieve a perfect fit to the data, because it is required to be monotonically decreasing and can therefore not always follow upward fluctuations at high energies.
In the bottom row we consider a different situation and use a different model to generate the mock data, namely MDM interactions with m χ = 500 GeV, but still use SI interactions with m χ = 50 GeV for the mapping to v min -space. The different shape of the recoil spectrum for MDM interactions implies that the SHM no longer gives a good fit to the data. Nevertheless, by considering arbitrary velocity integrals it is still possible to obtain a good fit. In other words, it is not possible in this example to rule out SI interactions in a halo-independent way even if DM scattering actually arises from MDM interactions.
Clearly, this situation changes significantly if two experiments are considered (see section 5 for details), because the two experiments will provide independent measurements of the rescaled velocity integral at the same values of v min . This is illustrated in figure 2. In the top row we again consider the case that the true model is identical to the fitted model. As expected, the measurements from both experiments are in good agreement with each other and a consistent fit can be found, even though the fit will be less than perfect in the presence of Poisson fluctuations. Using MDM interactions to generate the mock data, we however find strong incompatibility between the inferred values ofg(v min ) from the two experiments. The reason is essentially that, compared to xenon targets, iodine targets are more sensitive to MDM interactions than to SI interactions.
Obviously, it is difficult to quantify these statements using the information from figure 2 alone. For example, once Poisson fluctuations are included it becomes less trivial to distinguish the top row from the bottom row. In particular, it is difficult from these plots to appreciate how sensitive experimental predictions are to variations in the velocity integral. A discussion of how to construct confidence bands in v min -space can be found in [44]. Here we take a different approach and focus entirely on the likelihood constructed from the binned event rates in energy space, which is given in eq. (2.19). We now discuss how one can use this approach to make rigorous quantitative statements about the incompatibility of different measurements and to conclude whether or not the interactions under consideration can be distinguished.

Goodness-of-fit estimates
For each of the models discussed above, and given a value of the DM mass m χ and the DM velocity distribution f (v), we can calculate the predicted number of events in a given experiment. Once we have calculated the predicted number of events in each bin, we can apply Poisson fluctuations to generate mock data. Each of the models from above can then be fitted to the mock data in order to quantify the likelihood that the data corresponds to this particular model, using the definition of L given in equation (2.19). If the fitted model is identical to the true model (assumed for generating the mock data), we generally expect to find a good fit, i.e. a likelihood within the range expected by random fluctuations in the data. If the fitted model is different from the true model, the likelihood may be much smaller, depending on how well true model and fitted model can be distinguished by the experimental set-up.
For a given set of mock data and a chosen model, we can determine the value of the DM mass that maximises the likelihood (i.e. that minimises −2 log L). Ideally, the resulting minimum value of −2 log L (called x 0 ) could simply be compared to a χ 2 distribution with the appropriate number of degrees of freedom in order to quantify the goodness-of-fit of the best-fit point. Unfortunately, as discussed in [30], this simple approach does not work in our context. Although the distribution of x 0 does approximately follow a χ 2 distribution (when considering random fluctuations in the underlying data), we cannot easily determine the appropriate number of degrees of freedom. The reason is that finding the optimum velocity integral typically requires a very large number of fitted parameters, many of which simply saturate the boundary condition from the monotonicity requirement without actually improving the fit to the data. In other words, even if the number of fitted parameters significantly exceeds the number of bins, we typically do not obtain a perfect fit, although the naive number of degrees of freedom should be zero.
While it is thus not possible to determine analytically the expected distribution for x 0 , this distribution can be extracted from Monte Carlo simulations [31]. For this purpose, we can simply take the best-fit point of the fitted model as the basis for generating new sets of mock data and then repeat the same fitting procedure as above. If the typical values for −2 log L determined in this second step (called x 1 ) are comparable to the one found in the first step above, the conclusion would be that the assumed model gives a good fit to the data.  If, on the other hand, the first step gave a much larger value for x 0 than what is typically found when taking the best-fit point to generate new sets of mock data, the goodness-of-fit is bad. More precisely, if ζ(x 1 ) is the (normalised) distribution of x 1 obtained from Monte Carlo simulations of the best-fit point, we can define the p-value for the assumed model as The procedure outlined above is straight-forward to carry out in the case of a single set of experimental data. Lacking a conclusive DM signal, however, we are interested in the typical performance of future experiments. Such a study requires a large number of mock data sets. We therefore adopt the following procedure: 1. For a given choice of the true model (used to generate the mock data) and of the fitted model (used to fit the mock data), we generate mock data sets and determine the bestfit DM mass, as well as the corresponding value of x 0 and the predicted number of events in each bin for the fitted model.
2. Once we have determined the best-fit models for a large number of mock data sets, we find a typical prediction of the fitted model by taking the median of the predicted number of events in each bin.
3. We then take these typical predictions to generate more mock data sets and once again repeat the fitting procedure from step 1 to determine x 1 .
For each combination of true model and fitted model, we end up with two distributions: the distribution of x 0 , corresponding to mock data generated from the true model, and the distribution of x 1 , corresponding to mock data generated from the typical predictions of the fitted model. If the two distribution are very similar, the conclusion is that true model and fitted model cannot be easily distinguished with the given experimental set-up and the chosen test statistic. 8 If the two distributions are different, we can use the results to address a number of interesting questions. First, what is the p-value corresponding to a typical value for x 0 (taking for example the median of x 0 across all mock data sets)? And second, in what fraction of the simulations is the p-value smaller than 5%, such that we would be able to exclude the fitted model at 95% confidence level? We shall refer to the former number as similarity (S) and the latter number as the distinguishability (D) of the two models. The two quantities are illustrated in figure 3. For completely indistinguishable models (in particular the case that true model and fitted model are identical), one would obtain S ≈ 50% and D ≈ 5%. Conversely, if S < 5% or D > 50% a typical realisation of the true model would allow to exclude the fitted model with at least 95% confidence level (C.L.).
Clearly it is an approximation to calculate the distribution ζ(x 1 ) only for a typical prediction of the fitted model (defined in step 2 of the procedure presented above), rather than for each set of mock data separately. In principle, one should perform a separate Monte Carlo simulation for each best-fit point in order to determine the p-value of each fitted model. One would then obtain a distribution of p-values, which can be used to determine the similarity (i.e. median p-value of the fitted model) and the distinguishability (i.e. the fraction of fitted models with a p-value smaller than 5%).
While it would be unfeasible to perform a separate Monte Carlo simulation for each realisation of the mock data and each combination of true model and fitted model, we have performed one such study for the case that the true model is DMDM and the fitted model is SI, employing mock data for a xenon-based experiment (see the following section for more details). The resulting distribution of p-values is shown in figure 4. The median p-value is found to be 0.079%, which agrees well with the value S = 0.1% inferred from the procedure defined above within the uncertainties resulting from limited Monte Carlo statistics. Conversely, the fraction of simulations with a p-value below 5% is 95.5%, compared to the value D = 94% obtained via the approximate method. The conclusion is that the approximation to use a typical prediction rather than each individual prediction is sufficient, because the spread of the predictions from the best-fit points is not too large.  Table 1: Specification of the assumed performances of the future direct detection experiments employed in this work. We assume uniform acceptance within the energy ranges given in the third column of the table.

Constraining DM coupling structures with a single experiment
In this section we consider a single experiment, namely a liquid-xenon-based ton-scale detector similar to XENON1T [4]. Details on the assumed exposure, energy window and detector resolution are given in table 1. We assume the contribution from backgrounds to be negligible. 9 Given current constraints on the interactions of DM, we can expect up to 200 DM scattering events in such an experiment. We will see that even in this minimal case there is some halo-independent sensitivity to the coupling structure of DM.
The results of the analysis discussed above are shown in figure 5. The different rows correspond to different choices of the true model, while the different columns correspond to different fitted models. In each case, the blue line shows the distribution of x 0 , while the red line shows the distribution of x 1 . We also show the resulting values for S and D in each panel.

Discussion
We can make a number of interesting observations from figure 5. First of all, we find SI interactions and SD interactions to be completely indistinguishable in our set-up. This result is unsurprising given that both types of interactions have the same momentum and velocity dependence and the differences arising from nuclear form factors can easily be compensated by adjusting the velocity integral accordingly [27,29,69].
More surprisingly, we find that both SI and SD are also completely indistinguishable from anapole interactions, even though these interactions depend on both the DM velocity and the momentum transfer in a more complicated way. In other words, if DM scattering is due to an anapole moment, astrophysical uncertainties may lead us to misidentify such an interaction as SI or SD scattering (and vice versa).
This result can be understood as follows. Anapole interactions can be split into two separate parts (see section 3.1): one part proportional to the second velocity integral h(v) and one part proportional to the first velocity integral g(v) and the momentum transfer q 2 . The former contribution predicts a monotonically decreasing spectrum, which differs from SI interactions only because it depends on h(v) rather than g(v), while the latter contribution predicts a spectrum that peaks at finite recoil energy and therefore should lead to strikingly different predictions. It turns out, however, that for a xenon target the relevant nuclear form factor for the second contribution (which is approximately proportional to the magnetic   dipole moment of the nucleus) is much smaller than the SI form factor relevant for the first contribution, so that experimental predictions are completely dominated by the first contribution [53]. For a xenon-based experiment anapole and SI interactions therefore only differ in their dependence on the DM velocity and cannot be distinguished without making an assumption on the velocity integrals g(v) and h(v). As we will see below, this conclusion can change for different target materials.
For magnetic dipole interactions a more non-trivial result is obtained. We find that, if DM scattering is described by SI, SD or anapole interactions, it will be possible to exclude scattering via a magnetic dipole moment at high confidence level. The distinction does however not work in the opposite direction. If DM scattering is due to a magnetic dipole moment, standard SI or SD interactions would still give a good fit to the data (when allowing the DM velocity distribution to vary).
The reason for this observation is that magnetic dipole interactions are enhanced for small recoil energies proportional to inverse powers of the momentum transfer (because they result from the exchange of massless photons). As a result, magnetic dipole interactions fall off more steeply towards larger momentum transfer than SI or SD interactions (for similar velocity integrals). It is therefore possible to interpret scattering from MDM interactions in terms of SI or SD interactions by considering a more steeply falling velocity integral (possibly combined with a smaller value of the DM mass).
To fit a recoil spectrum from SI or SD interactions with MDM interactions, on the other hand, would require a flatter velocity integral and a larger value of the DM mass. However, the monotonicity requirement of the velocity integral implies that the recoil spectrum for MDM interactions cannot be arbitrarily flat. In fact, even if the velocity integral is as flat as possible (i.e. effectively constant), the predicted recoil spectrum is still too steep to give a good fit to the assumed signal. This expectation is confirmed explicitly in the left panel of figure 6, where we show the binned event rates predicted for SI interactions compared to the typical best fit obtained from MDM interactions (as defined above).
For scattering via a dark magnetic dipole moment, we find the converse situation: If DMDM corresponds to the true model of DM, all the other kinds of interaction considered in our study can be ruled out at high confidence level. The reason for this observation is obvious: As shown in the right panel of figure 6, DMDM interactions predict a non-monotonic recoil spectrum, i.e. a peak at finite recoil energy (see figure 6). Given the monotonicity requirement of the velocity integral, such a spectrum can never be fitted with the other interactions that we consider here (see also [33]). Our study shows that the non-monotonic nature of the spectrum can be established with a relatively small number of observed events, i.e. within the level of sensitivity achievable by upcoming direct detection experiments. In this particular scenario, it would therefore be possible to rule out standard SI or SD interactions at high statistical significance without the need to make any assumptions on the astrophysical distributions of DM.
We emphasise that taking DMDM interactions as the fitted model can nevertheless give a good fit to experimental data resulting from other kinds of interactions. The reason is that it is still possible to obtain a monotonically falling recoil spectrum from DMDM interactions, provided the velocity integral falls very steeply. This explains why no distinguishability is found in the final column of figure 5.
To conclude this discussion we note that, although not shown in figure 5, we have also considered DM interactions due to an electric dipole moment or a dark electric dipole moment (i.e. a new heavy mediator coupling to a DM electric dipole moment). We have found these cases to be indistinguishable from MDM and DMDM interactions, respectively, and to yield the same conclusions concerning their distinguishability from other kinds of interactions. As these interactions offer no additional insights, we omit them from figure 5 and similar figures below.

Different couplings to protons and neutrons
So far we have only considered very specific coupling structures. As described in section 3.1, we have always assumed isoscalar couplings (f p = f n ) for the SI case, isovector couplings (g p = −g d ) for the SD case and photon-like couplings to the electromagnetic charge for the other cases. Nevertheless, there are many models (such as DM interacting via a Z with kinetic mixing and mass mixing [70]), where the interactions of DM can have a more complicated structure. To conclude this section, we therefore relax our assumptions on the coupling structures.
Ref. [30] has studied the question to what degree the ratio f n /f p can be determined from future experiments in a halo-independent way. The conclusion was that, by combining information from several experiments, it may be possible to distinguish between SI interactions mediated by the Higgs (f p = f n ) from SI interactions mediated by the Z-boson (f p ∼ 0). Here we want to focus on a different question: Could DM interactions via a non-standard DM operator be misidentified as SI interactions with non-isoscalar couplings? In other words, are there situations where standard SI interactions give a bad fit to the data, but an acceptable fit could be obtained for non-isoscalar couplings?
A particularly interesting case to consider in this context is the one where DM has DMDM interactions. As discussed above, in this case it should be possible to rule out standard SI interactions at high significance. Let us now repeat the analysis from above by treating the ratio f n /f p as an additional free parameter.
At first sight, changing the ratio of the couplings to protons and neutrons should not significantly alter the shape of the recoil spectrum and should therefore not improve the fit to a non-monotonic recoil spectrum. The crucial observation, however, is that in heavy  nuclei like xenon the distribution of protons and neutrons in the nucleus differ slightly, with neutrons having a higher density close to the surface (the so-called neutron skin) [71]. As a result, the form factor for protons falls off more slowly towards larger momentum transfer than the one for neutrons. The difference between the form factors is completely negligible for f p = f n , but it plays an important role in the case of destructive interference between the two contributions (sometimes called isospin-violating DM [72]). Due to the different form factors, the amount of interference then depends sensitively on the nuclear recoil energy. For example, for f n /f p ≈ −0.7, the cancellation between the two contributions is maximal for zero recoil energy and becomes less important for larger recoil energies, so that the predicted recoil spectrum obtains a non-standard shape [71]. In other words, the form factor differences can potentially lead to a non-monotonic recoil spectrum from SI interactions with f n /f p < 0.
Our results are shown in figure 7. We find that, allowing variable f n /f p , it is indeed possible to fit a non-monotonic recoil spectrum like the one obtained from a dark magnetic dipole moment with SI interactions. As expected, the best-fit point is then found to be close to f n /f p ≈ −0.7. We conclude that these two models are essentially indistinguishable and that DMDM interactions may therefore incorrectly be identified as isospin-violating DM (or vice versa).
It should be clear however that such a confusion can only occur if DM is only observed in a single experiment. If the non-monotonic spectrum does indeed result from an almost maximal destructive interference, much larger event rates would be expected in other target materials with different ratios of protons to neutrons, enabling us to easily test such a scenario. We will quantitatively confirm this statement in section 5.2.

Constraining DM coupling structures with several experiments
We have seen in the previous section that in many cases a single liquid-xenon based experiment is insufficient to distinguish different coupling scenarios. In this section, we study whether better discrimination can be achieved if a DM signal is seen in more than one kind of experiment. In principle, the presence of several experiments is expected to lead to significant improvements, as in general different target materials have different form factors and therefore both the shape and the normalisation of the spectrum will depend on the coupling structure in different ways for different experiments. 10 A trivial example would be an argonbased target, which has no sensitivity to SD interactions and could therefore allow perfect discrimination between SI and SD scattering.
A more interesting situation can occur for more complicated models like DM with a magnetic dipole moment, where more than one operator contributes in the non-relativistic limit. In this case, the relative contribution of the different operators may vary from target to target, potentially making it possible to distinguish such a scenario from standard interactions. An example for such a case was already discussed in section 3.2 (see figure 2). We now use the quantitative methods introduced in section 3.3 to study this and other interesting cases more closely.
In addition to the liquid-xenon-based experiment discussed above, we will consider two further target materials representative of next-generation experiments: a germanium (Ge) semiconductor experiment such as EDELWEISS-III [9] or SuperCDMS [8] and a lowbackground iodine (I) crystal scintillator such as SABRE [7]. 11 Our implementation of both types of experiments is summarised in table 1. As before, we make the simplifying assumption that background contributions can be neglected.

Combination of xenon and germanium
Let us first discuss to what extent different coupling scenarios can be distinguished in a halo-independent way by combining a xenon-and a germanium-based experiment. To this end we show in figure 8 the results of our analysis for the five benchmark models introduced in section 3.1, using the experimental set-ups given in table 1. Analogously to figure 5, the blue lines in each panel show the distribution of x 0 , defined as the minimal χ 2 obtained in the fit of a given fitted model (corresponding to the different columns) to the mock data generated by a particular true model (corresponding to the different rows). Furthermore, the distribution of x 1 (shown in red) represents the quality of the fit to the typical prediction of the fitted model (see section 3.3).
We can make two immediate observations from figure 8: Firstly, the fourth column shows that the MDM scenario does not provide a good fit to the data generated for any other model. Secondly, the last row implies that if the true model of DM is given by DMDM interactions, one can distinguish it with high confidence level from any of the other scenarios, in particular from the standard SI or SD interactions. Comparing with figure 5, we observe that these pairs of models could already be distinguished with a xenon-based experiment alone (albeit in some cases with a smaller value of the distinguishability parameter D and thus with larger similarity S), and we refer to section 4.1 for a physical interpretation of these cases. Most importantly, we find that adding the germanium target on top of the xenonbased experiment does not allow for the discrimination of any additional pair of models. At first sight, this is unexpected: in general, the individual target nuclei have different relative sensitivities to different models, and hence combining the information from two or more experiments should break any degeneracy between them. For the combination of xenon and germanium targets, however, we find that even increasing the germanium exposure by a factor of 10 (such that the number of observed events in xenon and germanium are comparable) the distinguishability does not improve significantly. Let us therefore pause our discussion of figure 8 to discuss a simple argument to make sense of this observation, which is based on the method for halo-independent comparison introduced in [34].

Combination of xenon and germanium
Assuming that for a given DM mass the two experiments based on xenon and germanium are sensitive to a similar range of v min values, the rates expected for SI couplings can schematically be written as Similarly, the rates induced by SD couplings are given by In a region of v min -space accessible to both experiments, we can define the ratio of rates for SD interactions as and analogously RR SI (Xe/Ge) for SI interactions.
In order for SI and SD interactions to be distinguishable in a halo-independent way, it is therefore necessary that RR SI (Xe/Ge) is sufficiently different from RR SD (Xe/Ge), so that only one of these scenarios can give a good fit to the data for the assumed DM mass. 12 It turns out, however, that numerically the ratios of the sensitivities are, by coincidence, very similar: where η T is the natural abundance of the target isotope T , S n are the expectation values of the spin-content of the proton and neutron group in the nucleus T , respectively, and J T is its total spin. We thus find While this argument is of course simplified, as it does not take into account effects arising e.g. from Poisson fluctuations in the data, from the form factor suppression at finite q 2 , or from the fact that usually one has several bins with different degree of overlap in the v min space, equation (5.6) still gives a qualitative explanation for the similarity of the SI and SD interactions observed in figure 8.
Returning to the discussion of figure 8, we observe that also in the cases where the true model is given by AM or MDM interactions, a good (halo-independent) fit to the data can be obtained for standard SI or SD interactions. Consequently, one cannot discriminate between these scenarios. The reason is similar to the one discussed in section 4.1 for the xenon-based experiment alone. For scattering off germanium the cross section induced by the anapole moment is dominated by the v 2 ⊥ term, with a dependence on the target nucleus similar to the SI case. Hence, RR AM (Xe/Ge) is very similar to the corresponding ratios for SI or SD interactions. The different dependence on the DM velocity does not increase the distinguishability, as they can be fully compensated by choosing an appropriate velocity distribution in the fit. A similar argument holds when replacing the AM interaction by the MDM interaction as the true model. Again, we have confirmed that our conclusions do not change substantially when increasing the exposure of the germanium experiment by a factor of 10, and are hence "intrinsic" to this set of target nuclei.
Lastly, we observe that all models under discussion can be well fitted to by the DMDM model, as shown in the final column of figure 8. This observation seems surprising, because for DMDM interactions event rates are proportional to q 4 = 16 v 4 min µ 4 T χ , so that scattering is enhances for heavy targets: RR DMDM (Xe/Ge) RR SI (Xe/Ge). One would thus expect the information from the two different targets to be highly complementary. It turns out, however, that in this particular case the difference in the ratio of rates can be fully compensated by changing the assumed DM mass. For example, for true model SI and fitted model DMDM, a good fit to the data from both xenon and germanium can typically be obtained for m χ ≈ 10 GeV (compared to an assumed true DM mass of 50 GeV). As for the case of a xenon target alone, the additional momentum dependence of the DMDM interaction is then compensated by a steeply falling velocity integral.

Different couplings to protons and neutrons
Let us now turn to an example where adding the information from a germanium target makes an important difference for the discrimination of different models. For this purpose we return to the question of whether a non-standard DM-nucleon interaction such as DMDM could be confused with an SI scenario with variable neutron-to-proton coupling ratio f n /f p . In section 4.2 we found that due to the slightly different form factors for protons and neutrons, the recoil spectrum in a single experiment can be non-monotonic even for SI interactions, provided that f n /f p is close to the value leading to maximal destructive interference. In particular we concluded from figure 7 that an SI interaction with variable neutron-to-proton coupling ratio can provide a good (halo-independent) fit to mock data generated for DMDM interactions (considering only a xenon experiment).
The value of f n /f p leading to maximally destructive interference does however differ from one target nucleus to the other: for example, it is approximately −0.7 for xenon, but −0.79 for germanium. In particular, we find that there is no value of the neutron-toproton coupling ratio for which the recoil spectrum expected from SI interactions is nonmonotonic both at a xenon-and a germanium-based experiment. Consequently, we expect the combination of these two experiments to be able to rule out an SI interaction scenario with variable f n /f p , when the true model is DMDM. This expectation is confirmed in figure 9, where we show the corresponding distributions of x 0 and x 1 . Clearly there is a large haloindependent distinguishability of the two interaction scenarios.
An intuitive way to understand this result is presented in the right panel of figure 9 using the mapping to v min -space introduced in section 3.2. For this figure we have generated mock data (without Poisson fluctuations) based on the DMDM model and then performed the mapping to v min -space for SI interactions with f n /f p = −0.7. For the xenon-based experiment, this reconstruction yields monotonically decreasing measurements of the velocity integral, allowing for a consistent fit of all data points. The values of the velocity integral inferred from germanium, however, can be seen to increase towards larger values of v min (although for the assumed exposures this increase is not statistically significant). More importantly, since the destructive interference is maximal in xenon, but less pronounced in germanium, the velocity integrals inferred from the two experiments are in strong tension with each other. This tension ultimately allows to exclude the assumed value of f n /f p used for the reconstruction in a halo-independent way.

Combination of xenon and iodine
Let us now turn to the discussion of the results obtained from combining the information of a xenon-and an iodine-based experiment (figure 10). The upper left 2 × 2 panels indicate that this set of experiments also cannot distinguish between standard SI and SD interactions with large significance in a halo-independent way. As above, this can be qualitatively explained by estimating the ratios of rates Interestingly, we find that if the true model is given by AM interactions, adding the information from the iodine experiment helps substantially in discriminating this scenario from the standard SI or SD interaction. This observation can be explained by the fact that for AM scattering off iodine the term ∝ q 2 in the scattering cross section dominates due to the large magnetic dipole moment of the iodine nucleus [73]. This term gives rise to a peak in the recoil spectrum (see also figure 9 in [53]), which cannot be compensated by adjusting the velocity distribution, leading to a worse fit for models without such a peak. Furthermore, as the q 2 contribution is only important for iodine, the ratio of the rates for scattering off xenon and iodine is significantly different for AM and for SI (or SD) interactions, leading to a (halo-independent) tension between the experiments. Together these two observations explain the large distinguishability for the case that the true model has AM interactions and the fitted model is SI or SD.
If, on the other hand, the true model is given by MDM interactions the additional 1/q 2 term in the cross section induced by the photon propagator always leads to a monotonically decreasing spectrum. Nevertheless, the second argument from above remains valid. The different magnitude of the nuclear magnetic dipole moment for iodine and xenon still leads to a halo-independent tension when incorrectly trying to fit MDM interactions with SI or SD interactions (see figure 2). Hence, there is still some discrimination power (D = 0.79 and 0.64, respectively) for this case.

Conclusions
Near-future experiments for the direct detection of DM promise significant improvements in sensitivity over existing searches. In anticipation of these improvements it is timely to develop and refine strategies for extracting the particle physics properties of DM from data. Because of the strong dependence of experimental event rates on the essentially unknown Galactic velocity distribution of DM it is of particular importance to devise halo-independent methods, i.e. methods that enable us to deduce properties of the DM particle from a positive detection without the need to make any assumptions on the astrophysical distributions.
In this work we have presented a very general framework for analysing the data from one or several direct detection experiments independent of astrophysical assumptions and for quantifying the goodness-of-fit of a given particle physics scenario. Following [30,31], we perform a halo-independent fit to a given set of data by parametrizing the velocity integral as a piecewise constant function with a sufficiently large number of steps. As we determine the best-fit velocity distribution directly from data instead of fixing it to e.g. a Maxwell-Boltzmann distribution, our approach can be employed in order to test whether there exists any velocity distribution for which a given particle physics scenario of DM is compatible with the data. Importantly, we have shown that this formalism can be applied universally to any combination of non-relativistic operators describing the interactions of DM with nuclei, in particular to models with non-standard dependence of the scattering cross section on the momentum exchange and the DM velocity. This enables us to consider a much broader set of particle physics models compared to many existing halo-independent analyses in the literature.
To demonstrate how this method can be applied to the case of a positive detection of DM in one or several experiments, we have studied mock data generated for realistic near-future direct detection experiments. To this end we have considered a range of different possibilities for the true particle physics nature of DM, including the standard spin-independent and spindependent interaction as well as scenarios involving a non-trivial velocity and momentum dependence of the DM-nucleon scattering cross section. For each mock data set we then perform halo-independent fits to the data for various different assumptions on the properties of the DM particle (see figures 1 and 2).
By repeating this fitting procedure for a large number of Poisson realisations of mock data sets, we consistently take into account the statistical noise expected in the observed data. We can then determine which of the interaction scenarios can be distinguished without specifying the velocity distribution. To quantify the similarity and distinguishability of different DM interaction scenarios we have developed a simplified procedure based on the typical predictions of the fitted model ( figure 3). We have confirmed the validity of this approximation using explicit Monte Carlo simulations ( figure 4).
Interestingly we find that already a single experiment like XENON1T with a moderate exposure of two ton-years can be employed for inferring non-trivial information about the particle physics properties of DM in a fully halo-independent way (figures 5 and 6). For example, our results show that a model of DM predicting a non-monotonic recoil spectrum (realized e.g. for a dark magnetic dipole moment) can be clearly distinguished from the standard spin-independent or spin-dependent interaction scenario with a XENON1T-like experiment. Also, if DM interacts via the standard spin-independent or spin-dependent coupling, it is possible to rule out in a halo-independent way long-range interactions of DM with nucleons, which could be induced e.g. by a magnetic dipole moment of DM.
Lastly we have studied to what extent a detection of DM in more than one experiment helps in discriminating different particle physics scenarios. We have found that adding a germanium target to the xenon-based experiment does not increase significantly the haloindependent distinguishability of the scenarios discussed in this work (figure 8). We have illuminated the numerical results by a simplified analytical argument explaining why the complementarity of a xenon and germanium target is strongly limited. An exception to this statement is the case in which the data of both experiments are fitted assuming a spinindependent interaction with variable neutron-to-proton coupling ratio f n /f p , which can be constrained much more tightly by two experiments with different target nuclei (figures 7 and 9). On the other hand, we have shown that the presence of a DM signal in an iodinebased experiment in addition to the xenon-based detector yields additional discrimination power ( figure 10). In particular, if scattering is induced by an anapole or magnetic dipole moment of DM, it should be possible to rule out the standard assumption of spin-independent or spin-dependent interactions by combining the information from both experiments, again without referring to a particular velocity distribution of DM in the galaxy.
Our analysis is based on rather modest assumptions on the exposures of upcoming experiments. Clearly, stronger discrimination power can be achieved by more ambitious experiments and by combining the information from more than two experiments. Nevertheless, hundred DM events in a liquid xenon experiment and ten events in an additional experiment may already be sufficient to extract highly non-trivial information on the particle physics nature of DM in a completely halo-independent way. While the conclusion may simply be that some more exotic models of DM can be ruled out, future experiments may also point in the opposite direction and tell us that the interactions of DM are much more complicated than what is usually assumed.
Strings and the Early Universe, as well as by the Studienstiftung des Deutschen Volkes and the TUM Graduate School.

A Non-relativistic effective interactions
Here and below, we follow largely the conventions and notation from [57]. It has been shown in [49] that in the non-relativistic limit the interactions between the DM particle and quarks can be written as a linear combination of 15 effective operators where τ = 0 for isoscalar interactions (i.e. equal couplings to protons and neutrons) and τ = 1 for isovector interactions (opposite couplings to protons and neutrons). 13 Additional operators can be constructed by multiplying this basic set with (inverse) powers of q 2 . Such composite operators can arise for example from long-range interactions or from derivative interactions and are typically present if DM interacts via an electric or magnetic dipole moment. We include the three simplest cases in our discussion by expanding the coefficients c where m N ≡ m p ≈ m n is the nucleon mass. The interactions of DM with nucleons can therefore be fully described by the 14 × 2 × 3 coefficients z (τ ) κ,λ , which we combine into a coefficient vector z for simplicity. We note that our formalism assumes that the mediator(s) of DM scattering are either very light or very heavy compared to the typical momentum transfer. We do not discuss the case of a mediator with m 2 med ∼ q 2 . In order to calculate the differential scattering cross section dσ/dE R for a given target nucleus, the next step is to calculate the so-called DM response functions S (τ,τ ) α for α = 1, . . . , 8, which depend on the coefficient vector z in a known way. The DM response functions then need to be multiplied with the corresponding nuclear form factors W where the sum over S is also referred to as the transition probability. We will now discuss the two ingredients for calculating this quantity in detail.
A.1 Dark matter response functions S (τ,τ ) α The DM response functions S (τ,τ ) α are defined in appendix A of [57]. The index α = 1, . . . , 8 runs over the eight different response functions in the order as they are defined in equation (A.1) of [57]. The response functions depend on the coefficient vector z, the relative velocity v and momentum transfer q. The crucial observation for our purposes is that it is possible to decompose the response functions as a sum of a velocity-independent term and a term proportional to v 2 : The functions S (τ,τ ) β,α q 2 (with β ∈ {1, 2}) can furthermore be written as a linear combination of integer powers of q 2 (see [57]): where we have introduced the dimensionless quantity and implicitly defined the coefficients σ (τ,τ ) β,α,l for l = −2, . . . , 5. These 2×8×8×2×2 coefficients fully characterise the DM response and can be calculated in terms of the coefficient vector z. This calculation, although cumbersome, does not pose any fundamental difficulties and will therefore not be discussed further. For the remainder of the discussion it will be sufficient to simply think of σ (τ,τ ) β,α,l (z) as the quantity describing the particle physics properties of DM.
A.2 Nuclear form factors W (τ,τ ) α Apart from the DM response functions, the differential cross section defined in equation (A.3) also depends on the nuclear form factors W (τ,τ ) α , which are provided in [50]. 14 These form factors depend on q 2 via the combination y = 1 The nuclear form factors are given as fit functions of the form Defining y(a) = 1 4 b 2 A m 2 N · a, this can be rewritten in the form 14 For convenience we define W The decomposition of the differential cross section into different powers of v 2 is essential for the halo-independent formalism discussed in this paper. The additional decomposition into different powers of a will be extremely convenient for reducing the number of numerical integrations necessary for a concrete implementation of this formalism, as discussed in below.

A.4 Efficient numerical implementation
In order to determine the matrix D ij for a given experiment, we need to calculate the matrices In principle, the necessary numerical integrations are straight-forward once the details of the DM interactions have been specified. Nevertheless, these calculation can be computationally rather expensive, in particular if N s is very large. Fortunately it is possible to significantly reduce the number of necessary integrations by splitting the calculation into several separate steps. Comparing these expressions for G ij and H ij to equation (A.11) and noting that the coefficients σ (τ,τ ) β,α,l (z) and w (τ,τ ) α,k do not depend on the recoil energy, it becomes clear that we first need to calculate the functions where a = q 2 /m 2 N = 2 m T E R /m 2 N . These functions can be thought of as the predicted number of events in the ith bin for a hypothetical scattering cross section with simplified energy dependence of the form dσ/dE R ∼ exp(−E R )E p R and a velocity integral of the form g(v min ) = Θ(a − a(v min )), where a(v min ) = 4µ 2 T χ v 2 min /m 2 N . We therefore refer to these functions as reduced event rates. The reduced event rates contain all relevant information on the experiments under consideration (such as the exposure, the energy-dependent acceptance and the energy resolution). For a given target, the functions d i (a) need to be calculated for p ∈ {−2, −1, . . . , 15}. Crucially, they are independent of the DM mass, the coefficient vector z and of the binning of v min space. Hence they only have to be calculated once for every experiment (for a given binning in recoil energy) and can then be used for arbitrary DM mass, combination of effective operators and binning of v min -space.
Once the particle physics properties of DM have been specified in terms of the coefficients σ (τ,τ ) β,α,l introduced above, it is possible to combine the reduced event rates d where the coefficients w (τ,τ ) α,k contain the details of the nuclear response (see above). Finally, once the binning v 1 , . . . , v Ns of v min -space has been specified, one simply needs to calculate the corresponding values a j ≡ a(v j ). The matrix G ij is then given by 15 16) and the matrix H ij in analogy. The key point is that the second and third step require no further numerical integration and can therefore be performed very fast even for complicated non-standard interactions. Splitting the calculation of G ij and H ij into the three steps discussed above therefore allows to scan over the parameters characterising the particle physics properties of DM in a very efficient way.
B Calculating event rates from the velocity integral

B.1 The second velocity integral
We consider the case that the velocity integral g(v min ) is piecewise constant, so that it can be written as where Θ(x) is the Heaviside step function and l j ≡ g j − g j+1 (with g Ns+1 = 0). In order to calculate h(v min ) for given v min , we choose k in such a way that v k−1 < v min < v k . One then obtains: We observe that h(v min ) is also a piecewise-constant function, i.e. we can write h( The relationship between h j and g j can be written as h j = F jj g j with F jj = v 2 j+1 , F jj = 0 for j > j and F jj = v 2 j +1 − v 2 j for j > j. Using the matrix F jj and the matrices G ij and H ij defined in eqs. (2.15) and (2.16), we can then construct the matrix D ij that relates the velocity integral to the predicted number of events in a given bin:

B.2 The impact of high-velocity bins
Let us now discuss how to determine the range of v min -space that should be used to calculate the matrices G ij and H ij . We note that both G ij and H ij are non-zero only if the energy E i of the bin under consideration is comparable to the physical recoil energy E j corresponding to the jth bin in v min -space. How close E i and E j have to be in order for G ij and H ij to be non-zero depends on the energy resolution ∆E R , which determines the probability that a physical recoil energy E j leads to an observed energy E i . The important point is that there is always a v min and a v max such that G ij ≈ 0 and H ij ≈ 0 for v j < v min or v j > v max . This observation suggests that one can simply take v 1 = v min and v Ns+1 = v max . However, things become more complicated when the matrices G ij and H ij are combined into the matrix D ij = G ij + j H ij F j j . To illustrate the problem, let us consider a range of v min -space [v 1 , v Ns+1 ] such that E R (v Ns+1 ) is significantly larger than v max . In this case there is a j max such that for J > j max one obtains G iJ = 0 and H iJ = 0 for all bins i. For such a J it then follows that because F j j = v 2 j+1 − v 2 j independent of j for j > j . Crucially we find that D iJ can be non-zero even for J > j max .
The calculation above shows that arbitrarily high bins in v min -space can influence the predictions even for the lowest bins in reconstructed energy. This observation seems paradoxical at first, but it is a direct consequence of the relation between the velocity integrals g(v min ) and h(v min ), which imply that a local variation in g(v min ) leads to a global variation in h(v min ), see also [39]. To be more specific, if we vary g(v min ) only in the interval [v j , v j+1 ] and keep it fixed everywhere else, h(v min ) will vary for all velocities v min < v j+1 and the matrix F ij encodes this behaviour.
At first sight, the observation that D ij remains non-zero for arbitrarily large v j seems to imply that an arbitrarily large number of steps is needed for a fully general treatment, which would pose a significant problem for our approach. The simplest solution to this problem would be to impose some upper limit on the allowed velocity range, for example motivated by the observed galactic escape velocity or simply by the requirement that DM ought to be non-relativistic. As we will now show, however, it is possible to solve this problem more elegantly.
The important observation is that D iJ takes a particularly simple form. Defining H i as the sum of all entries in the ith row of H ij , i.e. H i ≡ j H ij , we obtain In other words, for fixed J > j max the column vector D iJ is simply proportional to the vector H i . As a result, for any vector (g j ) = (g 1 , . . . , g Ns ) in v min -space we can define a new vector (g j ) = (g 1 , . . . , g jmax ,g, 0, . . . , 0) with such that j D ij g j = j D ijgj . While (g j ) is assumed to be monotonically decreasing (i.e. g j+1 ≤ g j ), this no longer needs to be true for (g j ), because it is possible thatg > g jmax . In other words, adding a large number of monotonically decreasing bins beyond j max is equivalent to adding a single bin which does not need to satisfy the monotonicity requirement. This observation implies that it is not necessary to consider arbitrarily large values of v min in order to find the best-fit velocity integral.