# Studying generalised dark matter interactions with extended halo-independent
methods

Felix Kahlhoefer  and Sebastian Wild

###### Abstract

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.

##  1 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.

Answering these questions is complicated by the fact that event rates in
direct detection experiments depend in complicated ways on the velocity
distribution of DM particles in the Galactic halo, which is subject to large
uncertainties [15, 16, 17, 18]. Many studies have investigated the impact of
these uncertainties on our ability to employ the results of future direct
detection experiments in order to infer the mass of the DM particle or to
discriminate between different DM models (e.g. by distinguishing between spin-
independent and spin-dependent interactions or by determining separately the
DM-proton and the DM-neutron coupling) [19, 20, 21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31]. In particular, a range of so-called halo-independent methods
have been developed to reduce or even eliminate completely the impact of
astrophysical uncertainties on the DM properties that can be inferred from
direct detection experiments [32, 33, 34, 35, 36, 37, 38, 39, 40, 30, 31, 41,
42, 43, 44, 45, 46, 47, 48].

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 [51, 52, 53, 54, 55, 56,
57, 58, 46, 59, 60]. The possible importance of the additional DM-nucleon
scattering operators makes the issue of astrophysical uncertainties even more
pressing, because a non-standard 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
non-relativistic 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 $i$th 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.

##  2 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
$\mathrm{d}\sigma/\mathrm{d}E_{\mathrm{R}}\propto 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.

###  2.1 Calculating event rates from velocity integrals

The differential event rate with respect to recoil energy in a direct
detection experiment is given by

$\frac{\text{d}R}{\text{d}E_{\text{R}}}=\frac{\rho}{m_{T}\,m_{\chi}}\int_{v_{\text{min}}}^{\infty}vf(\boldsymbol{v}+\boldsymbol{v}_{\text{E}}(t))\frac{\text{d}\sigma}{\mathrm{d}E_{\text{R}}}\text{d}^{3}v\;,$
(2.1)

where $\rho$ is the local DM density, $m_{\chi}$ and $m_{T}$ are the DM and
target nucleus mass, $f(\mathbf{v})$ is the local DM velocity distribution
evaluated in the Galactic rest frame, $\boldsymbol{v}_{\text{E}}(t)$ is the
velocity of the Earth relative to the Galactic rest frame and
$v=|\boldsymbol{v}|$. For elastic scattering the _minimum velocity_ required
for a DM particle to transfer the energy $E_{\text{R}}$ to a given nucleus $T$
is

$v_{\text{min}}(E_{\text{R}})=\sqrt{\frac{m_{T}E_{\text{R}}}{2\,\mu_{T\chi}^{2}}}\;,$
(2.2)

where $\mu_{T\chi}=m_{T}\,m_{\chi}/(m_{T}+m_{\chi})$ is the reduced mass of
the DM-nucleus system.

The differential cross section $\text{d}\sigma/\mathrm{d}E_{\text{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:

$\frac{\text{d}\sigma}{\mathrm{d}E_{\text{R}}}=\frac{\text{d}\sigma_{1}}{\mathrm{d}E_{\text{R}}}\frac{1}{v^{2}}+\frac{\text{d}\sigma_{2}}{\mathrm{d}E_{\text{R}}}\;.$
(2.3)

Defining the velocity integrals111Here we make two approximations. First, we
neglect the small time dependence of $\mathbf{v}_{\text{E}}$ and second we
assume that $\text{d}\sigma/\mathrm{d}E_{\text{R}}$ only depends on $v^{2}$
and not on the direction of $\mathbf{v}$ (e.g. because the target is
unpolarised and the direction of recoil tracks cannot be distinguished).

$g(v_{\text{min}})=\int_{v_{\text{min}}}^{\infty}\frac{1}{v}\,f(\boldsymbol{v}+\boldsymbol{v}_{\text{E}}(t))\,\text{d}^{3}v\;,\quad
h(v_{\text{min}})=\int_{v_{\text{min}}}^{\infty}v\,f(\boldsymbol{v}+\boldsymbol{v}_{\text{E}}(t))\,\text{d}^{3}v$
(2.4)

we can then simply write the differential event rate as

$\frac{\text{d}R}{\text{d}E_{\text{R}}}=\frac{\rho}{m_{T}\,m_{\chi}}\left(\frac{\text{d}\sigma_{1}}{\mathrm{d}E_{\text{R}}}g(v_{\text{min}})+\frac{\text{d}\sigma_{2}}{\mathrm{d}E_{\text{R}}}h(v_{\text{min}})\right)\;.$
(2.5)

If the second term of this equation is absent, all information on the DM
velocity distribution is encoded in the velocity integral $g(v_{\text{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
$\left[E_{i},E_{i+1}\right]$. In order to study how $R_{i}$ depends on the
velocity integral, it was proposed in [30] to parametrize $g(v_{\text{min}})$
as a piecewise constant function with $N_{\mathrm{s}}$ steps.222Note that this
is equivalent to writing the speed distribution $f(v)$ as a sum of
$\delta$-functions. In other words, we decompose the velocity distribution
into a sum of streams with fixed velocities and negligible velocity
dispersion. By choosing $N_{\mathrm{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_{\text{min}}$-space probed by
the experiments under consideration into $N_{\mathrm{s}}$ intervals of the
form $\left[v_{j},\,v_{j+1}\right]$, where the number of steps
$N_{\mathrm{s}}$ can be as large as 50. In each interval, the velocity
integral is taken to be constant:

$\displaystyle g(v_{\text{min}})=g_{j}\quad\text{for
}\,v_{\text{min}}\in\left[v_{j},v_{j+1}\right]\,.$ (2.6)

The fact that the velocity integral is both positive and monotonically
decreasing leads to the requirement $0\leq g_{j}\leq 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}=\sum_{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.

###  2.2 Halo-independent methods for generalised dark matter interactions

At first sight, the situation becomes more complicated if both
$g(v_{\text{min}})$ and $h(v_{\text{min}})$ contribute to the differential
event rate. The crucial observation is however that $g(v_{\text{min}})$ and
$h(v_{\text{min}})$ are not independent. In fact, we can make use of the fact
that $g^{\prime}(v)=-f(v)/v$, where $f(v)=\int
v^{2}\,f(\mathbf{v}+\mathbf{v}_{\text{E}})\,\mathrm{d}\Omega_{v}$, to write
[39, 47]

$\displaystyle h(v)$
$\displaystyle=-\int_{v_{\text{min}}}^{\infty}v^{2}\,g^{\prime}(v)\mathrm{d}v$
$\displaystyle=\left[-g(v)\,v^{2}\right]_{v_{\text{min}}}^{\infty}+\int_{v_{\text{min}}}^{\infty}2\,v\,g(v)\mathrm{d}v\;.$
(2.7)

In other words, it is possible to calculate $h(v_{\text{min}})$ directly from
the usual velocity integral $g(v_{\text{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

$g(v_{\text{min}})=\sum_{j=1}^{N_{\mathrm{s}}}l_{j}\Theta(v_{j+1}-v_{\text{min}})\;,$
(2.8)

where $\Theta(x)$ is the Heaviside step function and $l_{j}\equiv
g_{j}-g_{j+1}$, with the $g_{j}$ defined via equation (2.6), setting in
addition $g_{N_{\mathrm{s}}+1}\equiv 0$. We can then use equation (2.7) to
calculate $h(v_{\text{min}})$:

$h(v_{\text{min}})=\sum_{j=1}^{N_{\mathrm{s}}}l_{j}\,v_{j+1}^{2}\,\Theta(v_{j+1}-v_{\text{min}})\;,$
(2.9)

where details of the calculation are provided in appendix B.1. Most
importantly, we find that $h(v_{\text{min}})$ is also piecewise constant and
monotonically decreasing, i.e. we can write $h(v_{\text{min}})=h_{j}$ for
$v_{\text{min}}\in\left[v_{j},\,v_{j+1}\right]$ with $0\leq h_{j}\leq
h_{j-1}$. The parameters $g_{j}$ and $h_{j}$ are then related by a simple
linear transformation:

$h_{j}=\sum_{j^{\prime}}F_{jj^{\prime}}g_{j^{\prime}}\;,$ (2.10)

with the matrix $F$ given by $F_{jj}=v_{j+1}^{2}$, $F_{jj^{\prime}}=0$ for
$j>j^{\prime}$ and $F_{jj^{\prime}}=v_{j^{\prime}+1}^{2}-v_{j^{\prime}}^{2}$
for $j^{\prime}>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_{\mathrm{D}};E_{\mathrm{R}})$
that a DM interaction with true (i.e. physical) recoil energy $E_{\mathrm{R}}$
leads to a signal with detected energy $E_{\mathrm{D}}$. In terms of this
function, and using our ansatz for $g(v_{\text{min}})$, we find

$\displaystyle\frac{\mathrm{d}R}{\mathrm{d}E_{\mathrm{D}}}$
$\displaystyle=\int_{0}^{\infty}\frac{\text{d}R}{\text{d}E_{\text{R}}}\,p(E_{\mathrm{D}};E_{\mathrm{R}})\,\mathrm{d}E_{\mathrm{R}}$
$\displaystyle=\frac{\rho}{m_{T}\,m_{\chi}}\int_{0}^{\infty}\left(\frac{\text{d}\sigma_{1}}{\mathrm{d}E_{\text{R}}}g(v_{\text{min}})+\frac{\text{d}\sigma_{2}}{\mathrm{d}E_{\text{R}}}h(v_{\text{min}})\right)\,p(E_{\mathrm{D}};E_{\mathrm{R}})\,\mathrm{d}E_{\mathrm{R}}$
$\displaystyle=\frac{\rho}{m_{T}\,m_{\chi}}\sum_{j}\left[g_{j}\int_{E_{j}}^{E_{j+1}}\frac{\text{d}\sigma_{1}}{\mathrm{d}E_{\text{R}}}\,p(E_{\mathrm{D}};E_{\mathrm{R}})\,\mathrm{d}E_{\mathrm{R}}+h_{j}\int_{E_{j}}^{E_{j+1}}\frac{\text{d}\sigma_{2}}{\mathrm{d}E_{\text{R}}}\,p(E_{\mathrm{D}};E_{\mathrm{R}})\,\mathrm{d}E_{\mathrm{R}}\right]\;,$
(2.11)

where $E_{j}$ is defined implicitly by $v_{j}=v_{\text{min}}(E_{j})$.

If we are interested in the number of events $R_{i}$ predicted in a given bin,
we simply need to multiply $\mathrm{d}R/\mathrm{d}E_{\mathrm{D}}$ with the
total exposure $\kappa$ and integrate over $E_{\mathrm{D}}$:

$R_{i}=\kappa\int_{E_{i}}^{E_{i+1}}\frac{\mathrm{d}R}{\mathrm{d}E_{\mathrm{D}}}\,\mathrm{d}E_{\mathrm{D}}\;.$
(2.12)

We now make the simplifying assumption that the detection efficiency depends
only on the physical recoil energy $E_{\mathrm{R}}$ (and not on
$E_{\mathrm{D}}$) and that fluctuations can be approximated by a Gaussian
distribution with standard deviation given by $\Delta E(E_{\mathrm{R}})$. In
this case we can write

$p(E_{\mathrm{D}};E_{\mathrm{R}})=\epsilon(E_{\mathrm{R}})\,\frac{1}{\sqrt{2\pi}\Delta
E(E_{\mathrm{R}})}\exp\left[-\frac{(E_{\mathrm{D}}-E_{\mathrm{R}})^{2}}{2\Delta
E(E_{\mathrm{R}})^{2}}\right]\;.$ (2.13)

The integration over $E_{\mathrm{D}}$ can then be performed analytically, and
one obtains

$R_{i}=\sum_{j}G_{ij}g_{j}+\sum_{k}H_{ik}h_{k}\;,$ (2.14)

where

$\displaystyle G_{ij}$
$\displaystyle=\frac{\kappa\,\rho}{2m_{T}\,m_{\chi}}\int_{E_{j}}^{E_{j+1}}\frac{\text{d}\sigma_{1}}{\mathrm{d}E_{\text{R}}}\,\epsilon(E_{\mathrm{R}})\left[\text{erf}\left(\frac{E_{i+1}-E_{\text{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)-\text{erf}\left(\frac{E_{i}-E_{\text{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)\right]\,\mathrm{d}E_{\mathrm{R}}\,,$ (2.15)
$\displaystyle H_{ij}$
$\displaystyle=\frac{\kappa\,\rho}{2m_{T}\,m_{\chi}}\int_{E_{j}}^{E_{j+1}}\frac{\text{d}\sigma_{2}}{\mathrm{d}E_{\text{R}}}\,\epsilon(E_{\mathrm{R}})\left[\text{erf}\left(\frac{E_{i+1}-E_{\text{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)-\text{erf}\left(\frac{E_{i}-E_{\text{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)\right]\,\mathrm{d}E_{\mathrm{R}}\,,$ (2.16)

with $\text{erf}(x)=(2/\sqrt{\pi})\int_{0}^{x}\exp(-t^{2})\,\mathrm{d}t$.

We can further simplify equation (2.14) by making use of the relation
$h_{k}=\sum_{j}F_{kj}g_{j}$ derived above. It is then possible to combine the
two matrices $G_{ij}$ and $H_{ij}$ into one matrix

$D_{ij}=G_{ij}+\sum_{k}F_{kj}H_{ik}\;,$ (2.17)

such that

$R_{i}=\sum_{j}D_{ij}g_{j}\;.$ (2.18)

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

$-2\log\mathcal{L}=2\sum\left[R_{i}+B_{i}-N_{i}+N_{i}\log\frac{N_{i}}{R_{i}+B_{i}}\right]\;.$
(2.19)

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_{\mathrm{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.333We emphasize that there may in principle
be several degenerate minima, i.e. different velocity distributions that yield
the same likelihood. This does not pose a difficulty for our approach, as we
are interested only in the minimum value of the likelihood and not in the
underlying velocity distribution. We refer to [30] for a more detailed
discussion.

The equations above generalise directly to the case of several experiments
(labelled by $\alpha=1,\ldots,N_{\text{exp}}$). In this case, the matrix
$D^{(\alpha)}_{ij}$ must be calculated separately for each experiment.444Note
that the number of bins and therefore the range of the index $i$ may be
different for each experiment. However, the number of steps in
$v_{\text{min}}$-space and hence the range of the index $j$ is the same for
all experiments. One can then calculate the likelihood
$\mathcal{L}^{(\alpha)}$ for each experiment given the expected number of
background events $B_{i}^{(\alpha)}$ and the total number of observed events
$N_{i}^{(\alpha)}$. 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_{\text{min}}$-space for DM interactions with non-standard
velocity dependence.

##  3 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.

###  3.1 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 DM-nucleon Lagrangian, expressed via the non-relativistic operators
$\mathcal{O}_{\lambda}^{\tau}$ [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_{\text{DM}}=50\,$GeV, as well as a standard Maxwell-Boltzmann
velocity distribution, with the most probable speed
$v_{0}=220\,\text{km}/\text{s}$, the Galactic escape velocity
$v_{0}=544\,\text{km}/\text{s}$ and the mean velocity of the Earth
$v_{\text{E}}=230\,\text{km}/\text{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_{N}^{\text{(SI)}}\bar{\chi}\chi\,\bar{N}N$
and
$f_{N}^{\text{(SD)}}\bar{\chi}\gamma^{\mu}\gamma^{5}\chi\,\bar{N}\gamma_{\mu}\gamma^{5}N$,
respectively. In the notation of [57], the resulting non-relativistic
Lagrangians read

$\displaystyle\mathcal{L}_{\text{eff}}=\sum_{N=p,n}f_{N}^{\text{(SI)}}\mathcal{O}_{1}^{(N)}\quad\text{and}\quad\mathcal{L}_{\text{eff}}=-4\sum_{N=p,n}f_{N}^{\text{(SD)}}\mathcal{O}_{4}^{(N)}\,,$
(3.1)

respectively.555We refer to appendix A for the conversion of such a linear
combination of non-relativistic operators into a scattering cross section
$\text{d}\sigma/\text{d}E_{\mathrm{R}}$. 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 $\sigma_{p}^{\text{(SI)}}=5\cdot
10^{-46}\,\text{cm}^{2}$. For spin-dependent scattering we consider a model
with $f_{n}^{\text{(SD)}}/f_{p}^{\text{(SD)}}=-1$ (motivated by a scenario of
Majorana DM interacting with nucleons via the Standard Model $Z$-boson) with a
total cross section $\sigma_{p}^{\text{(SD)}}=1.2\cdot
10^{-40}\,\text{cm}^{2}$. 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
$\mathcal{A}\,\bar{\chi}\gamma^{\mu}\gamma^{5}\chi\,\partial^{\nu}F_{\mu\nu}$,
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]

$\displaystyle\mathcal{L}_{\text{eff}}=2\mathcal{A}e\sum_{N=p,n}\left(Q_{N}\mathcal{O}_{8}^{(N)}+\widetilde{\mu}_{N}\mathcal{O}_{9}^{(N)}\right)\,,$
(3.2)

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 $\widetilde{\mu}_{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_{\perp}^{2}\equiv
v^{2}-q^{2}/(4\mu_{T}^{2})$, 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_{\text{DM}}=50\,$GeV, as
well as an anapole moment $\mathcal{A}=3.5\cdot 10^{-5}\,\text{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
$\mathcal{D}_{\text{magn}}\,\bar{\chi}\sigma^{\mu\nu}\chi F_{\mu\nu}$ (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]

$\displaystyle\mathcal{L}_{\text{eff}}=\frac{2\mathcal{D}_{\text{magn}}e}{q^{2}}\sum_{N=p,n}\left[Q_{N}\left(m_{N}\mathcal{O}_{5}^{(N)}-\frac{q^{2}}{4m_{\chi}}\mathcal{O}_{1}^{(N)}\right)+\widetilde{\mu}_{N}\left(m_{N}\mathcal{O}_{6}^{(N)}-\frac{q^{2}}{m_{N}}\mathcal{O}_{4}^{(N)}\right)\right]\,.$
(3.3)

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_{\text{R}}$ at small recoil energies, implying a steeper
recoil spectrum. For our benchmark scenario, we choose
$\mathcal{D}_{\text{magn}}=9.3\cdot 10^{-7}\,e\,\text{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}\rightarrow 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
$\mathcal{D}_{\text{DMDM}}$ (which has units GeV-3), we obtain

$\displaystyle\mathcal{L}_{\text{eff}}=2\mathcal{D}_{\text{DMDM}}\sum_{N=p,n}\left[Q_{N}\left(m_{N}\mathcal{O}_{5}^{(N)}-\frac{q^{2}}{4m_{\chi}}\mathcal{O}_{1}^{(N)}\right)+\widetilde{\mu}_{N}\left(m_{N}\mathcal{O}_{6}^{(N)}-\frac{q^{2}}{m_{N}}\mathcal{O}_{4}^{(N)}\right)\right]\,.$
(3.4)

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 spin-independent or spin-dependent scenarios. For
the generation of the mock data corresponding to this benchmark model, we
assume $\mathcal{D}_{\text{DMDM}}=1.6\cdot 10^{-4}\,\text{GeV}^{-3}$.

###  3.2 Mapping to $v_{\text{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_{\text{min}}$-space [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 $\mathrm{d}R/\mathrm{d}E_{\mathrm{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(v_{\text{min,0}})$ for $v_{\text{min,0}}=v_{\text{min}}(E_{0})$. For
example, for SI interactions with $f_{n}/f_{p}=1$, one finds

$\frac{\text{d}\sigma}{\text{d}E_{\text{R}}}=A^{2}\,F^{2}(E_{\text{R}})\frac{m_{T}\,\sigma_{n}}{2\,\mu_{n\chi}^{2}\,v^{2}}\;,$
(3.5)

where $A$ is the mass number of the nucleus, $F(E_{\text{R}})$ is the standard
SI form factor and $\sigma_{n}$ is the DM-neutron scattering cross section at
zero momentum transfer. Substituting this expression into equation (2.5), one
finds

$\frac{\text{d}R}{\text{d}E_{\text{R}}}=\frac{\rho\,A^{2}\,F^{2}(E_{\text{R}})\,\sigma_{n}}{2\,\mu_{n\chi}^{2}\,m_{\chi}}g(v_{\text{min}})\;.$
(3.6)

To simplify this expression, we absorb all quantities that do not depend on
the target material into the _rescaled velocity integral_

$\tilde{g}(v_{\text{min}})\equiv\frac{\rho\,\sigma_{n}}{m_{\chi}}g(v_{\text{min}})\;.$
(3.7)

The rescaled velocity integral can now be calculated from the measured
differential event rate by inverting eq. (3.6):

$\tilde{g}(v_{\text{min,0}})=\frac{2\,\mu_{n\chi}^{2}}{A^{2}\,F^{2}(E_{\text{R}})}\frac{\text{d}R}{\text{d}E_{\text{R}}}(E_{0})\;.$
(3.8)

This translation can be performed for further measurements of the differential
event rate at different recoil energies and for measurements made by
additional experiments.666We emphasise that in practice both the measurement
of $\mathrm{d}R/\mathrm{d}E_{\mathrm{R}}$ and of the corresponding recoil
energy $E_{0}$ will have some experimental error due to limited statistics and
finite energy resolution of the detector. Although we do not discuss these
complications in detail here, they are taken into account in the translations
that we show. For further details, we refer to [37]. We furthermore note that
the mapping to $v_{\text{min}}$-space becomes more complicated if the
differential cross section has a non-trivial dependence on the DM velocity
[39]. To avoid this complication, we only consider the mapping for SI
interactions. Since all experiments probe the same velocity integral, it
should be possible to fit the inferred values for $\tilde{g}(v_{\text{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.

Figure 1: Reconstruction of the velocity integral in comparison with the
underlying mock data mapped to $v_{\text{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.

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_{\chi}=50\>\text{GeV}$, to generate the mock data and to
perform the mapping to $v_{\text{min}}$-space. In the absence of Poisson
fluctuations, one can see that the inferred values of
$\tilde{g}(v_{\text{min}})$ agree perfectly with the SHM used to generate the
mock data (indicated by the red dashed line).777To study the case with no
Poisson fluctuations we employ our definition of the likelihood in eq. (2.19)
also for non-integer values of $N_{i}$. The best-fit velocity integral for
$N_{\mathrm{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_{\chi}=500\>\text{GeV}$, but still use SI interactions with
$m_{\chi}=50\>\text{GeV}$ for the mapping to $v_{\text{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.

Figure 2: Same as figure 1, but using the combined information from a xenon
experiment and an iodine experiment for the reconstruction.

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_{\text{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 of $\tilde{g}(v_{\text{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_{\text{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.

###  3.3 Goodness-of-fit estimates

For each of the models discussed above, and given a value of the DM mass
$m_{\chi}$ 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 $\mathcal{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\mathcal{L}$). Ideally, the resulting minimum value of
$-2\log\mathcal{L}$ (called $x_{0}$) could simply be compared to a $\chi^{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 $\chi^{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\mathcal{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 $\zeta(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

$p=\int_{x_{0}}^{\infty}\zeta(x_{1})\mathrm{d}x_{1}\;.$ (3.9)

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. 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 best-fit 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. 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. 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.

Figure 3: Illustration of our definition of the similarity $S$ (left panel)
and the distinguishability $D$ (right panel). See the text for details.

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.888We note that it is still possible that
the distributions of $x_{0}$ and $x_{1}$ would look different for a different
test statistic. In particular, using an unbinned instead of a binned
likelihood would retain additional information which may lead to stronger
discrimination power. The use of such unbinned likelihoods is however
computationally very expensive and therefore left for future work. 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\approx 50\%$ and $D\approx 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 $\zeta(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%).

Figure 4: Distribution of $p$-values obtained directly via Monte Carlo
simulations, using for concreteness the true model DMDM and the fitted model
SI. As indicated in the figure and explained in the text, in this example the
similarity $S$ and the distinguishability $D$ obtained via our approximate
method using a _typical prediction_ of the fitted model is in good agreement
with the direct evaluation of these quantities from the distribution of
$p$-values.

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.

##  4 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.999In order to achieve
background-free data, a detector based on liquid xenon typically needs to
reject nuclear recoils that look too similar to electron recoils, leading to a
reduction in the total exposure. We do not consider this effect, nor the
potential background from solar neutrinos. 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.

Target | Exposure [kg yr] | Energy range [keV] | Number of bins | Energy resolution [keV]
---|---|---|---|---
Xe | 2000 | 4–50 | 7 | 0.6 $\sqrt{E_{\mathrm{R}}/1\>\text{keV}}$
Ge | 200 | 0.3–50 | 9 | 0.06
I | 100 | 10–100 | 9 | 0.15 $\sqrt{E_{\mathrm{R}}/1\>\text{keV}}$
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.

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.

Considering only a xenon experiment

Figure 5: Distribution of $x_{0}$ (blue) and $x_{1}$ (red) for different
combinations of true models (rows) and fitted models (columns), using only the
information from a future xenon-based experiment. Following section 3.1, we
consider spin-independent (SI) and spin-dependent (SD) interactions, as well
as DM with an anapole moment (AM), a magnetic dipole moment (MDM), or a dark
magnetic dipole moment (DMDM). A pair of models can be distinguished in a
halo-independent way if the distributions of $x_{0}$ and $x_{1}$ are
sufficiently different, leading to a small (large) value of the similarity $S$
(the distinguishability $D$). We indicate in green the scenarios in which
$S<0.05$ or $D>0.5$, corresponding to an exclusion of the fitted model in a
typical realisation of the true model with at least $95\%$ C.L.

###  4.1 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 [69, 27, 29].

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).

Figure 6: Comparison of the predicted event rates and the best-fit event rates
for two different combinations of true model and fitted model. For the
predictions, we show error bars of size $\sqrt{N}$, where $N$ is the predicted
number of events. In these plots no Poisson fluctuations are included, so
prediction and best-fit model should agree perfectly if the true model were
identical to the fitted model.

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.

###  4.2 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^{\prime}$
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}\sim
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}\approx-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$.

Figure 7: Left: Same as figure 5 (last row, first column), but allowing
varying $f_{n}/f_{p}$. Right: Same as figure 6 (right), but allowing varying
$f_{n}/f_{p}$.

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}\approx-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.

##  5 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.101010We note that, furthermore, considering more than one
experiment tremendously improves the ability to constrain the DM mass in a
halo-independent way, see e.g. [22, 24]. Although we do not consider the
reconstruction of the DM mass in the present context, we emphasise that this
reason alone makes the combination of data from several experiments a very
important undertaking. A trivial example would be an argon-based 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 low-background iodine (I) crystal scintillator such
as SABRE [7].111111As discussed above, an argon-based experiment would
trivially allow to distinguish between SI and SD interactions. Due to its
large threshold however it does not improve the distinguishability for cases
with non-standard momentum and velocity dependence and is therefore of limited
interest in the present context. Similarly, we also do not consider fluorine
experiments, since the spectral information is of crucial information in order
to distinguish between different models and it is very difficult for bubble-
chamber experiments to infer the nuclear recoil energy. 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.

###  5.1 Combination of xenon and germanium

Combination of xenon and germanium

Figure 8: Same as figure 5, but considering data from both a xenon and a
germanium target. Again, a combination of a true and fitted model can be
distinguished in a halo-independent way if the distributions of $x_{0}$ (blue)
and $x_{1}$ (red) are sufficiently different, quantified by a small similarity
$S$ and/or a large distinguishability $D$.

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
$\chi^{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
xenon-based 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].

Assuming that for a given DM mass the two experiments based on xenon and
germanium are sensitive to a similar range of $v_{\text{min}}$ values, the
rates expected for SI couplings can schematically be written as

$\displaystyle\text{Rate}_{\text{SI,Xe}}\sim(\text{SI
sensitivity})_{\text{Xe}}\,$ $\displaystyle\times\,g(v_{\text{min}})\,,$
$\displaystyle\text{Rate}_{\text{SI,Ge}}\sim(\text{SI
sensitivity})_{\text{Ge}}\,$ $\displaystyle\times\,g(v_{\text{min}})\,.$ (5.1)

Similarly, the rates induced by SD couplings are given by

$\displaystyle\text{Rate}_{\text{SD,Xe}}\sim(\text{SD
sensitivity})_{\text{Xe}}\,$ $\displaystyle\times\,g(v_{\text{min}})\,,$
$\displaystyle\text{Rate}_{\text{SD,Ge}}\sim(\text{SD
sensitivity})_{\text{Ge}}\,$ $\displaystyle\times\,g(v_{\text{min}})\,.$ (5.2)

In a region of $v_{\text{min}}$-space accessible to both experiments, we can
define the ratio of rates for SD interactions as

$\mathrm{RR}_{\mathrm{SD}}(\text{Xe}/\text{Ge})\equiv\frac{\text{Rate}_{\text{SD,Xe}}}{\text{Rate}_{\text{SD,Ge}}}=\frac{(\text{SD
sensitivity})_{\text{Xe}}}{(\text{SD sensitivity})_{\text{Ge}}}$ (5.3)

and analogously $\mathrm{RR}_{\mathrm{SI}}(\text{Xe}/\text{Ge})$ for SI
interactions.

In order for SI and SD interactions to be distinguishable in a halo-
independent way, it is therefore necessary that
$\mathrm{RR}_{\mathrm{SI}}(\text{Xe}/\text{Ge})$ is sufficiently different
from $\mathrm{RR}_{\mathrm{SD}}(\text{Xe}/\text{Ge})$, so that only one of
these scenarios can give a good fit to the data for the assumed DM
mass.121212Note that even if this is the case, the second scenario may still
give a good fit to the data for a different choice of the DM mass, which can
change both the ratio of rates and the overlap in $v_{\text{min}}$-space. In
other words, having different ratios of rates is a necessary but not
sufficient condition for two interactions to be distinguishable in a halo-
independent way. It turns out, however, that numerically the ratios of the
sensitivities are, by coincidence, very similar:

$\displaystyle\mathrm{RR}_{\mathrm{SI}}(\text{Xe}/\text{Ge})$
$\displaystyle\simeq\frac{A_{\text{Xe}}^{2}\mu_{\text{Xe}}^{2}}{A_{\text{Ge}}^{2}\mu_{\text{Ge}}^{2}}\approx\
3.2\frac{\mu_{\text{Xe}}^{2}}{\mu_{\text{Ge}}^{2}}$ (5.4)
$\displaystyle\mathrm{RR}_{\mathrm{SD}}(\text{Xe}/\text{Ge})$
$\displaystyle\simeq\frac{\sum_{T\in\left\\{129\text{Xe},131\text{Xe}\right\\}}\eta_{T}(S_{p}^{(T)}-S_{n}^{(T)})^{2}(J_{T}+1)/J_{T}}{\eta_{73\text{Ge}}(S_{p}^{(73\text{Ge})}-S_{n}^{(73\text{Ge})})^{2}(J_{73\text{Ge}}+1)/J_{73\text{Ge}}}\approx
2.9\frac{\mu_{\text{Xe}}^{2}}{\mu_{\text{Ge}}^{2}}\,,$ (5.5)

where $\eta_{T}$ is the natural abundance of the target isotope $T$,
$S_{p}^{(T)}$ and $S_{n}^{(T)}$ 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

$\mathrm{RR}_{\mathrm{SI}}(\text{Xe}/\text{Ge})\approx
1.1\,\mathrm{RR}_{\mathrm{SD}}(\text{Xe}/\text{Ge})\;.$ (5.6)

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_{\text{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_{\perp}^{2}$ term, with a dependence
on the target nucleus similar to the SI case. Hence,
$\mathrm{RR}_{\mathrm{AM}}(\text{Xe}/\text{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_{\text{min}}^{4}\,\mu_{T\chi}^{4}$, so that scattering is
enhances for heavy targets:
$\mathrm{RR}_{\mathrm{DMDM}}(\text{Xe}/\text{Ge})\gg\mathrm{RR}_{\mathrm{SI}}(\text{Xe}/\text{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_{\chi}\approx 10\>\text{GeV}$ (compared to an assumed true DM mass of
$50\>\text{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.

###  5.2 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).

Figure 9: Left: The distributions of $x_{0}$ and $x_{1}$ for the case that the
true model corresponds to DMDM interactions and the fitted model is SI with
variable coupling ratio $f_{n}/f_{p}$, using the combined information from a
xenon and a germanium experiment. Right: Illustration of the best-fit velocity
integral obtained for $f_{n}/f_{p}=-0.7$. The mapping to
$v_{\text{min}}$-space is described in section 3.2 (see also figures 1 and 2).

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-to-proton coupling ratio for which
the recoil spectrum expected from SI interactions is non-monotonic _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 halo-
independent 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_{\text{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_{\text{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_{\text{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.

###  5.3 Combination of xenon and iodine

Combination of xenon and iodine

Figure 10: Same as figures 5 and 8, but considering data from both a xenon and
a iodine target.

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\times 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

$\displaystyle\mathrm{RR}_{\mathrm{SI}}(\text{Xe}/\text{I})\approx
1.1\frac{\mu_{\text{Xe}}^{2}}{\mu_{\text{I}}^{2}}\quad\text{and}\quad\mathrm{RR}_{\mathrm{SD}}(\text{Xe}/\text{I})\approx
1.2\frac{\mu_{\text{Xe}}^{2}}{\mu_{\text{I}}^{2}}\,,$ (5.7)

so that

$\mathrm{RR}_{\mathrm{SI}}(\text{Xe}/\text{I})\approx
0.9\,\mathrm{RR}_{\mathrm{SD}}(\text{Xe}/\text{I})\;.$ (5.8)

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 $\propto 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.

##  6 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 spin-
dependent 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 halo-independent 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 spin-independent 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 iodine-based
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.

## Acknowledgments

We thank Julien Billard, David Cerdeno, Brian Feldstein, Mattia Fornasa,
Alejandro Ibarra, Miguel Peiro and Andreas Rappelt for discussions. This work
is supported by the German Science Foundation (DFG) under the Collaborative
Research Center (SFB) 676 Particles, Strings and the Early Universe, as well
as by the Studienstiftung des Deutschen Volkes and the TUM Graduate School.

##  Appendix 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

$\displaystyle\mathcal{L}_{\text{eff}}=\sum_{\lambda=1}^{15}\sum_{\tau=0}^{1}c_{\lambda}^{(\tau)}\mathcal{O}_{\lambda}^{(\tau)}\;,$
(A.1)

where $\tau=0$ for isoscalar interactions (i.e. equal couplings to protons and
neutrons) and $\tau=1$ for isovector interactions (opposite couplings to
protons and neutrons).131313The operator $\mathcal{O}_{2}$ is typically
discarded, because it does not arise at leading order from a relativistic UV
completion.

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_{\lambda}^{(\tau)}$ from equation (A.1) as

$\displaystyle
c_{\lambda}^{(\tau)}\equiv\frac{m_{N}^{2}}{q^{2}}\,z_{0,\lambda}^{(\tau)}+z_{1,\lambda}^{(\tau)}+\frac{q^{2}}{m_{N}^{2}}\,z_{2,\lambda}^{(\tau)}\,,$
(A.2)

where $m_{N}\equiv m_{p}\approx m_{n}$ is the nucleon mass. The interactions
of DM with nucleons can therefore be fully described by the $14\times 2\times
3$ coefficients $z_{\kappa,\lambda}^{(\tau)}$, which we combine into a
coefficient vector $\mathbf{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_{\text{med}}^{2}\sim q^{2}$.

In order to calculate the differential scattering cross section
$\text{d}\sigma/\text{d}E_{\mathrm{R}}$ for a given target nucleus, the next
step is to calculate the so-called DM response functions
$S_{\alpha}^{(\tau,\tau^{\prime})}$ for $\alpha=1,\ldots,8$, which depend on
the coefficient vector $\mathbf{z}$ in a known way. The DM response functions
then need to be multiplied with the corresponding nuclear form factors
$\widetilde{W}_{\alpha}^{(\tau,\tau^{\prime})}$:

$\frac{\text{d}\sigma}{\text{d}E_{\text{R}}}=\frac{m_{T}}{2\pi
v^{2}}\sum_{\alpha=1}^{8}\sum_{\tau=0}^{1}\sum_{\tau^{\prime}=0}^{1}\,S_{\alpha}^{(\tau,\tau^{\prime})}\left(\textbf{z},v^{2},q^{2}\right)\,\,\widetilde{W}_{\alpha}^{(\tau,\tau^{\prime})}(q^{2})\,,$
(A.3)

where the sum over
$S_{\alpha}^{(\tau,\tau^{\prime})}\,W_{\alpha}^{(\tau,\tau^{\prime})}$ 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
$\boldsymbol{S_{\alpha}^{(\tau,\tau^{\prime})}}$

The DM response functions $S_{\alpha}^{(\tau,\tau^{\prime})}$ are defined in
appendix A of [57]. The index $\alpha=1,\dots,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}$:

$\displaystyle
S_{\alpha}^{(\tau,\tau^{\prime})}\left(\textbf{z},v^{2},q^{2}\right)\equiv
S_{1,\alpha}^{(\tau,\tau^{\prime})}\left(\textbf{z},q^{2}\right)+v^{2}\cdot
S_{2,\alpha}^{(\tau,\tau^{\prime})}\left(\textbf{z},q^{2}\right)\,.$ (A.4)

The functions $S_{\beta,\alpha}^{(\tau,\tau^{\prime})}\left(q^{2}\right)$
(with $\beta\in\left\\{1,2\right\\}$) can furthermore be written as a linear
combination of integer powers of $q^{2}$ (see [57]):

$\displaystyle
S_{\beta,\alpha}^{(\tau,\tau^{\prime})}\left(\textbf{z},q^{2}\right)=\sum_{l=-2}^{5}\sigma_{\beta,\alpha,l}^{(\tau,\tau^{\prime})}\left(\textbf{z}\right)\cdot
a^{l}\,,$ (A.5)

where we have introduced the dimensionless quantity

$\displaystyle a\equiv\frac{q^{2}}{m_{N}^{2}}$ (A.6)

and implicitly defined the coefficients
$\sigma_{\beta,\alpha,l}^{(\tau,\tau^{\prime})}$ for $l=-2,\ldots,5$. These
$2\times 8\times 8\times 2\times 2$ coefficients fully characterise the DM
response and can be calculated in terms of the coefficient vector
$\mathbf{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
$\sigma_{\beta,\alpha,l}^{(\tau,\tau^{\prime})}(\mathbf{z})$ as the quantity
describing the particle physics properties of DM.

###  A.2 Nuclear form factors
$\boldsymbol{\widetilde{W}_{\alpha}^{(\tau,\tau^{\prime})}}$

Apart from the DM response functions, the differential cross section defined
in equation (A.3) also depends on the nuclear form factors
$\widetilde{W}_{\alpha}^{(\tau,\tau^{\prime})}$, which are provided in
[50].141414For convenience we define
$\widetilde{W}_{\alpha}^{(\tau,\tau^{\prime})}\equiv\frac{4\pi}{2J_{A}+1}W_{\alpha}^{(\tau,\tau^{\prime})}$,
where $W_{\alpha}^{(\tau,\tau^{\prime})}$ are the functions provided in [50].
These form factors depend on $q^{2}$ via the combination
$y=\frac{1}{4}b_{A}^{2}q^{2}$, with

$b_{A}=1\>\text{fm}\times\sqrt{41.467/\left(45A^{-1/3}-25A^{-2/3}\right)}$
(A.7)

The nuclear form factors are given as fit functions of the form

$\displaystyle\widetilde{W}_{\alpha}^{(\tau,\tau^{\prime})}(y)=\frac{4\pi}{2J_{A}+1}\,\exp(-2y)\,\sum_{k=0}^{10}\overline{w}_{\alpha,k}^{(\tau,\tau^{\prime})}\cdot
y^{k}\,.$ (A.8)

Defining $y(a)=\frac{1}{4}b_{A}^{2}m_{N}^{2}\cdot a$, this can be rewritten in
the form

$\displaystyle\widetilde{W}_{\alpha}^{(\tau,\tau^{\prime})}(a)=\exp(-2y(a))\,\sum_{k=0}^{10}w_{\alpha,k}^{(\tau,\tau^{\prime})}\cdot
a^{k}\,,$ (A.9)

with
$w_{\alpha,k}^{(\tau,\tau^{\prime})}=\frac{4\pi}{2J_{A}+1}\left(\frac{1}{4}b_{A}^{2}m_{N}^{2}\right)^{k}\overline{w}_{\alpha,k}^{(\tau,\tau^{\prime})}$.

###  A.3 Decomposing the differential cross section

Substituting equations (A.4), (A.5), and (A.9) into equation (A.3), it becomes
clear that $\text{d}\sigma/\text{d}E_{\text{R}}$ can be decomposed into

$\frac{\text{d}\sigma}{\text{d}E_{\text{R}}}=\frac{\text{d}\sigma_{1}}{\text{d}E_{\text{R}}}\frac{1}{v^{2}}+\frac{\text{d}\sigma_{2}}{\text{d}E_{\text{R}}}$
(A.10)

with

$\displaystyle\frac{\text{d}\sigma_{\beta}}{\text{d}E_{\text{R}}}$
$\displaystyle=\frac{m_{T}}{2\pi}\sum_{\alpha=1}^{8}\,\sum_{\tau=0}^{1}\,\sum_{\tau^{\prime}=0}^{1}\,S_{\beta,\alpha}^{(\tau,\tau^{\prime})}\left(\textbf{z},q^{2}\right)e^{-2y(a)}\sum_{k=0}^{10}w_{\alpha,k}^{(\tau,\tau^{\prime})}\cdot
a^{k}$
$\displaystyle=\frac{m_{T}}{2\pi}\sum_{\alpha=1}^{8}\,\sum_{\tau=0}^{1}\,\sum_{\tau^{\prime}=0}^{1}\,\sum_{k=0}^{10}\,\sum_{l=-2}^{5}\sigma_{\beta,\alpha,l}^{(\tau,\tau^{\prime})}\left(\textbf{z}\right)\cdot
w_{\alpha,k}^{(\tau,\tau^{\prime})}\cdot e^{-2y(a)}a^{k+l}\;.$ (A.11)

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

$\displaystyle G_{ij}$
$\displaystyle=\frac{\kappa\,\rho}{2\,m_{T}\,m_{\chi}}\int_{E_{j}}^{E_{j+1}}\frac{\text{d}\sigma_{1}}{\mathrm{d}E_{\text{R}}}\,\epsilon(E_{\mathrm{R}})\left[\text{erf}\left(\frac{E_{i+1}-E_{\text{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)-\text{erf}\left(\frac{E_{i}-E_{\text{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)\right]\,\mathrm{d}E_{\mathrm{R}}\;,$ (A.12)
$\displaystyle H_{ij}$
$\displaystyle=\frac{\kappa\,\rho}{2\,m_{T}\,m_{\chi}}\int_{E_{j}}^{E_{j+1}}\frac{\text{d}\sigma_{2}}{\mathrm{d}E_{\text{R}}}\,\epsilon(E_{\mathrm{R}})\left[\text{erf}\left(\frac{E_{i+1}-E_{\text{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)-\text{erf}\left(\frac{E_{i}-E_{\text{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)\right]\,\mathrm{d}E_{\mathrm{R}}\;.$ (A.13)

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_{\mathrm{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
$\sigma_{\beta,\alpha,l}^{(\tau,\tau^{\prime})}\left(\textbf{z}\right)$ and
$w_{\alpha,k}^{(\tau,\tau^{\prime})}$ do not depend on the recoil energy, it
becomes clear that we first need to calculate the functions

$\displaystyle d_{i}^{(p)}(a)$
$\displaystyle\equiv\frac{1}{2}\int_{0}^{E_{\text{R}}(a)}\text{d}E_{\mathrm{R}}\,e^{-2\,y(a(E_{\mathrm{R}}))}\,a(E_{\mathrm{R}})^{p}\,\epsilon(E_{\mathrm{R}})\left[\text{erf}\left(\frac{E^{\prime}_{i+1}-E_{\mathrm{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)-\text{erf}\left(\frac{E^{\prime}_{i}-E_{\mathrm{R}}}{\sqrt{2}\Delta
E_{\text{R}}}\right)\right]$
$\displaystyle=\frac{m_{N}^{2}}{4\,m_{T}}\int_{0}^{a}\text{d}a^{\prime}\,e^{-2y(a^{\prime})}\,a^{\prime
p}\,\epsilon\Bigl{(}\frac{a^{\prime}\,m_{N}^{2}}{2\,m_{T}}\Bigr{)}\Biggl{[}\text{erf}\biggl{(}\frac{E^{\prime}_{i+1}-\frac{a^{\prime}\,m_{N}^{2}}{2\,m_{T}}}{\sqrt{2}\Delta
E_{\text{R}}}\biggr{)}-\text{erf}\biggl{(}\frac{E^{\prime}_{i}-\frac{a^{\prime}\,m_{N}^{2}}{2\,m_{T}}}{\sqrt{2}\Delta
E_{\text{R}}}\biggr{)}\Biggr{]}\,,$ (A.14)

where $a=q^{2}/m_{N}^{2}=2\,m_{T}\,E_{\text{R}}/m_{N}^{2}$.

These functions can be thought of as the predicted number of events in the
$i$th bin for a hypothetical scattering cross section with simplified energy
dependence of the form
$\mathrm{d}\sigma/\mathrm{d}E_{\mathrm{R}}\sim\exp(-E_{\mathrm{R}})E_{\mathrm{R}}^{p}$
and a velocity integral of the form
$g(v_{\text{min}})=\Theta(a-a(v_{\text{min}}))$, where
$a(v_{\text{min}})=4\mu_{T\chi}^{2}v_{\text{min}}^{2}/m_{N}^{2}$. 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}^{(p)}(a)$
need to be calculated for $p\in\left\\{-2,-1,\dots,15\right\\}$. Crucially,
they are independent of the DM mass, the coefficient vector $\mathbf{z}$ and
of the binning of $v_{\text{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_{\text{min}}$-space.

Once the particle physics properties of DM have been specified in terms of the
coefficients $\sigma_{\beta,\alpha,l}^{(\tau,\tau^{\prime})}$ introduced
above, it is possible to combine the reduced event rates $d_{i}^{(p)}(a)$ into
two function for each bin:

$\displaystyle G_{i}(a)$
$\displaystyle=\frac{\kappa\,\rho}{2\pi\,m_{\chi}}\sum_{\alpha=1}^{8}\,\sum_{\tau=0}^{1}\,\sum_{\tau^{\prime}=0}^{1}\,\sum_{l=-2}^{5}\,\sum_{k=0}^{10}\,\sigma_{1,\alpha,l}^{(\tau,\tau^{\prime})}\left(\textbf{z}\right)\,w_{\alpha,k}^{(\tau,\tau^{\prime})}\,\cdot\,d_{i}^{(l+k)}(a)\,,$
$\displaystyle H_{i}(a)$
$\displaystyle=\frac{\kappa\,\rho}{2\pi\,m_{\chi}}\sum_{\alpha=1}^{8}\,\sum_{\tau=0}^{1}\,\sum_{\tau^{\prime}=0}^{1}\,\sum_{l=-2}^{5}\,\sum_{k=0}^{10}\,\sigma_{2,\alpha,l}^{(\tau,\tau^{\prime})}\left(\textbf{z}\right)\,w_{\alpha,k}^{(\tau,\tau^{\prime})}\,\cdot\,d_{i}^{(l+k)}(a)\,,$
(A.15)

where the coefficients $w_{\alpha,k}^{(\tau,\tau^{\prime})}$ contain the
details of the nuclear response (see above).

Finally, once the binning $v_{1},\ldots,v_{N_{\mathrm{s}}}$ of
$v_{\text{min}}$-space has been specified, one simply needs to calculate the
corresponding values $a_{j}\equiv a(v_{j})$. The matrix $G_{ij}$ is then given
by151515Note that for a target consisting of different isotopes (or in fact
different elements), the matrix $G_{ij}$ needs to be calculated for each
isotope. The different matrices can then simply be multiplied with the mass
fraction of each isotope and summed up.

$\displaystyle G_{ij}=G_{i}\left(a_{j+1}\right)-G_{i}\left(a_{j}\right)$
(A.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.

##  Appendix B Calculating event rates from the velocity integral

###  B.1 The second velocity integral

We consider the case that the velocity integral $g(v_{\text{min}})$ is
piecewise constant, so that it can be written as

$g(v_{\text{min}})=\sum_{j=1}^{N_{\mathrm{s}}}l_{j}\,\Theta(v_{j+1}-v_{\text{min}})\;,$
(B.1)

where $\Theta(x)$ is the Heaviside step function and $l_{j}\equiv
g_{j}-g_{j+1}$ (with $g_{N_{\mathrm{s}}+1}=0$). In order to calculate
$h(v_{\text{min}})$ for given $v_{\text{min}}$, we choose $k$ in such a way
that $v_{k-1}<v_{\text{min}}<v_{k}$. One then obtains:

$\displaystyle h(v_{\text{min}})$
$\displaystyle=\left[-g(v)\,v^{2}\right]_{v_{\text{min}}}^{\infty}+\int_{v_{\text{min}}}^{\infty}2\,v\,g(v)\mathrm{d}v$
$\displaystyle=\sum_{j=1}^{N_{\mathrm{s}}}l_{j}\left[-\Theta(v_{j+1}-v)\,v^{2}\right]_{v_{\text{min}}}^{\infty}+\sum_{j=1}^{N_{\mathrm{s}}}l_{j}\int_{v_{\text{min}}}^{\infty}2\,v\,\Theta(v_{j+1}-v)\mathrm{d}v$
$\displaystyle=\sum_{j=k-1}^{N_{\mathrm{s}}}l_{j}v_{\text{min}}^{2}+\sum_{j=k-1}^{N_{\mathrm{s}}}l_{j}\int_{v_{\text{min}}}^{v_{j+1}}2v\,\mathrm{d}v$
$\displaystyle=\sum_{j=k-1}^{N_{\mathrm{s}}}l_{j}v_{\text{min}}^{2}+\sum_{j=k-1}^{N_{\mathrm{s}}}l_{j}\left(v_{j+1}^{2}-v_{\text{min}}^{2}\right)$
$\displaystyle=\sum_{j=k-1}^{N_{\mathrm{s}}}l_{j}\,v_{j+1}^{2}\;.$ (B.2)

We observe that $h(v_{\text{min}})$ is also a piecewise-constant function,
i.e. we can write $h(v_{\text{min}})=h_{j}$ for
$v_{\text{min}}\in\left[v_{j},\,v_{j+1}\right]$ with

$h_{j}=\sum_{j^{\prime}=j}^{N_{\mathrm{s}}}l_{j^{\prime}}\,v_{j^{\prime}+1}^{2}=\sum_{j^{\prime}=j}^{N_{\mathrm{s}}}(g_{j^{\prime}}-g_{j^{\prime}+1})\,v_{j^{\prime}+1}^{2}\;.$
(B.3)

The relationship between $h_{j}$ and $g_{j^{\prime}}$ can be written as
$h_{j}=F_{jj^{\prime}}g_{j^{\prime}}$ with $F_{jj}=v_{j+1}^{2}$,
$F_{jj^{\prime}}=0$ for $j>j^{\prime}$ and
$F_{jj^{\prime}}=v_{j^{\prime}+1}^{2}-v_{j^{\prime}}^{2}$ for $j^{\prime}>j$.

Using the matrix $F_{jj^{\prime}}$ 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:

$R_{i}=\sum D_{ij}\,g_{j}\quad\text{with}\quad
D_{ij}=G_{ij}+\sum_{k}F_{kj}H_{ik}\;.$ (B.4)

###  B.2 The impact of high-velocity bins

Let us now discuss how to determine the range of $v_{\text{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 $j$th bin in $v_{\text{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 $\Delta E_{\mathrm{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_{\text{min}}$ and a
$v_{\text{max}}$ such that $G_{ij}\approx 0$ and $H_{ij}\approx 0$ for
$v_{j}<v_{\text{min}}$ or $v_{j}>v_{\text{max}}$. This observation suggests
that one can simply take $v_{1}=v_{\text{min}}$ and
$v_{N_{\mathrm{s}}+1}=v_{\text{max}}$.

However, things become more complicated when the matrices $G_{ij}$ and
$H_{ij}$ are combined into the matrix
$D_{ij}=G_{ij}+\sum_{j^{\prime}}H_{ij^{\prime}}F_{j^{\prime}j}$. To illustrate
the problem, let us consider a range of $v_{\text{min}}$-space
$[v_{1},v_{N_{\mathrm{s}}+1}]$ such that
$E_{\mathrm{R}}(v_{N_{\mathrm{s}}+1})$ is significantly larger than
$v_{\text{max}}$. In this case there is a $j_{\text{max}}$ such that for
$J>j_{\text{max}}$ one obtains $G_{iJ}=0$ and $H_{iJ}=0$ for all bins $i$. For
such a $J$ it then follows that

$D_{iJ}=\sum_{j^{\prime}=1}^{N_{\mathrm{s}}}H_{ij^{\prime}}F_{j^{\prime}J}=\sum_{j^{\prime}=1}^{j_{\text{max}}}H_{ij^{\prime}}F_{j^{\prime}J}=(v_{J+1}^{2}-v_{J}^{2})\sum_{j^{\prime}=1}^{j_{\text{max}}}H_{ij^{\prime}}\;,$
(B.5)

because $F_{j^{\prime}j}=v_{j+1}^{2}-v_{j}^{2}$ independent of $j^{\prime}$
for $j>j^{\prime}$. Crucially we find that $D_{iJ}$ can be non-zero even for
$J>j_{\text{max}}$.

The calculation above shows that arbitrarily high bins in
$v_{\text{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_{\text{min}})$ and $h(v_{\text{min}})$, which imply that a local
variation in $g(v_{\text{min}})$ leads to a global variation in
$h(v_{\text{min}})$, see also [39]. To be more specific, if we vary
$g(v_{\text{min}})$ only in the interval $\left[v_{j},v_{j+1}\right]$ and keep
it fixed everywhere else, $h(v_{\text{min}})$ will vary for all velocities
$v_{\text{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 $i$th row of $H_{ij}$, i.e.
$H_{i}\equiv\sum_{j}H_{ij}$, we obtain

$D_{iJ}\equiv(v_{J+1}^{2}-v_{J}^{2})H_{i}\;.$ (B.6)

In other words, for fixed $J>j_{\text{max}}$ the column vector $D_{iJ}$ is
simply proportional to the vector $H_{i}$.

As a result, for any vector
$(g_{j})=\left(g_{1},\ldots,g_{N_{\mathrm{s}}}\right)$ in
$v_{\text{min}}$-space we can define a new vector
$(\tilde{g}_{j})=\left(g_{1},\ldots,g_{j_{\text{max}}},\tilde{g},0,\ldots,0\right)$
with

$\tilde{g}=\frac{1}{(v_{j_{\text{max}}+2}^{2}-v_{j_{\text{max}}+1}^{2})}\sum_{J=j_{\text{max}}+1}^{N_{\mathrm{s}}}(v_{J+1}^{2}-v_{J}^{2})g_{J}$
(B.7)

such that $\sum_{j}D_{ij}g_{j}=\sum_{j}D_{ij}\tilde{g}_{j}$. While $(g_{j})$
is assumed to be monotonically decreasing (i.e. $g_{j+1}\leq g_{j}$), this no
longer needs to be true for $(\tilde{g}_{j})$, because it is possible that
$\tilde{g}>g_{j_{\text{max}}}$. In other words, adding a large number of
monotonically decreasing bins beyond $j_{\text{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_{\text{min}}$ in order to find the best-fit velocity integral.

## References

  * [1] P. Cushman et al., Snowmass CF1 Summary: WIMP Dark Matter Direct Detection,  1310.8327.
  * [2] L. Baudis, WIMP Dark Matter Direct-Detection Searches in Noble Gases,  Phys. Dark Univ. 4 (2014) 50–59, [1408.4371].
  * [3] LZ Collaboration, D. S. Akerib et al., LUX-ZEPLIN (LZ) Conceptual Design Report, 1509.02910.
  * [4] XENON Collaboration, E. Aprile et al., Physics reach of the XENON1T dark matter experiment,  JCAP 1604 (2016) 027, [1512.07501].
  * [5] DARWIN Collaboration, J. Aalbers et al., DARWIN: towards the ultimate dark matter detector, 1606.07001.
  * [6] PICO collaboration webpage, accessed on July 6, 2016.  http://www.picoexperiment.com/pico250.php.
  * [7] E. Shields, J. Xu, and F. Calaprice, SABRE: A New NaI(T1) Dark Matter Direct Detection Experiment,  Phys. Procedia 61 (2015) 169–178.
  * [8] SuperCDMS Collaboration, R. Calkins, The SuperCDMS Soudan high threshold WIMP search and the planned SuperCDMS SNOLAB experiment,  J. Phys. Conf. Ser. 718 (2016) 042009.
  * [9] EDELWEISS Collaboration, Q. Arnaud, J. Billard, and A. Juillard, Optimizing Cryogenic Detectors for Low-Mass WIMP Searches,  J. Low. Temp. Phys. 184 (2016) 308–315.
  * [10] LUX Collaboration, D. S. Akerib et al., Improved Limits on Scattering of Weakly Interacting Massive Particles from Reanalysis of 2013 LUX Data,  Phys. Rev. Lett. 116 (2016) 161301, [1512.03506].
  * [11] XENON100 Collaboration, E. Aprile et al., Dark Matter Results from 225 Live Days of XENON100 Data,  Phys. Rev. Lett. 109 (2012) 181301, [1207.5988].
  * [12] PICO Collaboration, C. Amole et al., Improved dark matter search results from PICO-2L Run 2,  Phys. Rev. D93 (2016) 061101, [1601.03729].
  * [13] SuperCDMS Collaboration, R. Agnese et al., Search for Low-Mass Weakly Interacting Massive Particles with SuperCDMS,  Phys. Rev. Lett. 112 (2014) 241302, [1402.7137].
  * [14] CRESST Collaboration, G. Angloher et al., Results on light dark matter particles with a low-threshold CRESST-II detector,  Eur. Phys. J. C76 (2016) 25, [1509.01515].
  * [15] M. Kuhlen, N. Weiner, J. Diemand, P. Madau, B. Moore, et al., Dark Matter Direct Detection with Non-Maxwellian Velocity Structure,  JCAP 1002 (2010) 030, [0912.2358].
  * [16] M. Lisanti, L. E. Strigari, J. G. Wacker, and R. H. Wechsler, The Dark Matter at the End of the Galaxy,  Phys. Rev. D83 (2011) 023519, [1010.4300].
  * [17] C. McCabe, The Astrophysical Uncertainties Of Dark Matter Direct Detection Experiments,  Phys. Rev. D82 (2010) 023530, [1005.0579].
  * [18] Y.-Y. Mao, L. E. Strigari, R. H. Wechsler, H.-Y. Wu, and O. Hahn, Halo-to-Halo Similarity and Scatter in the Velocity Distribution of Dark Matter,  Astrophys. J. 764 (2013) 35, [1210.2721].
  * [19] L. E. Strigari and R. Trotta, Reconstructing WIMP Properties in Direct Detection Experiments Including Galactic Dark Matter Distribution Uncertainties,  JCAP 0911 (2009) 019, [0906.5361].
  * [20] A. H. G. Peter, Getting the astrophysics and particle physics of dark matter out of next-generation direct detection experiments,  Phys. Rev. D81 (2010) 087301, [0910.4765].
  * [21] A. M. Green, Dependence of direct detection signals on the WIMP velocity distribution,  JCAP 1010 (2010) 034, [1009.0916].
  * [22] A. H. G. Peter, WIMP astronomy and particle physics with liquid-noble and cryogenic direct-detection experiments,  Phys. Rev. D83 (2011) 125029, [1103.5145].
  * [23] M. Pato, L. E. Strigari, R. Trotta, and G. Bertone, Taming astrophysical bias in direct dark matter searches,  JCAP 1302 (2013) 041, [1211.7063].
  * [24] B. J. Kavanagh and A. M. Green, Improved determination of the WIMP mass from direct detection data,  Phys. Rev. D86 (2012) 065027, [1207.2039].
  * [25] M. T. Frandsen, F. Kahlhoefer, C. McCabe, S. Sarkar, and K. Schmidt-Hoberg, The unbearable lightness of being: CDMS versus XENON,  JCAP 1307 (2013) 023, [1304.6066].
  * [26] B. J. Kavanagh and A. M. Green, Model independent determination of the dark matter mass from direct detection experiments,  Phys. Rev. Lett. 111 (2013) 031302, [1303.6868].
  * [27] A. H. G. Peter, V. Gluscevic, A. M. Green, B. J. Kavanagh, and S. K. Lee, WIMP physics with ensembles of direct-detection experiments,  Phys. Dark Univ. 5-6 (2014) 45–74, [1310.7039].
  * [28] L. E. Strigari, Galactic Searches for Dark Matter,  Phys. Rept. 531 (2013) 1–88, [1211.7090].
  * [29] B. J. Kavanagh, M. Fornasa, and A. M. Green, Probing WIMP particle physics and astrophysics with direct detection and neutrino telescope data,  Phys. Rev. D91 (2015) 103533, [1410.8051].
  * [30] B. Feldstein and F. Kahlhoefer, A new halo-independent approach to dark matter direct detection analysis,  JCAP 1408 (2014) 065, [1403.4606].
  * [31] B. Feldstein and F. Kahlhoefer, Quantifying (dis)agreement between direct detection experiments in a halo-independent way,  JCAP 1412 (2014) 052, [1409.5446].
  * [32] M. Drees and C.-L. Shan, Model-Independent Determination of the WIMP Mass from Direct Dark Matter Detection Data,  JCAP 0806 (2008) 012, [0803.4477].
  * [33] P. J. Fox, G. D. Kribs, and T. M. P. Tait, Interpreting Dark Matter Direct Detection Independently of the Local Velocity and Density Distribution,  Phys. Rev. D83 (2011) 034007, [1011.1910].
  * [34] P. J. Fox, J. Liu, and N. Weiner, Integrating Out Astrophysical Uncertainties,  Phys. Rev. D83 (2011) 103514, [1011.1915].
  * [35] C. McCabe, DAMA and CoGeNT without astrophysical uncertainties,  Phys. Rev. D84 (2011) 043525, [1107.0741].
  * [36] M. T. Frandsen, F. Kahlhoefer, C. McCabe, S. Sarkar, and K. Schmidt-Hoberg, Resolving astrophysical uncertainties in dark matter direct detection,  JCAP 1201 (2012) 024, [1111.0292].
  * [37] P. Gondolo and G. B. Gelmini, Halo independent comparison of direct dark matter detection data,  JCAP 1212 (2012) 015, [1202.6359].
  * [38] J. Herrero-Garcia, T. Schwetz, and J. Zupan, Astrophysics independent bounds on the annual modulation of dark matter signals,  Phys. Rev. Lett. 109 (2012) 141301, [1205.0134].
  * [39] E. Del Nobile, G. Gelmini, P. Gondolo, and J.-H. Huh, Generalized Halo Independent Comparison of Direct Dark Matter Detection Data,  JCAP 1310 (2013) 048, [1306.5273].
  * [40] P. J. Fox, Y. Kahn, and M. McCullough, Taking Halo-Independent Dark Matter Methods Out of the Bin,  JCAP 1410 (2014) 076, [1403.6830].
  * [41] N. Bozorgnia and T. Schwetz, What is the probability that direct detection experiments have observed Dark Matter?,  JCAP 1412 (2014) 015, [1410.6160].
  * [42] J. F. Cherry, M. T. Frandsen, and I. M. Shoemaker, Halo Independent Direct Detection of Momentum-Dependent Dark Matter,  JCAP 1410 (2014) 022, [1405.1420].
  * [43] A. J. Anderson, P. J. Fox, Y. Kahn, and M. McCullough, Halo-Independent Direct Detection Analyses Without Mass Assumptions,  JCAP 1510 (2015) 012, [1504.03333].
  * [44] G. B. Gelmini, A. Georgescu, P. Gondolo, and J.-H. Huh, Extended Maximum Likelihood Halo-independent Analysis of Dark Matter Direct Detection Data,  JCAP 1511 (2015) 038, [1507.03902].
  * [45] F. Ferrer, A. Ibarra, and S. Wild, A novel approach to derive halo-independent limits on dark matter properties,  JCAP 1509 (2015) 052, [1506.03386].
  * [46] S. Scopel, K.-H. Yoon, and J.-H. Yoon, Generalized spin-dependent WIMP-nucleus interactions and the DAMA modulation effect,  JCAP 1507 (2015), no. 07 041, [1505.01926].
  * [47] E. Del Nobile, G. B. Gelmini, and S. J. Witte, Prospects for detection of target-dependent annual modulation in direct dark matter searches,  JCAP 1602 (2016), no. 02 009, [1512.03961].
  * [48] G. B. Gelmini, J.-H. Huh, and S. J. Witte, Assessing Compatibility of Direct Detection Data: Halo-Independent Global Likelihood Analyses, 1607.02445.
  * [49] A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers, and Y. Xu, The Effective Field Theory of Dark Matter Direct Detection,  JCAP 1302 (2013) 004, [1203.3542].
  * [50] N. Anand, A. L. Fitzpatrick, and W. C. Haxton, Weakly interacting massive particle-nucleus elastic scattering response,  Phys. Rev. C89 (2014) 065501, [1308.6288].
  * [51] A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers, and Y. Xu, Model Independent Direct Detection Analyses, 1211.2818.
  * [52] R. Catena and P. Gondolo, Global fits of the dark matter-nucleon effective interactions,  JCAP 1409 (2014) 045, [1405.2637].
  * [53] M. I. Gresham and K. M. Zurek, Effect of nuclear response functions in dark matter direct detection,  Phys. Rev. D89 (2014) 123521, [1401.3739].
  * [54] R. Catena, Analysis of the theoretical bias in dark matter direct detection,  JCAP 1409 (2014) 049, [1407.0127].
  * [55] R. Catena, Prospects for direct detection of dark matter in an effective theory approach,  JCAP 1407 (2014) 055, [1406.0524].
  * [56] V. Gluscevic and A. H. G. Peter, Understanding WIMP-baryon interactions with direct detection: A Roadmap,  JCAP 1409 (2014) 040, [1406.7008].
  * [57] R. Catena and P. Gondolo, Global limits and interference patterns in dark matter direct detection,  JCAP 1508 (2015) 022, [1504.06554].
  * [58] V. Gluscevic, M. I. Gresham, S. D. McDermott, A. H. G. Peter, and K. M. Zurek, Identifying the Theory of Dark Matter with Direct Detection,  JCAP 1512 (2015) 057, [1506.04454].
  * [59] J. B. Dent, L. M. Krauss, J. L. Newstead, and S. Sabharwal, General analysis of direct dark matter detection: From microphysics to observational signatures,  Phys. Rev. D92 (2015) 063515, [1505.03117].
  * [60] R. Catena, A. Ibarra, and S. Wild, DAMA confronts null searches in the effective theory of dark matter-nucleon interactions,  JCAP 1605 (2016) 039, [1602.04074].
  * [61] M. Felizardo et al., Final Analysis and Results of the Phase II SIMPLE Dark Matter Search,  Phys. Rev. Lett. 108 (2012) 201302, [1106.3014].
  * [62] COUPP Collaboration, E. Behnke et al., First Dark Matter Search Results from a 4-kg CF3I Bubble Chamber Operated in a Deep Underground Site,  Phys. Rev. D86 (2012) 052001, [1204.3094]. [Erratum: Phys. Rev.D90,no.7,079902(2014)].
  * [63] PICASSO Collaboration, S. Archambault et al., Constraints on Low-Mass WIMP Interactions on ${}^{19}F$ from PICASSO,  Phys. Lett. B711 (2012) 153–161, [1202.1240].
  * [64] S. Chang, R. Edezhath, J. Hutchinson, and M. Luty, Leptophilic Effective WIMPs,  Phys. Rev. D90 (2014) 015011, [1402.7358].
  * [65] A. Ibarra and S. Wild, Dirac dark matter with a charged mediator: a comprehensive one-loop analysis of the direct detection phenomenology,  JCAP 1505 (2015) 047, [1503.03382].
  * [66] V. Barger, W.-Y. Keung, and D. Marfatia, Electromagnetic properties of dark matter: Dipole moments and charge form factor,  Phys. Lett. B696 (2011) 74–78, [1007.4345].
  * [67] T. Banks, J.-F. Fortin, and S. Thomas, Direct Detection of Dark Matter Electromagnetic Dipole Moments, 1007.5515.
  * [68] A. L. Fitzpatrick and K. M. Zurek, Dark Moments and the DAMA-CoGeNT Puzzle,  Phys.Rev. D82 (2010) 075004, [1007.5325].
  * [69] M. Pato, What can(not) be measured with ton-scale dark matter direct detection experiments,  JCAP 1110 (2011) 035, [1106.0743].
  * [70] M. T. Frandsen, F. Kahlhoefer, S. Sarkar, and K. Schmidt-Hoberg, Direct detection of dark matter in models with a light Z’,  JHEP 09 (2011) 128, [1107.2118].
  * [71] H. Zheng, Z. Zhang, and L.-W. Chen, Form Factor Effects in the Direct Detection of Isospin-Violating Dark Matter,  JCAP 1408 (2014) 011, [1403.5134].
  * [72] J. L. Feng, J. Kumar, D. Marfatia, and D. Sanford, Isospin-Violating Dark Matter,  Phys. Lett. B703 (2011) 124–127, [1102.4331].
  * [73] S. Chang, N. Weiner, and I. Yavin, Magnetic Inelastic Dark Matter,  Phys. Rev. D82 (2010) 125011, [1007.4200].