Exploring light mediators with low-threshold direct detection experiments
Felix Kahlhoefer Suchita Kulkarni Sebastian Wild
hep-ph astro-ph.CO hep-ex

Polarization Transfer Observables in Elastic Electron-Proton Scattering at 𝐐^2= 2.5, 5.2, 6.8 and 8.5 GeV^2
    L. Zhu
    December 30, 2023
===========================================================================================================





§ INTRODUCTION


Dark matter (DM) direct detection experiments are conventionally interpreted in terms of effective operators, which parametrize the details of the underlying interactions of the DM particle in terms of an effective suppression scale Λ <cit.>. This is a good approximation as long as the particle mediating the interaction is heavy compared to the typical momentum transfer in direct detection experiments, which is of order 1–100 MeV. Nevertheless, it is perfectly conceivable that the mediator mass is close to or below this scale, so that the effective operator description is no longer valid <cit.>. This scenario has for example been considered in the context of self-interacting DM <cit.>. In these models the exchange of light mediators induces large DM self-scattering rates <cit.>, which can potentially resolve the tension between the predictions of collisionless cold DM and observations on small astrophysical scales <cit.>.

Current experimental limits e.g. from searches for rare decays <cit.> put strong constraints on the mediator coupling to Standard Model (SM) particles. Nevertheless, the smallness of the mediator mass leads to a huge enhancement of direct detection cross sections, so that an observation of DM scattering may be possible in spite of the small couplings. In fact, if the mediator has sizeable couplings to DM, direct detection experiments can probe regions of parameter space inaccessible to other low-energy searches. As the recoil spectrum depends sensitively on both the DM and mediator mass, given an observation of DM scattering via a light mediator with sufficient statistics at a direct detection experiment, it may be possible to reconstruct both masses.

Cryogenic direct detection experiments, such as SuperCDMS <cit.>, CRESST <cit.> and EDELWEISS <cit.>, are particularly well-suited for this task due to their excellent energy resolution and low threshold <cit.>. In fact, a low energy threshold is even more important in the case of a light mediator than for a heavy one, because the recoil spectrum falls even more steeply and therefore the sensitivity can be considerably improved by lowering the threshold <cit.>. The excellent energy resolution, on the other hand, makes it possible to extract the maximum amount of information on the particle physics properties of DM from a successful discovery. In other words, cryogenic detectors are not only well-suited to explore models with light DM particles (see e.g. <cit.>), but also to probe light mediators.

The projected progress for the low-threshold technology implies that parameter points that are currently consistent with all experimental constraints may predict up to thousands of events in near-future detectors. In this paper we study the amount of information that can be extracted from such a signal, taking into account background uncertainties, astrophysical uncertainties and degeneracies with other particle physics parameters. We demonstrate that cryogenic experiments can probe the mediator mass precisely in the regions of parameter space relevant for DM self-interactions, potentially enabling us to infer the behaviour of DM on astrophysical scales with laboratory experiments.

Direct detection experiments in the context of self-interacting DM have been studied previously <cit.>, most notably in ref. <cit.>. Our work differs from these earlier studies in that we do not attempt to derive existing constraints but rather to explore the potential of future low-threshold detectors to infer the properties of self-interacting DM. For this purpose, we implement several present and future direct detection experiments in a realistic and efficient manner, in order to perform parameter reconstruction with a number of nuisance parameters. For similar studies in the context of effective operators see refs. <cit.>.

This paper is structured as follows. In section <ref> we discuss the phenomenology of direct detection experiments in the presence of light mediators. We review current and proposed low-threshold experiments and calculate their sensitivity to long-range interactions in comparison to conventional direct detection experiments. Section <ref> focusses on the potential of low-threshold experiments to determine the particle physics parameters of the DM particle and its interactions. We discuss the impact of experimental, theoretical and astrophysical uncertainties, introduce suitable nuisance parameters to represent them and assess their impact on our results. Finally, in section <ref> we connect our results to the idea of self-interacting DM. Additional details are provided in appendices <ref> and <ref>.



§ DIRECT DETECTION WITH LIGHT MEDIATORS



We consider a DM particle of mass m_DM scattering off nuclei via the exchange of a mediator with mass m_med. Throughout this paper we will focus on the case that the mediator has spin-independent couplings to both nucleons and DM. The differential event rate with respect to recoil energy E_R for DM scattering off a given target isotope T with mass m_T and mass fraction ξ_T is then given by

    dR_T/dE_R = ρ_0  ξ_T/2 π  m_DMg^2   F_T^2(E_R)/( 2   m_T   E_R + m_med^2 )^2 η ( v_min (E_R))

with ρ_0 = 0.3 GeV/cm^3 being the local DM density.

As long as the assumption of spin-independent interactions holds, the functional form of the differential event rate does not depend on the spin of the DM particle or the mediator nor on whether or not the DM particle is its own anti-particle. The numerical pre-factors, however, may differ for these different scenarios. We assume that these pre-factors have been absorbed into the definition of the effective low-energy coupling g, i.e. we take g to be defined via eq. (<ref>). The precise definition of g in terms of the fundamental parameters of the specific models that we discuss will be provided below.

In eq. (<ref>) the factor F_T^2(E_R) denotes the nuclear response function, which depends on the ratio of the mediator couplings to neutrons and protons, f_n / f_p. We parametrise this ratio via θ≡arctan f_n / f_p, such that models with f_n = 0 (such as a dark photon with kinetic mixing <cit.>) have θ = 0, whereas models with f_n = f_p (such as a light scalar mixing with the Higgs boson <cit.>) have θ = π / 4. In the limit of zero momentum transfer the nuclear response function is then given by F_T^2(0) = (Z_T cosθ + (A_T - Z_T) sinθ)^2, where Z_T and A_T denote the charge and mass number of the target isotope T. For non-zero momentum transfer, F_T^2(E_R) decreases due to a loss of coherence, which we parametrise by the standard Helm form factors <cit.>.

The final factor in eq. (<ref>) denotes the velocity integral, which is given by

    η(v_min) = ∫_v_min^∞d^3v f(𝐯)/v   ,

where v_min = √(m_T E_R/2  μ_T^2). Unless stated otherwise, we assume the velocity distribution f(𝐯) to be given by an isotropic Maxwell-Boltzmann distribution in the Galactic rest frame with v_0 = 220  km/s, cut off at the Galactic escape velocity v_esc = 544  km/s and transformed into the solar rest frame with v_obs =232  km/s.

For the purpose of this work, we will be most interested in the impact of the term ( 2   m_T   E_R + m_med^2 )^2 = (q^2 + m_med^2)^2 in the denominator, where q denotes the momentum transfer in the scattering process. This momentum transfer is bounded by , where  in terms of the low-energy threshold E_th and  with . If q^max is small compared to m_med, we recover the limit of contact interactions conventionally considered in the effective operator approach. Conversely, if m_med is negligible compared to q^min, scattering with large momentum transfer is suppressed by an additional factor of 1/q^4, leading to very steeply falling recoil spectra. In the intermediate regime we may hope to infer the properties of the mediator from the detailed shape of the recoil spectra. These interesting mediator masses are typically in the MeV range. For example, for m_DM = 4 GeV a germanium experiment (m_T ≈ 68 GeV) with E_th = 100 eV has q^min≈ 3.7 MeV and q^max≈ 20.7 MeV.

In real direct detection experiments, the detected recoil energy E_D may differ from the true physical recoil energy E_R due to the finite energy resolution of the detector. Assuming this resolution to be described by a Gaussian distribution with energy-dependent standard deviation σ(E_R), we can calculate the probability for a scattering event with energy E_R to be observed in the interval E_1 ≤ E_D≤ E_2:

    p(E_R, E_1, E_2) = 1/2[(E_2 - E_R/√(2)σ(E_R)) - (E_1 - E_R/√(2)σ(E_R))]   ,

where  denotes the error function. Given the total exposure of an experiment κ(E_R), which may again depend on the recoil energy, we can then calculate the total number of events expected in the interval [E_1, E_2]:

    N(E_1, E_2) = ∑_T ∫ p_T(E_R, E_1, E_2)  κ_T(E_R)  dR_T/dE_RdE_R  ,

where the sum is over all isotopes T in the target, weighted with appropriate factors ξ_T, and we allow both p and κ to depend on the isotope.

The two most sensitive cryogenic direct detection experiments are CRESST-II <cit.> and CDMSlite <cit.>. In the near future CRESST-III <cit.> and SuperCDMS SNOLAB <cit.> plan to significantly improve sensitivity.[A proposal for a similar effort with optimized EDELWEISS-III detectors was very recently put forward in ref. <cit.>. The suggested approach is very similar to the one employed by the SuperCDMS collaboration and we expect qualitatively similar results.] We describe our implementation of these experiments in appendix <ref>. For comparison, we also consider bounds from Xenon1T <cit.>, as well as projections for the future sensitivity of LZ <cit.>. We note that all projections have similar time scales and correspond to the sensitivity that may be achievable within the next five to ten years.



As a validation of our implementation, we show in figure <ref> the (projected) constraints on standard spin-independent interactions that we obtain from present and future direct detection experiments. For this purpose, we take m_med≫ q^max, set θ = π / 4 and g^2 = 2 π σ_p   m_med^4 / μ_p^2, where μ_p ≡ m_DM  m_p / (m_DM + m_p) and m_p denotes the proton mass. Our recalculation of existing constraints is in good agreement with the respective published exclusion limits. For SuperCDMS and CRESST-III we have chosen a rather conservative low-energy threshold of 100 eV in both cases, so that our projected exclusion limits are somewhat weaker than the official projections for small DM masses. Nevertheless, we find good agreement at larger DM masses, where the precise value of the threshold does not matter.



Having validated our implementation, we now show the corresponding constraints on the effective coupling g for smaller mediator masses (figure <ref>). As expected, we observe that with decreasing mediator mass direct detection experiments become sensitive to smaller values of the effective coupling g, with particularly large enhancement factors found for low-threshold experiments. Another important observation is that for m_med≪ q^min the shape of the recoil spectra (and hence the exclusion limit) becomes independent of m_med. We will show in more detail in section <ref> that the ability of direct detection experiments to reconstruct the mediator mass is limited to the range q^min≲ m_med≲ q^max.

We note that the constraints shown in figure <ref> depend on the assumed value of the coupling ratio θ. Here we have chosen θ = 0 corresponding for example to a vector mediator with kinetic mixing. In terms of the kinetic mixing parameter ϵ and the mediator-DM coupling g_DM the effective coupling g is then simply given by g = e  ϵ  g_DM, where e = √(4 πα) is the electromagnetic coupling. Different values of θ would typically enhance the sensitivity of heavy targets like tungsten relative to light targets like oxygen, except for specific values of θ that lead to destructive interference between proton and neutron contributions (see section <ref>).



§ RECONSTRUCTING PARTICLE PHYSICS PARAMETERS



From figure <ref> we make two central observations: first, if DM-nucleon scattering is due to the exchange of light mediators, cryogenic experiments will have the best sensitivity to such interactions up to DM masses of around 10 GeV. And second, compared to current bounds this sensitivity will improve by up to two orders of magnitude in terms of the effective coupling g, corresponding to up to four orders of magnitude in terms of the scattering rate. These observations immediately raise the question what we can hope to learn from a DM signal in low-threshold direct detection experiments. In this section we will answer this question by generating mock data and using this data to perform a parameter reconstruction.

To determine the compatibility of a specific particle physics hypothesis (characterized by a set of parameters 𝐱) with a given set of data, we construct a likelihood function ℒ(𝐱) as follows. For each individual experiment α, we calculate the Poisson likelihood

    -2 logℒ^α(𝐱, 𝐲) = 2 ∑_i [ R^α_i(𝐱, 𝐲) + B^α_i(𝐲) - N^α_i + N^α_i logN^α_i/R^α_i(𝐱, 𝐲) + B^α_i(𝐲)]   ,

where the sum is over all bins, and R^α_i, B^α_i, and N^α_i denote the number of predicted signal events, predicted background events and observed events, respectively. In addition to the particle physics parameters 𝐱 we have introduced a number of nuisance parameters 𝐲, which represent for example astrophysical or experimental uncertainties and may affect both signal and background predictions. These nuisance parameters may be constrained by an additional likelihood function ℒ^n. The total profile likelihood is then given by the product of the individual likelihoods maximised with respect to the nuisance parameters:

    ℒ(𝐱) = max_𝐲ℒ^n(𝐲) ∏_αℒ^α(𝐱, 𝐲)   .


For the purpose of parameter estimation, the next step is to determine the value of the particle physics parameters 𝐱 that maximise the profile likelihood, called 𝐱_0. We can then construct the likelihood ratio ℛ(𝐱) = ℒ(𝐱) / ℒ(𝐱_0), which by definition is smaller than unity. Under random fluctuations in the data, the quantity -2 logℛ(𝐱) is expected to follow a χ^2 distribution with number of degrees of freedom n equal to the number of parameters 𝐱. We can therefore exclude a hypothetical set of parameters 𝐱 at confidence level 1-p if

    1 - CDF_χ^2(n, -2 logℛ(𝐱)) < p   ,

where CDF_χ^2(n, x) denotes the cumulative distribution function for the χ^2 distribution with n degrees of freedom. For the case of two parameters, the 95% confidence level (CL) bound is therefore given by -2 logℛ < 5.99.

In the following we will focus on m_DM≲ 5 GeV, where cryogenic detectors have better sensitivity than liquid xenon experiments (see figure <ref>). We will first focus on one specific benchmark case, namely m_DM = 2 GeV, m_med = 3 MeV and θ = 0, and then discuss alternative benchmarks in section <ref>. The assumed value of g is chosen to be compatible with existing direct detection constraints. Choosing g close to current exclusion limits will lead to an optimistic scenario, in which thousands of events can be observed in future experiments, whereas smaller values of g imply smaller statistics and less precise parameter reconstruction. In the following, we will consider two alternative choices, namely g = 2 · 10^-11 (referred to as the low-statistics case) and g = 6 · 10^-11 (the high-statistics case). For our benchmark scenario, these choices correspond to around 900 and 8000 events across the set of future experiments that we consider (with SuperCDMS SNOLAB predicted to observe about four times as many events as CRESST-III).

We can now generate mock data for our benchmark scenario and the two possible coupling choices and then determine which alternative choices of m_DM and m_med are compatible with this data. For the purpose of parameter reconstruction it is sufficient to consider mock data sets without Poisson fluctuations. Although in this case the best-fit point will have a very high likelihood, ℒ≈ 1, we nevertheless obtain reasonable estimates of the likelihood ratio ℛ(𝐱) expected in a typical realization of the experiments. We will return to the issue of Poisson fluctuations in the context of goodness-of-fit estimates in section <ref>.



Figure <ref> shows the regions of parameter space compatible with the mock data generated for our benchmark scenario. For the purpose of these plots we are not interested in reconstructing the effective coupling g, i.e. we will simply treat it as a nuisance parameter and fix it to the value that maximises the likelihood. Nevertheless, the assumed value of g does play an important role as it determines the number of events that we expect to observe. The left (right) panel corresponds to the low-statistics (high-statistics) case. Red and blue contours correspond to the parameter reconstruction based only on data from SuperCDMS SNOLAB and CRESST-III respectively, while the grey region indicates the combined constraints. Note that in these plots we do not yet take into account nuisance parameters related to background or astrophysical uncertainties; these will be discussed later in this section.

A striking feature in the left panel of figure <ref> is the accuracy of the parameter reconstruction from SuperCDMS SNOLAB as compared to CRESST-III. This happens because of two reasons: first, SuperCDMS SNOLAB is predicted to observe about four times more events than CRESST-III and hence has much better statistics. Second, several target elements contribute to the observed event rates in CRESST-III, leading to different ways in which a good fit to the observed data can be obtained. While for the benchmark case that we consider the event rate is dominated by scattering off oxygen (because tungsten recoils are below threshold), very similar recoil spectra are obtained for heavier masses and scattering off tungsten. This observation also explains the two different `branches' found for heavy mediator masses. With sufficient statistics it becomes possible to distinguish between the two possible scenarios and reject the solution with dominant scattering off tungsten (see right panel on figure <ref>).

Another interesting feature is that all reconstructed parameter regions exhibit a characteristic `tilt' in the sense that lighter mediators are necessary for fitting heavier DM masses and vice versa. The origin of this shape is that heavier DM masses predict flatter recoil spectra, while lighter mediators predict steeper recoil spectra. Increasing the DM mass while decreasing the mediator mass and the effective coupling g may therefore leave the recoil spectra approximately unchanged. This degeneracy disappears once the mediator becomes so light (or so heavy) that direct detection experiments are effectively in the massless mediator limit (or the contact interaction limit). The recoil spectra then no longer depend on the precise value of the mediator mass.

Finally we make the crucial observation that combining data from SuperCDMS SNOLAB and CRESST-III allows for a much more precise reconstruction of the mediator mass than considering the individual experiments. The primary reason behind this is that the degeneracy between DM mass, mediator mass and effective coupling strength g discussed above depends on the target element and therefore on the experiments (see eq. (<ref>)). This degeneracy is therefore effectively removed when combining data from several experiments. Nevertheless, it is of course conceivable that the degeneracy will reappear (or new degeneracies will arise) once we include various sources of uncertainties. We will therefore now discuss such uncertainties in detail and assess their impact on our results.



 §.§ Background uncertainties



In the parameter reconstruction performed above we have assumed exact knowledge of the shape and normalization of the experimental background(s). This is an overly optimistic assumption and we will now discuss a more conservative approach. Of course, if we were to allow arbitrary background shapes and normalizations, any kind of DM signal could be absorbed into the background, making it impossible to claim anything but exclusion limits. Any parameter reconstruction therefore necessarily requires some knowledge on the distribution of backgrounds. Here we assume that the shape of the background in each experiment is known, but we keep the normalization completely free. In other words, we introduce a nuisance parameter y^α for each experiment α such that the background predictions in eq. (<ref>) are given by B^α_i(𝐲) = y^α B^α_i. We do not impose any constraints on the parameters y^α apart from the trivial requirement that they must be strictly positive. As explained in more detail in appendix <ref>, we assume the backgrounds to be flat in energy both for CRESST-III and SuperCDMS SNOLAB. However, this assumption can easily be modified within our framework once more detailed informations about the future experiments are available.



We show the impact of including background uncertainties in figure <ref>. As expected, these additional nuisance parameters visibly increase the size of the allowed parameter regions, in particular for SuperCDMS SNOLAB, where only a relatively small range of recoil energies is used to constrain the DM properties. For CRESST-III we observe that the second branch with scattering dominantly on tungsten now appears also in the high-statistics case. In principle it may be possible to distinguish these two branches, because the light yield of tungsten recoils in CRESST-III differs from the one for recoils on oxygen and calcium. Including this information (for example by constructing a two-dimensional likelihood in recoil energy and light yield) may hence make it possible to remove the second branch and obtain only one best-fit region. We leave this interesting possibility for future work.

In principle, given additional information on the different background contributions, it would be straight-forward to extend our approach and introduce individual nuisance parameters for each background source, such that not only the normalization of the background but also its shape can be varied. However, such detailed information is not presently available for the future projections that we consider. Nevertheless, we note in passing that the likelihood function defined in eq. (<ref>) has an interesting property: as long as we only introduce nuisance parameters that rescale the signal prediction or (parts of) the background prediction, -2 logℒ is a convex function of these parameters, so that any local minimum is necessarily a global one <cit.>. It is therefore numerically trivial to maximize the likelihood with respect to these nuisance parameters.



 §.§ Degeneracies with coupling ratios



So far we have made the assumption that only the mass of the mediator is unknown but its coupling structure is fixed. However, even in the simplest models of light mediators there are a number of different possible coupling structures <cit.>. In fact, it is well-motivated to assume that the ratio of couplings to neutrons and protons is essentially a free parameter <cit.>. If the effect of varying the mediator mass can be compensated by changing the coupling ratio, our lack of knowledge concerning the coupling structure of the mediator may affect our ability to determine its mass. In this section we discuss how the parameter reconstruction is affected if we do not make any assumptions on θ.

For experiments consisting only of a single target element with charge number Z_T and mass number A_T, the differential event rate is approximately proportional to[Strictly speaking, this relation only holds exactly for zero momentum transfer, because at finite momentum transfer there may be differences in the form factors for protons and neutrons <cit.>, but taking these effects into account is beyond the scope of this work.] (Z_T cosθ + (A_T - Z_T) sinθ)^2. This rescaling factor takes different values for different elements, so that varying θ will affect the combination of information from different experiments. In fact, varying θ can even change the shape of the recoil spectrum in a single experiment that consists of several different elements.

In contrast to the nuisance parameters introduced to parametrize the background uncertainties, the parameter θ enters into the likelihood in a more complicated way. The reason is that for θ < 0 there will be destructive interference between proton and neutron contributions, which for a given experiment will be maximal if θ = arctan(- Z_T / (A_T - Z_T)). Consequently, there will generally be a number of different values of θ that maximize ℒ(θ) locally. We therefore perform an explicit scan over the full range -π/2 < θ < π / 2 to determine the global maximum of ℒ(θ).



If we consider SuperCDMS SNOLAB alone, the likelihood is essentially independent of θ, because the recoil spectra for the different isotopes look so similar that a change in θ can almost entirely be compensated by a change in the effective coupling g. Indeed, we will see below that the allowed parameter regions for SuperCDMS SNOLAB remain essentially unchanged when including θ as a nuisance parameter. For CRESST-III, on the other hand, the situation is very different. Since the recoil spectrum for scattering of tungsten is typically much steeper than the recoil spectrum for scattering of oxygen and calcium, the shape of the recoil spectrum depends sensitively on θ. Figure <ref> shows the best-fit value of θ depending on the assumed values of m_DM and m_med. We can make a number of interesting observations:

 
  * The best-fit value for θ typically differs significantly from the value assumed to generate the mock data, covering almost the entire range - π/2 ≤θ≤π/2.
 
  * If the assumed DM mass is large compared to the DM mass used to generate the mock data, the predicted recoil spectrum will be too flat, so the best fit can be obtained if scattering occurs exclusively on tungsten, i.e. for θ = -π/4.
 
  * If the assumed mediator mass is small compared to the mediator mass used to generate the mock data, the predicted recoil spectrum will be too steep, so the best fit can be obtained if scattering occurs exclusively on oxygen and calcium, i.e. for θ∼ -0.59.


It should be clear from these observations that the ability of CRESST-III to reconstruct the mediator mass and the DM mass is significantly reduced when allowing for arbitrary values of θ. This is confirmed explicitly in figure <ref>. Clearly, in this case it is essential to combine the information from CRESST-III with data from SuperCDMS SNOLAB to perform a precise parameter reconstruction. In the absence of data from SuperCDMS SNOLAB it could also be interesting to attempt a discrimination between the different target elements in CRESST-III or to perform a combination of CRESST-III with the existing exclusion limits from Xenon1T, which would disfavour solutions with heavy DM scattering on tungsten.





 §.§ Astrophysical uncertainties



There are two kinds of astrophysical uncertainties that may affect the interpretation of direct detection experiments: uncertainties in the local DM density ρ_0 and uncertainties in the DM velocity distribution f(𝐯). The local DM density ρ_0 enters linearly into the differential event rate, so varying this quantity is equivalent to varying the effective coupling g. Since we are not interested in determining g and simply treat it as a nuisance parameter, our approach therefore already accounts also for uncertainties in the local DM density. The velocity distribution, on the other hand, enters in a more complicated way, giving rise to additional uncertainties that we will now discuss.

The differential event rate depends on the velocity distribution via the velocity integral η(v_min), which in turn depends on the DM mass and the recoil energy via v_min = √(m_T   E_R/2 μ_T^2). Changes in the velocity distribution may therefore change the shape of the recoil spectrum and thereby limit our ability to extract information on the particle physics parameters. One possible way to deal with astrophysical uncertainties are so-called halo-independent methods <cit.>, which aim to combine information from different experiments in such a way that the dependence on η(v_min) drops out. This approach has been successfully applied to many different models and in particular to models with light mediators <cit.>.

However, as pointed out in ref. <cit.>, there is a fundamental limitation of this approach in the case of low-mass DM. If m_DM≪ m_T for all target elements under consideration, we find that v_min≃√(m_T   E_R/2   m_DM^2) = q/2   m_DM. The velocity integral hence depends on the same combination of m_T and E_R that also enters in the factor for light mediator exchange, eq. (<ref>). In other words, for low-mass DM, any change in the mediator mass can be compensated for by an appropriate change in the velocity integral for all target elements simultaneously. As a result, it will not be possible to infer any information on the mediator mass if we allow for completely arbitrary velocity distributions.

We will therefore take a different approach and consider only velocity distributions of the Maxwell-Boltzmann form. This assumption is supported by recent studies involving numerical simulations of Milky Way like galaxies <cit.>.  We can then study the impact of astrophysical parameters by varying the underlying parameters v_0, v_esc and v_obs. In this case we can actually use the fact that m_DM≪ m_T to our advantage. As shown in appendix <ref>, if we simultaneously rescale all three velocities by a common factor z, this change is fully equivalent to rescaling the DM mass by a factor z. We therefore introduce a new nuisance parameter z and, rather than calculating the differential event rate as a function of m_DM, we calculate the differential event rate as a function of z   m_DM. We restrict z to lie in the range consistent with observations. At 95% CL v_0 is constrained to lie in the range [180 km   s^-1, 280 km   s^-1] while the range for v_esc is approximately [450 km   s^-1, 650 km   s^-1], see <cit.> and references therein. These ranges can be reproduced if we require 0.8 ≤ z ≤ 1.2 at 95% CL. We implement this by means of a likelihood function for z given by

    ℒ^z = 1/√(2π)σ_zexp(- (z - 1)^2/2  σ_z^2)

with σ_z = 0.1 and include this extra factor in the total likelihood.



As in the case of θ it is possible that the likelihood has several local maxima for different values of z, making it necessary to explicitly scan over all values of z within the relevant range. Nevertheless, the simple way in which ℒ depends on z means that it is not in fact necessary to perform a two-dimensional scan over both θ and z, but rather that two separate one-dimensional scans are sufficient. Figure <ref> shows the impact of including astrophysical uncertainties in addition to the uncertainties discussed above. As expected, the effect of varying z is essentially to reduce our ability to reconstruct the DM mass, while not affecting our ability to reconstruct the mediator mass. Figure <ref> constitutes our central result for the benchmark scenario: even when including a number of different nuisance parameters, an accurate reconstruction of the DM and mediator masses is possible given sufficient statistics.



 §.§ Alternative benchmark scenarios



In the discussion above we have introduced a number of nuisance parameters that should be taken into account for a realistic assessment of the power of future low-threshold direct detection experiments. In addition to the two parameters that we are interested in reconstructing (the DM mass and the mediator mass), we have introduced two particle physics nuisance parameters (the coupling strength g and the coupling ratio parameter θ), one astrophysics nuisance parameter (the rescaling factor z) and one experimental nuisance parameter for each experiment (the background normalizations y^α).[Even in the presence of these nuisance parameters it only takes a few seconds on a single CPU to calculate the profile likelihood for a single grid point and around an hour to perform the parameter reconstruction on a grid with 10^4 points. A significant amount of computing time can be saved by reusing the information from neighbouring grid points to profile out the nuisance parameter related to astrophysical uncertainties.]  In this section we present our results for a number of different hypotheses on the particle physics properties of DM and discuss the physics interpretation.



In figure <ref> we perform a parameter reconstruction of the DM and mediator masses for two different assumptions on the true DM mass, namely 1 GeV (left panel) and 3  GeV (right panel). In both cases we fix the mediator mass to m_med = 3MeV and the effective coupling to g = 6 · 10^-11.[It should be noted that this procedure leads to somewhat different numbers of events in the left and right panel. In particular the event rate in SuperCDMS is significantly suppressed for m_DM = 1 GeV.] For m_DM = 1 GeV all observed events are very close to the low-energy threshold (i.e. within the first two or three bins). As a result the parameter reconstruction becomes more difficult and neither of the two experiments can individually constrain the mass of the mediator. For CRESST-III one furthermore finds a second branch of solutions corresponding to scattering off tungsten. Combining the information from both experiments leads to a somewhat better reconstruction, but the allowed parameter region still extends to arbitrarily heavy mediators. For heavier DM masses, on the other hand, an accurate parameter reconstruction is possible (see the right panel of figure <ref>).



In figure <ref> we investigate the effect of varying the assumed mediator mass while keeping the DM mass fixed to 2 GeV. For each mediator mass we choose the value of the coupling g such that the predicted number of events is comparable to the high-statistics case discussed previously. In the top-left panel, the mediator mass is set to 0.1 MeV, i.e. effectively massless for the experiments under consideration. The combined fit to both experiments then places an upper bound on the mediator mass of about m_med < 0.8  MeV. As expected, the mock data is compatible with arbitrarily light mediators so that no lower bound can be placed. Conversely, if the assumed mediator mass is larger than about 10 MeV, it is no longer possible to distinguish our scenario from the case of contact interactions and the allowed parameter region extends up to arbitrarily high mediator masses (see bottom row of figure <ref>). An accurate reconstruction of the mediator mass is possible only if the mass falls between these two extremes, as illustrated in the top-right panel for m_med = 1 MeV.

To conclude this section we note that we have also studied the effect of making different assumptions on the value of θ used to generate mock data. If the effective coupling g is adjusted in such a way that the event numbers are comparable to the ones discussed above, we find very similar results for different choices of θ. An example with θ≠ 0 will be discussed in section <ref>.



 §.§ Goodness-of-fit estimates



So far we have focussed on the issue of parameter estimation, i.e. we have constructed likelihood ratios to determine the parameter regions compatible with a given set of mock data. Another interesting topic that we can study in our framework are goodness-of-fit estimates, i.e. the question whether a specific choice of parameters yields a good description of the data. To answer this question we need to consider the absolute value of the likelihood rather than a likelihood ratio. Clearly, it is then no longer possible to neglect Poisson fluctuations in the data, because doing so would exaggerate the likelihood, i.e. would suggest unrealistically good agreement between data and model. In this section we therefore briefly discuss the effect of Poisson fluctuations and give a few examples for questions that can be answered with goodness-of-fit estimates.

In the limit of large bin counts the likelihood function defined in eq. (<ref>) approaches a χ^2 test statistic. We therefore expect that for the true parameters of nature 𝐱_0 the likelihood ℒ(𝐱_0) follows a χ^2-distribution with n = n_b - n_y degrees of freedom, where n_b denotes the total number of bins across all experiments and n_y denotes the total number of nuisance parameters that have been profiled out. For the two experiments that we consider n_b = 29, so if we profile out the effective coupling strength g and the normalization of the background in both experiments, we expect to find n = 26. We confirm that this is indeed the case by performing a Monte Carlo (MC) simulation, i.e. by considering a large ensemble of mock data sets with Poisson fluctuations.

We can make use of this observation to study whether a given data set may enable us to exclude specific hypotheses about the particle physics nature of DM. For example, any DM signal observed in future direct detection experiments will first be interpreted under the assumption of a heavy mediator, i.e. contact interactions between DM and nuclei. A question of great interest would therefore be whether this hypothesis can be confidently excluded if the mediator is in fact light. To answer this question, we can use an ensemble of mock data sets generated under the assumption of a light mediator and calculate the likelihood under the incorrect assumption of a heavy mediator. For this purpose, we treat both m_DM and θ as nuisance parameters, i.e. we fix them to the value that maximizes the likelihood, such that n = 24.[Note that due to the degeneracy between the DM mass and the nuisance parameter z parametrizing astrophysical uncertainties, including z as an additional nuisance parameter would not increase the likelihood further and therefore does not reduce the relevant number of degrees of freedom.] For any specific mock data set, the heavy-mediator assumption can then be excluded at 95% CL if -2 logℒ > 36.4.

Considering the same benchmark scenario as before (m_DM = 2 GeV, m_med = 3 MeV and θ = 0) we find that in the low-statistics case it is typically not possible to exclude the heavy-mediator hypothesis at 95% CL, whereas in the high-statistics case a 95%-CL exclusion is possible for about 98% of the mock data sets. In the latter case it is also possible to exclude the hypothesis of a very light mediator (with m_med≪ q^min) at 95% CL for more than 99.9% of the samples.



§ CONNECTION TO SELF-INTERACTING DARK MATTER



The analysis performed in section <ref> applies to any model of DM interacting with nucleons via a light mediator, provided that the couplings of the mediator are spin-independent. In this section we apply these results to a particle physics model of particular interest, namely the case of self-interacting DM. 

It is well known that the presence of a light mediator can significantly enhance the rate of DM self-scattering in astrophysical systems, not only because of the long-range nature of the interactions but also because of additional non-perturbative effects due to multiple mediator exchange <cit.>. For a mediator with spin-independent interactions, these effects can be calculated by solving the non-relativistic Schroedinger equation for a Yukawa potential. The resulting scattering rate exhibits a characteristic velocity dependence such that the largest effects are expected on small scales, just as required in order to resolve the potential small-scale problems of collisionless cold DM <cit.>. The quantity of interest for astrophysical observables is the momentum transfer cross section σ_T, which should lie in the range 0.1  cm^2/g≲σ_T / m_DM≲ 10  cm^2/g in order to induce sizeable effects consistent with observations.

Out of the various possible particle physics realizations of this general scenario, the most attractive model consists of a fermionic DM particle ψ (here assumed to be a Dirac fermion) and a scalar mediator ϕ <cit.>. For this model, DM self-annihilation is suppressed at small velocities due to CP conservation, evading the strong constraints from the Cosmic Microwave Background on many DM models with light mediators <cit.>. The effective coupling g is then given by g ≈ 1.6 × 10^-3  y_ψ  y_SM, where y_ψ denotes the DM-mediator coupling and y_SM denotes the rescaled Yukawa coupling of the mediator to Standard Model (SM) fermions (i.e. the mixing of the scalar mediator with the SM Higgs boson).

A recent analysis of this model <cit.> concluded that for DM masses in the GeV region and mediator masses in the MeV region it is possible to reproduce the observed DM relic abundance and at the same time obtain large DM self-interactions if y_ψ∼ 0.1. Bounds from direct detection experiments then require y_SM≲ 10^-6 so that constraints from Higgs measurements and flavour physics are easily satisfied.[We note however that for such small values of y_SM the lifetime of the mediator will be so large that it only decays after the onset of primordial nucleosynthesis. Nevertheless, the mass of the mediator is so small that it can only decay electromagnetically and with an energy small compared to typical nuclear binding energies. The dominant constraint is therefore expected to come from the entropy injected by the mediator decays, which modify the effective relativistic degrees of freedom N_eff. The magnitude of this effect depends however on additional aspects of the model, which are not relevant for the direct detection phenomenology. A detailed analysis of constraints from primordial nucleosynthesis is challenging and will be discussed elsewhere <cit.>.] This mass range is precisely the parameter region that we have identified above to be the most interesting for cryogenic direct detection experiments. In fact, these experiments already place the strongest bounds on the SM-mediator coupling in this scenario. For the purpose of this section, we will therefore focus on scalar mediators and set θ = π/4 as expected from Higgs mixing <cit.>, both for the generation of mock data and for the subsequent parameter reconstruction.



While the DM momentum transfer cross section depends only on the DM-mediator coupling y_ψ, the DM-nucleon scattering cross section depends also on the mediator-SM coupling y_SM, which arises for example from mixing of ϕ with the SM Higgs boson. The effective coupling g relevant for direct detection experiments can therefore be treated as an independent parameter in the same way as we have done above. Thus, if we can use direct detection experiments to infer the mediator mass and the DM mass, we can effectively determine the DM momentum transfer cross section (under the assumption that DM is a thermal relic <cit.>). This is illustrated in figure <ref>. In the left panel we show a plot very similar to the ones shown in section <ref>, except that we also indicate the combinations of DM mass and mediator mass that lead to a momentum transfer cross section in the desired range <cit.>. In the right panel, we show the range of momentum transfer cross sections σ_T compatible with the inferred range of masses as a function of the DM relative velocity v. 

Since the DM self-interaction cross section depends very sensitively on the mediator mass, a precise determination of σ_T is clearly impossible. Nevertheless, figure <ref> shows that it may be possible to demonstrate that both the magnitude and the velocity dependence of the DM self-scattering rates are broadly in agreement with astrophysical requirements, at least for the specific assumptions that we have made. Direct detection experiments may hence help us to determine not only the properties of the DM particle itself, but also to learn about its role in the formation of small-scale structures.



§ CONCLUSIONS



Light mediators communicating the interactions between DM and nuclei offer an interesting alternative to the effective operator approach conventionally adopted for the analysis of direct detection experiments. In this work we have demonstrated that cryogenic direct detection experiments are particularly well-suited for exploring this scenario. Their low threshold enhances the sensitivity to steeply falling recoil spectra and their excellent energy resolution allows for a precise reconstruction of the underlying particle physics from a potential signal. Present searches based on this technology already place stringent constraints on models with light mediators and significant improvements are expected for the next few years.

To illustrate the potential of future low-threshold detectors, we have performed parameter reconstructions from mock data sets in two planned experiments: SuperCDMS SNOLAB and CRESST-III. A strong emphasis has been placed on a realistic implementation of the experimental details. We include detector effects such as energy resolution, nuisance parameters for the background contributions as well as additional uncertainties reflecting both the unknown particle physics properties of DM and its astrophysical distribution. Even when including all of these uncertainties, low-threshold direct detection experiments maintain a remarkable potential to reconstruct the properties of the mediator.

For mediator masses in the range 1–10 MeV, the possibility to reconstruct mediator masses is mostly limited by statistics. Given sufficient exposure, it should be possible to reconstruct the mediator mass to better than a factor of 2. If the mediator is either lighter than about 1 MeV (and hence indistinguishable from truly long-range interactions) or heavier than about 10 MeV (and hence indistinguishable from contact interactions), a precise reconstruction of its mass is typically not possible. Nevertheless, with enough statistics one can clearly rule out the hypothesis of contact interactions for a very light mediator (and vice versa) using goodness-of-fit estimates.

The ability to reconstruct mediator masses relies strongly on the combination of different target materials. We have focussed on the combination of experiments based on CaWO_4 and germanium targets, but we have checked that a similarly good reconstruction can be obtained by combining data from a germanium experiment with results based on silicon detectors, which are being developed by the SuperCDMS collaboration. Another interesting, albeit challenging, possibility would be to exploit the differences in light yield for the different target elements in CRESST-III to statistically distinguish between models with dominant scattering on tungsten and those with dominant scattering on oxygen and calcium.

Throughout this work we have focussed on mediators with spin-independent interactions. Spin-dependent interactions are more challenging, since neither oxygen nor calcium nuclei carry spin and therefore the sensitivity of CRESST-III is significantly suppressed. It may be possible to gain complementary information with bubble chamber experiments like PICO-500 <cit.>, but these lack the excellent energy resolution of cryogenic detectors. The most promising avenue appears to be the combination of germanium and silicon detectors, both of which possess at least some sensitivity to spin-dependent interactions. A detailed investigation of a wider range of interactions and a broader set of experimental technologies is left for future work.

Finally, we have pointed out that the range of mediator and DM masses that can be studied with cryogenic experiments coincides with parameter regions that have been considered in the context of self-interacting DM. This observation implies that, within specific model assumptions, one can translate between signals in direct detection experiments and astrophysical observables, such as core sizes in dwarf galaxies. By correlating these very different signatures, we can therefore hope to ultimately obtain a coherent picture of both the microscopic and macroscopic properties of DM.

We thank Achim Gütlein for early collaboration and David Cerdeño, Josef Jochum and Kai Schmidt-Hoberg for useful discussions. FK is supported by the
DFG Emmy Noether Grant No. KA 4662/1-1, SK is supported by the `New Frontiers' program of the Austrian Academy of Sciences and by FWF project number V592-N27, and SW is supported by the ERC Starting Grant `NewAve' (638528). SK thanks DESY for hospitality during the completion of this work.



§ IMPLEMENTATION OF PRESENT AND FUTURE EXPERIMENTS



In this appendix we describe the experimental information and the assumptions that we make to construct the likelihood functions for the various experiments. 



 §.§ CRESST-II


The latest CRESST results <cit.> are based on 52.2 kg days using the Lise detector module. We smear the physical recoil spectrum assuming a Gaussian energy resolution function with σ = 62eV <cit.>, considering only fluctuations up to 3σ in order to avoid unphysical results far below threshold. The cut survival probabilities as a function of the detected recoil energy are taken from the recent data release <cit.>. We then group the events observed in the acceptance region into 47 equidistant bins between 0.3 and 5.0 keV. For the calculation of the upper limits we conservatively set the number of expected background events to zero, which gives results similar to the optimum interval method <cit.> employed in ref. <cit.>.



 §.§ CRESST-III


Our implementation of the future CRESST detector is based on the projections for the final state of CRESST-III presented in ref. <cit.>. The expected exposure is 1000 kg days with an energy threshold of 100eV. We group the events in the energy range between 0.1 and 2keV into 19 bins of width 0.1 keV. Furthermore, we assume a background level of , corresponding to 3.5 events in each of the bins. The expected number of signal events is calculated assuming an efficiency of unity down to threshold, and using a Gaussian energy resolution of 20eV. We take into account only recoils with a true energy above 60eV in order to avoid unphysical upward fluctuations from very small recoil energies.



 §.§ CDMSLite


CDMSLite has analysed data from an exposure of 70.14 kg days <cit.>. We take the energy-dependent signal efficiency from <cit.> and follow the procedure outlined there to convert nuclear recoil energies (eVnr) into electron equivalent energies (eVee), using k = 0.157. We perform a fit to the width of various electron-capture peaks to determine the energy-dependent detector resolution. We then consider the energy range [60 eVee, 500 eVee], which we divide into 10 bins of increasing size, such that the number of observed events in each bin is greater than one. Following the analysis in ref. <cit.>, we do not assume a background model, such that all observed events can potentially be DM signals. This procedure allows us to calculate a likelihood function that yields a bound very similar to the one obtained from the optimum interval method. 



 §.§ SuperCDMS SNOLAB


We estimate the sensitivity of SuperCDMS SNOLAB following ref. <cit.>, focussing on the high-voltage germanium detectors. Although the SuperCDMS collaboration expects that a threshold as low as 40  eV (nuclear recoil energy) can be achieved, we consider a low-energy threshold of 100 eV, so that it is a good approximation to assume the signal efficiency to be constant (at 85%) and the background level to be flat (at ) <cit.>. We also limit ourselves to recoil energies below 300 eV in order to avoid backgrounds from electron capture lines. This restriction will reduce the sensitivity of SuperCDMS SNOLAB to heavy DM but has no effect in the low-mass region that we are interested in. We divide the search window into 10 equally-spaced bins and include an energy resolution of 10 keV. SuperCDMS expects to achieve a total exposure of 1.6 · 10^4  kg days over the course of five years of data taking. To compare CRESST-III and SuperCDMS on similar time-scales, however, we consider the sensitivity that can be achieved with a single year of data taking, corresponding to an exposure of approximately 3200  kg   days.[Incidentally, this is also the exposure for which backgrounds become non-negligible (we expect about 2 events per bin), so increasing the exposure by a factor of 5 will improve the sensitivity by significantly less than a factor of 5.]



 §.§ Xenon1T


The first results of Xenon1T are based upon an exposure of 35636 kg days <cit.>. To calculate the acceptance function, we simulate fluctuations in the S1 and S2 signal making use of the scintillation and ionization yields from ref. <cit.> and taking into account the anti-correlation between the two signals. Rather than attempting to model the detector response, we determine the light collection efficiency and the electron extraction probability by fitting to the nuclear recoil band shown in ref. <cit.>. We then consider events with S1 ≥ 3 phe and require the S2 signal to lie below the mean of the nuclear recoil band. For low recoil energies, this procedure essentially reproduces the acceptance function shown in ref. <cit.>, where no cut on the S2 signal is applied, whereas for large recoil energies we obtain an acceptance that is approximately a factor of 2 smaller. This result is in agreement with the expectation that the nuclear recoil cut has almost 100% signal acceptance close to the threshold, where the signal results exclusively from upward fluctuations in the S1 signal, and around 50% signal acceptance away from the threshold. The total number of expected background events below the mean of the nuclear recoil band is 0.36, and no events have been observed.
We therefore set a Poisson upper bound on the signal contribution of 1.94 events at 90% CL, noting that this procedure leads to a slightly weaker bound than the asymptotic limit assumed in ref. <cit.>.



 §.§ LZ


Our implementation of LZ is based upon the LZ conceptual design report <cit.>. Specifically, we assume a total exposure of 5.6 · 10^6  kg days and consider the search region 3 phe≤ S1 ≤ 30 phe. We calculate the energy-dependent acceptance by considering Poisson fluctuations in the S1 signal based on the effective scintillation yield for liquid xenon taken from the LUX analysis <cit.> and assuming a light collection efficiency of 7.5%. We multiply this acceptance function by an overall factor of 1/2 to account for the requirement that the S2 signal must lie below the mean of the nuclear recoil band. Finally, we calculate the expected sensitivity by assuming that the number of observed events is equal to the expected background (2.37 events).



§ ASTROPHYSICAL UNCERTAINTIES FOR LOW-MASS DARK MATTER



In the Galactic rest frame the standard Maxwell-Boltzmann velocity distribution is given by

    f(𝐯, v_0, v_esc) = 1/N[ exp(-v^2/v_0^2) - exp(-v^2/v_esc^2) ]   ,

where N is an appropriate normalization constant. In the laboratory frame we need to account for the velocity of the Earth 𝐯_E relative to the Galactic rest frame. Neglecting the motion of the Earth relative to the Sun, we take this velocity to be time-independent, 𝐯_E = 𝐯_obs. We then obtain

    f(v, θ_v, v_0, v_esc, v_obs) = 1/N[    exp(-v^2 + 2   v   v_obscosθ_v + v_obs^2/v_0^2) 
        - exp(-v^2 + 2   v   v_obscosθ_v + v_obs^2/v_esc^2) ]   ,

where θ_v denotes the angle between 𝐯_obs and the velocity vector 𝐯.

Let us now consider two different velocity distributions of the Maxwell-Boltzmann form, f(v) and f̃(v), which are related by

    ṽ_0 = v_0/z   ,     ṽ_esc = v_esc/z     ,     ṽ_obs = v_obs/z   ,

where z is an arbitrary rescaling factor.[We note that for an isothermal halo the velocity dispersion and the circular velocity are directly proportional to each other and it is therefore well-motivated to rescale v_0 and v_obs simultaneously. The escape velocity is introduced by hand and is therefore in principle an independent parameter. Since our results do not depend sensitively on v_esc we use the same rescaling factor for simplicity.] It then follows immediately from eq. (<ref>) that the two velocity distributions satisfy the relation

    f̃(v) = z^3 f(z v)   ,

where the pre-factor ensures that f̃(v) is normalized, ∫ f(v)  d^3 v  = ∫f̃(v)  d^3 v = 1. The velocity integral η̃(v_min) corresponding to f̃(v) is then given by

    η̃(v_min)     = ∫_v_min^∞f̃(v)/vd^3 v = ∫_v_min^∞ z^3 f(z v)/vd^3 v = ∫_ṽ_min^∞αf(ṽ)/ṽd^3 ṽ = z η(ṽ_min)  ,

where ṽ_min = z v_min.

The differential event rate resulting from the rescaled velocity distribution therefore becomes

    dR̃_T/dE_R = ρ_0  ξ_T/2 π  m_DMg^2   F_T^2(E_R)/( 2   m_T   E_R + m_med^2 )^2  z η (ṽ_min (E_R, m_DM))

The crucial observation is now that for m_DM≪ m_T the minimum velocity is given by v_min(E_R, m_DM) = 1/m_DM√(m_T   E_R/2)
and hence ṽ_min(E_R, m_DM) = v_min(E_R, m̃_DM) with m̃_DM = m_DM / z. In terms of this rescaled mass, the differential event rate becomes

    dR̃_T/dE_R = ρ_0  ξ_T/2 π m̃_DMg^2   F_T^2(E_R)/( 2   m_T   E_R + m_med^2 )^2 η (v_min(E_R, m̃_DM))   ,

which is identical to eq. (<ref>) with m_DM replaced by m̃_DM. In other words, as long as m_DM≪ m_T, the effect of an overall rescaling in the velocity distribution is equivalent to a rescaling in the DM mass.

10


Fan:2010gt
J. Fan, M. Reece, and L.-T. Wang,
  http://dx.doi.org/10.1088/1475-7516/2010/11/042Non-relativistic
  effective theory of dark matter direct detection,  JCAP 1011
  (2010) 042, [http://arxiv.org/abs/1008.15911008.1591].
  
  
Fitzpatrick:2012ix
A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers, and Y. Xu,
  http://dx.doi.org/10.1088/1475-7516/2013/02/004The Effective
  Field Theory of Dark Matter Direct Detection,  JCAP 1302
  (2013) 004, [http://arxiv.org/abs/1203.35421203.3542].

Anand:2013yka
N. Anand, A. L. Fitzpatrick, and W. C. Haxton,
  http://dx.doi.org/10.1103/PhysRevC.89.065501Weakly interacting
  massive particle-nucleus elastic scattering response,  Phys. Rev.
  C89 (2014), no. 6 065501, [http://arxiv.org/abs/1308.62881308.6288].

Dent:2015zpa
J. B. Dent, L. M. Krauss, J. L. Newstead, and S. Sabharwal,
  http://dx.doi.org/10.1103/PhysRevD.92.063515General analysis of
  direct dark matter detection: From microphysics to observational
  signatures,  Phys. Rev. D92 (2015), no. 6 063515,
  [http://arxiv.org/abs/1505.031171505.03117].

Foot:2004pa
R. Foot, http://dx.doi.org/10.1142/S0218271804006449Mirror
  matter-type dark matter,  Int. J. Mod. Phys. D13 (2004)
  2161–2192, [http://arxiv.org/abs/astro-ph/0407623astro-ph/0407623].

Fornengo:2011sz
N. Fornengo, P. Panci, and M. Regis,
  http://dx.doi.org/10.1103/PhysRevD.84.115002Long-Range Forces
  in Direct Dark Matter Searches,  Phys. Rev. D84 (2011) 115002,
  [http://arxiv.org/abs/1108.46611108.4661].

Chu:2011be
X. Chu, T. Hambye, and M. H. G. Tytgat,
  http://dx.doi.org/10.1088/1475-7516/2012/05/034The Four Basic
  Ways of Creating Dark Matter Through a Portal,  JCAP 1205
  (2012) 034, [http://arxiv.org/abs/1112.04931112.0493].

Hooper:2012cw
D. Hooper, N. Weiner, and W. Xue,
  http://dx.doi.org/10.1103/PhysRevD.86.056009Dark Forces and
  Light Dark Matter,  Phys. Rev. D86 (2012) 056009,
  [http://arxiv.org/abs/1206.29291206.2929].

Kaplinghat:2013yxa
M. Kaplinghat, S. Tulin, and H.-B. Yu,
  http://dx.doi.org/10.1103/PhysRevD.89.035009Direct Detection
  Portals for Self-interacting Dark Matter,  Phys. Rev. D89
  (2014), no. 3 035009, [http://arxiv.org/abs/1310.79451310.7945].

Li:2014vza
T. Li, S. Miao, and Y.-F. Zhou,
  http://dx.doi.org/10.1088/1475-7516/2015/03/032Light mediators
  in dark matter direct detections,  JCAP 1503 (2015), no. 03
  032, [http://arxiv.org/abs/1412.62201412.6220].

deLaix:1995vi
A. A. de Laix, R. J. Scherrer, and R. K. Schaefer,
  http://dx.doi.org/10.1086/176322Constraints of selfinteracting
  dark matter,  Astrophys. J. 452 (1995) 495,
  [http://arxiv.org/abs/astro-ph/9502087astro-ph/9502087].

Spergel:1999mh
D. N. Spergel and P. J. Steinhardt,
  http://dx.doi.org/10.1103/PhysRevLett.84.3760Observational
  evidence for selfinteracting cold dark matter,  Phys. Rev. Lett.
  84 (2000) 3760–3763,
  [http://arxiv.org/abs/astro-ph/9909386astro-ph/9909386].

Ackerman:mha
L. Ackerman, M. R. Buckley, S. M. Carroll, and M. Kamionkowski,
  http://dx.doi.org/10.1103/PhysRevD.79.023519,
  10.1142/9789814293792_0021Dark Matter and Dark Radiation,  Phys. Rev. D79 (2009) 023519,
  [http://arxiv.org/abs/0810.51260810.5126]. [,277(2008)].

Feng:2009mn
J. L. Feng, M. Kaplinghat, H. Tu, and H.-B. Yu,
  http://dx.doi.org/10.1088/1475-7516/2009/07/004Hidden Charged
  Dark Matter,  JCAP 0907 (2009) 004,
  [http://arxiv.org/abs/0905.30390905.3039].

Buckley:2009in
M. R. Buckley and P. J. Fox,
  http://dx.doi.org/10.1103/PhysRevD.81.083522Dark Matter
  Self-Interactions and Light Force Carriers,  Phys. Rev. D81
  (2010) 083522, [http://arxiv.org/abs/0911.38980911.3898].

Feng:2009hw
J. L. Feng, M. Kaplinghat, and H.-B. Yu,
  http://dx.doi.org/10.1103/PhysRevLett.104.151301Halo Shape and
  Relic Density Exclusions of Sommerfeld-Enhanced Dark Matter Explanations of
  Cosmic Ray Excesses,  Phys. Rev. Lett. 104 (2010) 151301,
  [http://arxiv.org/abs/0911.04220911.0422].

Loeb:2010gj
A. Loeb and N. Weiner,
  http://dx.doi.org/10.1103/PhysRevLett.106.171302Cores in Dwarf
  Galaxies from Dark Matter with a Yukawa Potential,  Phys. Rev. Lett.
  106 (2011) 171302, [http://arxiv.org/abs/1011.63741011.6374].

Aarssen:2012fx
L. G. van den Aarssen, T. Bringmann, and C. Pfrommer,
  http://dx.doi.org/10.1103/PhysRevLett.109.231301Is dark matter
  with long-range interactions a solution to all small-scale problems of
  Λ CDM cosmology?,  Phys. Rev. Lett. 109 (2012) 231301,
  [http://arxiv.org/abs/1205.58091205.5809].

Tulin:2012wi
S. Tulin, H.-B. Yu, and K. M. Zurek,
  http://dx.doi.org/10.1103/PhysRevLett.110.111301Resonant Dark
  Forces and Small Scale Structure,  Phys. Rev. Lett. 110
  (2013), no. 11 111301, [http://arxiv.org/abs/1210.09001210.0900].

Tulin:2013teo
S. Tulin, H.-B. Yu, and K. M. Zurek,
  http://dx.doi.org/10.1103/PhysRevD.87.115007Beyond
  Collisionless Dark Matter: Particle Physics Dynamics for Dark Matter Halo
  Structure,  Phys. Rev. D87 (2013), no. 11 115007,
  [http://arxiv.org/abs/1302.38981302.3898].

Kaplinghat:2015aga
M. Kaplinghat, S. Tulin, and H.-B. Yu,
  http://dx.doi.org/10.1103/PhysRevLett.116.041302Dark Matter
  Halos as Particle Colliders: Unified Solution to Small-Scale Structure
  Puzzles from Dwarfs to Clusters,  Phys. Rev. Lett. 116 (2016),
  no. 4 041302, [http://arxiv.org/abs/1508.033391508.03339].

Bringmann:2016din
T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg, and P. Walia,
  http://dx.doi.org/10.1103/PhysRevLett.118.141802Strong
  constraints on self-interacting dark matter with light mediators,  Phys. Rev. Lett. 118 (2017), no. 14 141802,
  [http://arxiv.org/abs/1612.008451612.00845].

Kahlhoefer:2017umn
F. Kahlhoefer, K. Schmidt-Hoberg, and S. Wild, Dark matter
  self-interactions from a general spin-0 mediator,
  http://arxiv.org/abs/1704.021491704.02149.

Zavala:2012us
J. Zavala, M. Vogelsberger, and M. G. Walker,
  http://dx.doi.org/10.1093/mnrasl/sls053Constraining
  Self-Interacting Dark Matter with the Milky Way's dwarf spheroidals,  Mon. Not. Roy. Astron. Soc. 431 (2013) L20–L24,
  [http://arxiv.org/abs/1211.64261211.6426].

Vogelsberger:2012ku
M. Vogelsberger, J. Zavala, and A. Loeb,
  http://dx.doi.org/10.1111/j.1365-2966.2012.21182.xSubhaloes in
  Self-Interacting Galactic Dark Matter Haloes,  Mon. Not. Roy. Astron.
  Soc. 423 (2012) 3740, [http://arxiv.org/abs/1201.58921201.5892].

Buckley:2014hja
M. R. Buckley, J. Zavala, F.-Y. Cyr-Racine, K. Sigurdson, and M. Vogelsberger,
  http://dx.doi.org/10.1103/PhysRevD.90.043524Scattering,
  Damping, and Acoustic Oscillations: Simulating the Structure of Dark Matter
  Halos with Relativistic Force Carriers,  Phys. Rev. D90
  (2014), no. 4 043524, [http://arxiv.org/abs/1405.20751405.2075].

Hewett:2012ns
J. L. Hewett et al.,
  http://inspirehep.net/record/1114323/files/arXiv:1205.2671.pdfFundamental Physics at the Intensity Frontier,  2012.
http://arxiv.org/abs/1205.26711205.2671.

Dolan:2014ska
M. J. Dolan, F. Kahlhoefer, C. McCabe, and K. Schmidt-Hoberg,
  http://dx.doi.org/10.1007/JHEP07(2015)103,
  10.1007/JHEP03(2015)171A taste of dark matter: Flavour constraints on
  pseudoscalar mediators,  JHEP 03 (2015) 171,
  [http://arxiv.org/abs/1412.51741412.5174]. [Erratum:
  JHEP07,103(2015)].

Krnjaic:2015mbs
G. Krnjaic, http://dx.doi.org/10.1103/PhysRevD.94.073009Probing
  Light Thermal Dark-Matter With a Higgs Portal Mediator,  Phys. Rev.
  D94 (2016), no. 7 073009, [http://arxiv.org/abs/1512.041191512.04119].

Agnese:2015nto
SuperCDMS, R. Agnese et al.,
  http://dx.doi.org/10.1103/PhysRevLett.116.071301New Results
  from the Search for Low-Mass Weakly Interacting Massive Particles with the
  CDMS Low Ionization Threshold Experiment,  Phys. Rev. Lett. 116 (2016), no. 7 071301, [http://arxiv.org/abs/1509.024481509.02448].

CDMSlitedata
CDMSlite online date release,
  <http://cdms.berkeley.edu/data_releases/CDMSlite_Run1/>.

Agnese:2016cpb
SuperCDMS, R. Agnese et al.,
  http://dx.doi.org/10.1103/PhysRevD.95.082002Projected
  Sensitivity of the SuperCDMS SNOLAB experiment,  Phys. Rev. D95 (2017), no. 8 082002, [http://arxiv.org/abs/1610.000061610.00006].

Angloher:2015ewa
CRESST, G. Angloher et al.,
  http://dx.doi.org/10.1140/epjc/s10052-016-3877-3Results on
  light dark matter particles with a low-threshold CRESST-II detector,  Eur. Phys. J. C76 (2016), no. 1 25,
  [http://arxiv.org/abs/1509.015151509.01515].

Angloher:2015eza
CRESST, G. Angloher et al., Probing low WIMP masses with the next
  generation of CRESST detector,  http://arxiv.org/abs/1503.080651503.08065.

Angloher:2017zkf
CRESST, G. Angloher et al., Description of CRESST-II data,
  http://arxiv.org/abs/1701.081571701.08157.

Strauss:2016sxp
R. Strauss et al., http://dx.doi.org/10.1088/1742-6596/718/4/042048The CRESST-III low-mass WIMP detector,  J. Phys. Conf. Ser. 718 (2016), no. 4 042048.

Armengaud:2016cvl
EDELWEISS, E. Armengaud et al.,
  http://dx.doi.org/10.1088/1475-7516/2016/05/019Constraints on
  low-mass WIMPs from the EDELWEISS-III dark matter search,  JCAP 1605 (2016), no. 05 019, [http://arxiv.org/abs/1603.051201603.05120].

Hehn:2016nll
EDELWEISS, L. Hehn et al.,
  http://dx.doi.org/10.1140/epjc/s10052-016-4388-yImproved
  EDELWEISS-III sensitivity for low-mass WIMPs using a profile likelihood
  approach,  Eur. Phys. J. C76 (2016), no. 10 548,
  [http://arxiv.org/abs/1607.033671607.03367].

Arnaud:2017usi
EDELWEISS, Q. Arnaud et al., Optimizing EDELWEISS detectors for
  low-mass WIMP searches,  http://arxiv.org/abs/1707.043081707.04308.

Peter:2013aha
A. H. G. Peter, V. Gluscevic, A. M. Green, B. J. Kavanagh, and S. K. Lee,
  http://dx.doi.org/10.1016/j.dark.2014.10.006WIMP physics with
  ensembles of direct-detection experiments,  Phys. Dark Univ. 5-6 (2014) 45–74, [http://arxiv.org/abs/1310.70391310.7039].

Cooley:2014aya
J. Cooley, http://dx.doi.org/10.1016/j.dark.2014.10.005Overview
  of Non-Liquid Noble Direct Detection Dark Matter Experiments,  Phys.
  Dark Univ. 4 (2014) 92–97,
  [http://arxiv.org/abs/1410.49601410.4960].

An:2014twa
H. An, M. Pospelov, J. Pradler, and A. Ritz,
  http://dx.doi.org/10.1016/j.physletb.2015.06.018Direct
  Detection Constraints on Dark Photon Dark Matter,  Phys. Lett. B747 (2015) 331–338, [http://arxiv.org/abs/1412.83781412.8378].

Gelmini:2016emn
G. B. Gelmini, Light WIMPs,
  http://arxiv.org/abs/1612.091371612.09137.

Vogelsberger:2012sa
M. Vogelsberger and J. Zavala,
  http://dx.doi.org/10.1093/mnras/sts712Direct detection of
  self-interacting dark matter,  Mon. Not. Roy. Astron. Soc. 430
  (2013) 1722–1735, [http://arxiv.org/abs/1211.13771211.1377].

Chen:2015bwa
C.-S. Chen, G.-L. Lin, and Y.-H. Lin,
  http://dx.doi.org/10.1088/1475-7516/2016/01/013Complementary
  Test of the Dark Matter Self-Interaction by Direct and Indirect Detections,
   JCAP 1601 (2016) 013,
  [http://arxiv.org/abs/1505.037811505.03781].

DelNobile:2015uua
E. Del Nobile, M. Kaplinghat, and H.-B. Yu,
  http://dx.doi.org/10.1088/1475-7516/2015/10/055Direct Detection
  Signatures of Self-Interacting Dark Matter with a Light Mediator,  JCAP 1510 (2015), no. 10 055,
  [http://arxiv.org/abs/1507.040071507.04007].

Strigari:2009zb
L. E. Strigari and R. Trotta,
  http://dx.doi.org/10.1088/1475-7516/2009/11/019Reconstructing
  WIMP Properties in Direct Detection Experiments Including Galactic Dark
  Matter Distribution Uncertainties,  JCAP 0911 (2009) 019,
  [http://arxiv.org/abs/0906.53610906.5361].

Pato:2011de
M. Pato, http://dx.doi.org/10.1088/1475-7516/2011/10/035What
  can(not) be measured with ton-scale dark matter direct detection
  experiments,  JCAP 1110 (2011) 035,
  [http://arxiv.org/abs/1106.07431106.0743].

Kavanagh:2013wba
B. J. Kavanagh and A. M. Green,
  http://dx.doi.org/10.1103/PhysRevLett.111.031302Model
  independent determination of the dark matter mass from direct detection
  experiments,  Phys. Rev. Lett. 111 (2013), no. 3 031302,
  [http://arxiv.org/abs/1303.68681303.6868].

Feldstein:2014gza
B. Feldstein and F. Kahlhoefer,
  http://dx.doi.org/10.1088/1475-7516/2014/08/065A new
  halo-independent approach to dark matter direct detection analysis,  JCAP 1408 (2014) 065, [http://arxiv.org/abs/1403.46061403.4606].

Gluscevic:2015sqa
V. Gluscevic, M. I. Gresham, S. D. McDermott, A. H. G. Peter, and K. M. Zurek,
  http://dx.doi.org/10.1088/1475-7516/2015/12/057Identifying the
  Theory of Dark Matter with Direct Detection,  JCAP 1512
  (2015), no. 12 057, [http://arxiv.org/abs/1506.044541506.04454].

Kahlhoefer:2016eds
F. Kahlhoefer and S. Wild,
  http://dx.doi.org/10.1088/1475-7516/2016/10/032Studying
  generalised dark matter interactions with extended halo-independent
  methods,  JCAP 1610 (2016), no. 10 032,
  [http://arxiv.org/abs/1607.044181607.04418].

Rogers:2016jrx
H. Rogers, D. G. Cerdeno, P. Cushman, F. Livet, and V. Mandic,
  http://dx.doi.org/10.1103/PhysRevD.95.082003Multidimensional
  effective field theory analysis for direct detection of dark matter,  Phys. Rev. D95 (2017), no. 8 082003,
  [http://arxiv.org/abs/1612.090381612.09038].

Roszkowski:2016bhs
L. Roszkowski, E. M. Sessolo, S. Trojanowski, and A. J. Williams,
  http://dx.doi.org/10.1088/1475-7516/2016/08/033Reconstructing
  WIMP properties through an interplay of signal measurements in direct
  detection, Fermi-LAT, and CTA searches for dark matter,  JCAP 1608 (2016), no. 08 033, [http://arxiv.org/abs/1603.065191603.06519].

Roszkowski:2017dou
L. Roszkowski, S. Trojanowski, and K. Turzynski, Towards understanding
  thermal history of the Universe through direct and indirect detection of dark
  matter,  http://arxiv.org/abs/1703.008411703.00841.

Kouvaris:2014uoa
C. Kouvaris, I. M. Shoemaker, and K. Tuominen,
  http://dx.doi.org/10.1103/PhysRevD.91.043519Self-Interacting
  Dark Matter through the Higgs Portal,  Phys. Rev. D91 (2015),
  no. 4 043519, [http://arxiv.org/abs/1411.37301411.3730].

Bernal:2015ova
N. Bernal, X. Chu, C. Garcia-Cely, T. Hambye, and B. Zaldivar,
  http://dx.doi.org/10.1088/1475-7516/2016/03/018Production
  Regimes for Self-Interacting Dark Matter,  JCAP 1603 (2016),
  no. 03 018, [http://arxiv.org/abs/1510.080631510.08063].

Lewin:1995rx
J. D. Lewin and P. F. Smith,
  http://dx.doi.org/10.1016/S0927-6505(96)00047-3Review of
  mathematics, numerical factors, and corrections for dark matter experiments
  based on elastic nuclear recoil,  Astropart. Phys. 6 (1996)
  87–112.

Aprile:2017iyp
XENON, E. Aprile et al., First Dark Matter Search Results from the
  XENON1T Experiment,  http://arxiv.org/abs/1705.066551705.06655.

Akerib:2015cja
LZ, D. S. Akerib et al., LUX-ZEPLIN (LZ) Conceptual Design
  Report,  http://arxiv.org/abs/1509.029101509.02910.

Feldstein:2014ufa
  B. Feldstein and F. Kahlhoefer,
  http://dx.doi.org/10.1088/1475-7516/2014/12/052Quantifying (dis)agreement between direct detection experiments in a halo-independent way,  
  JCAP 1412 (2014) no.12,  052, 
  [http://arxiv.org/abs/1409.54461409.5446].
  
Frandsen:2011cg
M. T. Frandsen, F. Kahlhoefer, S. Sarkar, and K. Schmidt-Hoberg,
  http://dx.doi.org/10.1007/JHEP09(2011)128Direct detection of
  dark matter in models with a light Z',  JHEP 09 (2011) 128,
  [http://arxiv.org/abs/1107.21181107.2118].

Zheng:2014nga
H. Zheng, Z. Zhang, and L.-W. Chen,
  http://dx.doi.org/10.1088/1475-7516/2014/08/011Form Factor
  Effects in the Direct Detection of Isospin-Violating Dark Matter,  JCAP 1408 (2014) 011, [http://arxiv.org/abs/1403.51341403.5134].

Fox:2010bz
P. J. Fox, J. Liu, and N. Weiner,
  http://dx.doi.org/10.1103/PhysRevD.83.103514Integrating Out
  Astrophysical Uncertainties,  Phys. Rev. D83 (2011) 103514,
  [http://arxiv.org/abs/1011.19151011.1915].

Frandsen:2011gi
M. T. Frandsen, F. Kahlhoefer, C. McCabe, S. Sarkar, and K. Schmidt-Hoberg,
  http://dx.doi.org/10.1088/1475-7516/2012/01/024Resolving
  astrophysical uncertainties in dark matter direct detection,  JCAP
  1201 (2012) 024, [http://arxiv.org/abs/1111.02921111.0292].

Gondolo:2012rs
P. Gondolo and G. B. Gelmini,
  http://dx.doi.org/10.1088/1475-7516/2012/12/015Halo independent
  comparison of direct dark matter detection data,  JCAP 1212
  (2012) 015, [http://arxiv.org/abs/1202.63591202.6359].

Cherry:2014wia
J. F. Cherry, M. T. Frandsen, and I. M. Shoemaker,
  http://dx.doi.org/10.1088/1475-7516/2014/10/022Halo Independent
  Direct Detection of Momentum-Dependent Dark Matter,  JCAP 1410
  (2014), no. 10 022, [http://arxiv.org/abs/1405.14201405.1420].

Frandsen:2013cna
M. T. Frandsen, F. Kahlhoefer, C. McCabe, S. Sarkar, and K. Schmidt-Hoberg,
  http://dx.doi.org/10.1088/1475-7516/2013/07/023The unbearable
  lightness of being: CDMS versus XENON,  JCAP 1307 (2013) 023,
  [http://arxiv.org/abs/1304.60661304.6066].

Bozorgnia:2016ogo
N. Bozorgnia, F. Calore, M. Schaller, M. Lovell, G. Bertone, et al.,
  http://dx.doi.org/10.1088/1475-7516/2016/05/024Simulated Milky
  Way analogues: implications for dark matter direct searches,  JCAP
  1605 (2016), no. 05 024, [http://arxiv.org/abs/1601.047071601.04707].

Kelso:2016qqj
C. Kelso, C. Savage, M. Valluri, K. Freese, G. S. Stinson, et al.,
  http://dx.doi.org/10.1088/1475-7516/2016/08/071The impact of
  baryons on the direct detection of dark matter,  JCAP 1608
  (2016) 071, [http://arxiv.org/abs/1601.047251601.04725].

Green:2017odb
A. M. Green, http://dx.doi.org/10.1088/1361-6471/aa7819Astrophysical uncertainties on the local dark matter distribution and direct
  detection experiments,  J. Phys. G44 (2017), no. 8 084001,
  [http://arxiv.org/abs/1703.101021703.10102].

Bellazzini:2013foa
B. Bellazzini, M. Cliche, and P. Tanedo,
  http://dx.doi.org/10.1103/PhysRevD.88.083506Effective theory of
  self-interacting dark matter,  Phys. Rev. D88 (2013), no. 8
  083506, [http://arxiv.org/abs/1307.11291307.1129].

Weinberg:2013aya
D. H. Weinberg, J. S. Bullock, F. Governato, R. Kuzio de Naray, and A. H. G.
  Peter, http://dx.doi.org/10.1073/pnas.1308716112Cold dark
  matter: controversies on small scales,  Proc. Nat. Acad. Sci. 112 (2014) 12249–12255, [http://arxiv.org/abs/1306.09131306.0913].

Tulin:2017ara
S. Tulin and H.-B. Yu, Dark Matter Self-interactions and Small Scale
  Structure,  http://arxiv.org/abs/1705.023581705.02358.

Sebitobepublished
M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, 2017.
In preparation.

PICO
PICO, O. Harris, Dark matter searches with PICO bubble chambers:
  An overview,  in APS April Meeting, 2017.

Yellin:2002xd
S. Yellin, http://dx.doi.org/10.1103/PhysRevD.66.032005Finding an
  upper limit in the presence of unknown background,  Phys. Rev. D66 (2002) 032005, [http://arxiv.org/abs/physics/0203002physics/0203002].

Akerib:2015rjg
LUX, D. S. Akerib et al.,
  http://dx.doi.org/10.1103/PhysRevLett.116.161301Improved Limits
  on Scattering of Weakly Interacting Massive Particles from Reanalysis of 2013
  LUX Data,  Phys. Rev. Lett. 116 (2016), no. 16 161301,
  [http://arxiv.org/abs/1512.035061512.03506].



