Direct Detection Signatures of Self-Interacting Dark Matter with a Light Mediator

Self-interacting dark matter (SIDM) is a simple and well-motivated scenario that could explain long-standing puzzles in structure formation on small scales. If the required self-interaction arises through a light mediator (with mass $\sim 10$ MeV) in the dark sector, this new particle must be unstable to avoid overclosing the universe. The decay of the light mediator could happen due to a weak coupling of the hidden and visible sectors, providing new signatures for direct detection experiments. The SIDM nuclear recoil spectrum is more peaked towards low energies compared to the usual case of contact interactions, because the mediator mass is comparable to the momentum transfer of nuclear recoils. We show that the SIDM signal could be distinguished from that of DM particles with contact interactions by considering the time-average energy spectrum in experiments employing different target materials, or the average and modulated spectra in a single experiment. Using current limits from LUX and SuperCDMS, we also derive strong bounds on the mixing parameter between hidden and visible sector.


Introduction
Dark matter (DM) accounts for 85% of the matter in the universe [1], but its particle physics nature remains elusive. In one of the well-studied models, DM candidates have a weakscale mass and carry weak-scale interactions. In this weakly interacting massive particle (WIMP) paradigm, the interaction strength between DM and the standard model (SM) sector is also weak-scale and the resulting relic density for WIMPs could be consistent with, or smaller than, the observed cosmological DM density. This coincidence of weak-scale and relic density has motivated many direct, indirect, and collider searches for WIMPs leading to significant regions of the model parameter space being ruled out. Alternatively, axions motivated by the strong CP problem could have the right relic density to be all of the DM and present laboratory searches and limits from stellar evolution have ruled out large regions of parameter space.
Both WIMPs and axions are consistent with the prevailing view in astrophysics that DM is a cold non-interacting particle. The cold DM paradigm has been thoroughly tested on large scales. On galactic and sub-galactic scales, however, the data strongly suggest that our ideas of structure formation are incomplete. In particular, the DM densities inferred in the central regions of DM-dominated galaxies are, on average, lower than expected from DM-only simulations [2,3]. Given this situation and the lack of a consistent signal in WIMP searches, it is important to thoroughly study DM candidates other than WIMPs. A compelling possibility is that DM may strongly interact with itself [4] (self-interacting dark matter or SIDM). DM self-interactions can transfer energy from the hotter, outer region to the inner region of DM halos, thereby reducing the central density and providing a possible solution to the smallscale puzzles. Recent numerical simulations have shown that SIDM halos are consistent with observations if DM particles scatter with each other with a nuclear-scale cross section within halos [2,[5][6][7][8][9][10].
The key ingredient of the majority of SIDM models is the existence of a light mediator with mass 100 MeV [11][12][13][14][15][16][17], which must decay in the early universe to avoid overclosure [17][18][19]. Unless additional massless degrees of freedom are introduced, the mediator naturally decays to SM particles through a mixing portal between the two sectors, which may lead to SIDM signals in indirect and direct detection experiments [18]. For example, -1 -

JCAP10(2015)055
SIDM particles may annihilate to mediators, which subsequently decay to electrons. These energetic electrons could generate gamma-ray signals in environments such as the galactic center through inverse Compton scattering on starlight [20]. It has also been shown that DM direct detection experiments are sensitive to the SIDM parameter space even if the mixing parameter between the two sectors is extremely feeble [18,[21][22][23].
An interesting feature of SIDM in direct detection is that DM interacts with nuclei with a long-range force, in contrast to a contact interaction in WIMP models. The usual tree-level t-channel DM-nucleus differential scattering cross section scales as (q 2 + m 2 φ ) −2 , where q is the momentum transfer in nuclear recoils and m φ is the mediator mass. The typical mediator mass in SIDM models is comparable to the typical momentum transfer in direct detection experiments and hence the SIDM event spectrum is more peaked toward low recoil energies, compared to WIMP models where m φ q. This novel feature of SIDM signals in direct detection was noted in ref. [18], where an approximation method was used to reinterpret direct detection bounds for WIMPs to constrain SIDM models. The effect of light mediators on direct detection bounds was also discussed recently in ref. [23] (see also [24,25]), with a focus on the bound dependence on the mediator mass and the determination of the DM mass from experimental data. In this paper we investigate direct detection of DM with light mediators in detail, paying special attention to their scattering spectra. We first derive a strong upper bound on the mixing parameter between the SIDM and SM sectors. Then we study the SIDM event spectrum in the light of SuperCDMS, LUX, and DAMA, taking into account realistic efficiency and energy resolution of the detectors. Considering both the time-average and modulated event rates, we show that direct detection experiments can potentially distinguish SIDM from WIMPs.
In the remainder of this work, we first present a simple particle physics model for SIDM, and discuss basics of DM direct detection in section 2. Our results are presented in section 3. Lastly, we conclude in section 4.

Particle physics model
We assume that the DM particle X, either a Dirac fermion or a complex scalar, interacts with the vector mediator φ of a dark U(1) X gauge interaction. In the non-relativistic limit, self-interactions between DM particles can be described by a Yukawa potential [11,[26][27][28][29][30][31][32] where α X ≡ g X /(4π) is the fine structure constant in the dark sector and m φ is the mediator mass. We fix α X = 0.01, motivated by the value of the electromagnetic fine structure constant in the SM. We also focus on the case of asymmetric SIDM in which only DM X, and not its anti-particle, is present in DM halos. Hence, DM self-scattering is purely repulsive and the "+" sign of the Yukawa interaction in eq. (2.1) must be considered. In general, the dark sector can couple to the SM through the kinetic mixing γ φ µν F µν [33], where γ is the mixing parameter, and φ µν and F µν are the field strength of the mediator φ and of the photon, respectively. The mixing induces a coupling of φ to SM fermions f at O( γ ) upon diagonalization: γ e f Q ff γ µ f φ µ , where Q f denotes the electric charge (in units of e) of the SM fermions. In this case, direct detection signals of SIDM arise from DM-proton scattering via φ exchange. Our analysis can be easily generalized to other cases (MeV)  A  1000  3  127  74  B  100  15  62  46  C  10  20  10  10  D  5  20  5  5   Table 1. SIDM benchmark models considered in this paper. In the two rightmost columns we indicate the typical values of the momentum transfer for recoils off xenon (relevant for LUX) and germanium (relevant for SuperCDMS) of a DM particle with typical speed in Earth's frame v = 232 km/s. The maximum attainable momentum transfer, assuming a maximum DM speed v max = v esc + v = 776 km/s, is 4.7 times higher than this typical value, while values that are lower than those shown here are always possible. such as φ-Z mixing portal, or Higgs portal for a scalar mediator [18]. Notice however that all these models have similar phenomenology at direct detection experiments; the main difference being that for the kinetic mixing case the DM interacts dominantly with protons, for the Zmixing case mostly with neutrons, and for the Higgs portal case equally with protons and neutrons.
The differential cross section for DM-nucleus scattering is [18,24] dσ XT where α em = 1/137 is the SM fine structure constant, Z is the number of protons in the nucleus, q is the momentum transfer, v is the speed of the DM particle in the nucleus rest frame and F T (q 2 ) is the nuclear form factor related to the charge density in the nucleus [34,35]. The nuclear recoil energy E R is related to the momentum transfer and the nuclear mass To investigate the signal spectrum of SIDM in SIDM-nucleus scattering, we choose four benchmark models as shown in table 1. Also shown in table 1 are the typical values of the momentum transfer for recoils off xenon (relevant for LUX) and germanium (relevant for SuperCDMS) of a DM particle with typical speed in Earth's frame v = 232 km/s, q ≈ √ 2µ T v with µ T the DM-nucleus reduced mass. While smaller values of the momentum transfer are always possible depending on the scattering angle, the maximum attainable momentum transfer is 2µ T v max , with v max the maximum possible DM speed in Earth's frame (discussed later). We can anticipate for which models the long-range nature of the interaction will be important by comparing m φ with the typical q. Figure 1 shows the DM self-scattering cross section of our four benchmark models as a function of the DM relative speed. In model A, the self-scattering cross section is suppressed significantly at large velocities, because DM self-scattering occurs in the Rutherford limit with σ XX ∝ 1/v 4 on large scales. For model B, DM self-interactions are important in dwarf galaxies, and mildly in Milky Way-sized galaxies, but are suppressed on cluster scales. On the other hand, σ XX /m X is relevant from dwarf to cluster scales for both model C and D. Repulsive self-interaction is assumed. Curves correspond to the SIDM models A, B, C, D summarized in Table 1.

Scattering rate in direct detection experiments
The differential recoil rate for a DM particle scattering off a target nucleus T , expressed in counts per day per kilogram per keV, is where ξ T /m T is the number of targets per detector mass in units of kg −1 , f (v, t) is the DM velocity distribution in Earth's frame and ρ is the local DM density, which we set to 0.3 GeV/cm 3 . v min (E R ) is the minimum speed a DM particle must have to impart a recoil energy E R to the target nucleus; for elastic scattering, While no information is available at present on the DM velocity distribution in our galaxy, we notice that DM self-interactions would thermalize the halo, and therefore it is justified to consider a thermal velocity distribution [36]. We assume the Standard Halo Model (SHM) for the DM halo, i.e. an isothermal sphere with an isotropic Maxwellian velocity distribution in the galactic frame: π a normalization factor that ensures that d 3 v f G (v) = 1. The distribution is truncated at the galactic escape speed v esc = 544 km/s [37], and we set v 0 equal to the speed of the Local Standard of Rest, v 0 = 220 km/s. The DM velocity distribution f (v, t) in Earth's frame is related to f G (v) by the Galilean transformation f (v, t) = f G (v + v obs (t)), with v obs (t) the velocity of Earth with respect to JCAP10(2015)055 the galactic frame. This is v obs (t) = v + v ⊕ (t), where v ⊕ (t) is Earth's rotational velocity around the Sun and v is the Sun's velocity with respect to the galactic frame. The average Earth's speed over a revolution period of one year equals | v obs | = v = 232 km/s. Notice that, because of the Galilean boost, the maximum speed attainable in Earth's frame is (in Given the dependency of the differential cross section in eq. (2.2) on the DM speed, dσ XT /dE R ∝ 1/v 2 , which is common to many non-relativistic DM-nucleus scattering processes, the relevant velocity integral entering the differential rate in eq. (2.3) is Notice that the DM velocity distribution being normalized implies that also η(v min , t) is at any time t. This can be seen by denoting /v, and then inverting the integration order.
If the DM halo is smooth and isotropic, 1 since v , v 0 v ⊕ = 30 km/s we can Taylorexpand the rate in eq. (2.3) in the dimensionless parameter and v are aligned, ω = 2π/yr and b ≈ 0.5 (see e.g. [42]). Therefore we have for the velocity integral in eq. (2.6) where we identified the unmodulated (time-average) and modulated components of the velocity integral, η 0 and η 1 respectively (notice that these are time-independent quantities). While η 0 is positive by definition, η 1 can also assume negative values. Eq. (2.7) implies as can be inferred by taking t = t 0 + π 2 . These equalities hold as long as the approximation in eq. (2.9) is valid. Eq. (2.9) implies that the differential rate depends on time approximately as The factor multiplying cos[ω(t − t 0 )] is the modulated (component of the) rate, while the first term is the average or unmodulated rate, for which we keep the same symbol as the total JCAP10(2015)055 rate, with a slight abuse of notation, since dR mod T /dE R is negligible in comparison. As with η 0 and η 1 , the modulated and the unmodulated rates are time-independent quantities, and dR mod T /dE R can assume negative values.
To actually compare the theoretical rate in eq. (2.3) with the event rate measured by the experiments, we now have to take into account detection efficiency, cuts acceptance and energy resolution of the detectors. We consider LUX and SuperCDMS for the measure of the unmodulated rate, and DAMA for the measure of the modulated rate.
For LUX, we determine the differential rate in S1 (primary scintillation light, measured in photoelectrons phe) following the procedure used by the similar experiment XENON100 [43]: the number n of photons created by a recoiling Xe nucleus in the liquid phase, which is Poisson distributed, is smeared at detection due to the finite resolution of the photomultiplier tubes (PMT) whose effect is parameterized by a Gaussian distribution. The measured differential rate is then where σ PMT = 0.37 is the average single photoelectron resolution of the PMT's [44], and T denotes the different Xe isotopes. Gauss(x|µ, σ) and Poiss(x|µ) denote respectively Gaussian and Poisson probability distributions for x with mean µ and, for the Gaussian, standard deviation σ. For the efficiency (S1) we take the dashed red line in the top panel of figure 1 of ref. [45]. We set a conservative lower threshold of 3 keV on the modeled signal [45]. ν(E R ), the average number of photons generated by a Xe nucleus recoiling with energy E R , is obtained from the spreadsheet available on the NEST website [46] using as input 2.888 g/cm 3 for the density of liquid xenon and 181 V/cm for the drift field [47,48]. The energy range where the signal S1 is measured is 2-30 phe. For SuperCDMS, the measured differential recoil rate in detected energy E is with the detector resolution σ(E ) = 0.293 2 + 0.056 2 E /keV adopted for the CDMS-II germanium detectors [49], and the efficiency (E R ) taken to be the red line in figure 1 of ref. [50]. T here denotes the different Ge isotopes. The energy range where the signal E is measured is 1.6-10 keV.
We also study the annual modulation signal of SIDM, taking the DAMA experiment as an example [51]. Our study can also be applied straightforwardly to future DAMAlike experiments as KIMS-NaI, ANAIS, DM-Ice17, and SABRE (see e.g. [52] and references therein). DAMA collected events in the energy range 2-20 keV for an exposure of 1.33 ton · yr. The Collaboration measures the modulated rate without attempting any background subtraction. The measured modulated rate is with the detector resolution σ(E) = 0.448 √ E + 0.0091E [53]. The sum over target nuclei counts sodium and iodine, with quenching factors Q Na = 0.3 and Q I = 0.09 respectively. 3 Direct detection of self-interacting dark matter

Current constraints from LUX and SuperCDMS
In the first LUX results [45], the collected data were found consistent with the backgroundonly hypothesis. With an exposure of ω = 85.3 day · 118.3 kg, one single candidate event was found below the mean of the nuclear recoil distribution in the S1-log 10 (S2/S1) plane (solid red line in figures 3 and 4 of ref. [45]), where 0.64 background events were estimated. We consider the DM signal in the same region by assuming a 50% efficiency for DM events below the nuclear recoil mean for heavy enough DM particles. We therefore compute the number of signal events as 0.50 × w To set a constraint on the kinetic mixing parameter γ for SIDM models, we require the p-value of the data to be smaller than the significance level α = 10%, assuming a Poisson distribution for the total number of events. This leads to an upper limit of 3.25 events on the total number of signal events in LUX.
In SuperCDMS, a blind analysis with an exposure 577 kg · day revealed 11 candidate events against an estimated background of 6.1 events [50]. We derive an upper bound on γ as above, yielding an upper bound of 10.5 signal events in SuperCDMS.
In deriving our bounds, we only use the total number of events in the signal range instead of using the full spectral information of the rate. While a spectral analysis is ultimately needed in order to distinguish SIDM from WIMPs (see next section), extracting information on the particle physics model of DM interactions relies on assuming a model for the DM distribution. The SHM, used in this work, although motivated in the framework of SIDM, is only a first guess. The effect of baryons, anisotropies in the DM distribution and DM substructure will affect the velocity distribution. Our bounds are therefore somewhat solid in this respect, in the sense that they allow for small variations of the DM distribution, that can modify the detected event spectrum, but such that the total number of events remains fixed. Figure 2 shows direct detection constraints on the SIDM parameter space. The blue lines in the left panel denote the portion of parameter space where SIDM could explain the small scale anomalies, for three different values of the self-interaction cross section per unit mass σ XX /m X = 0.1, 1, 10 cm 2 /g. It is clear that both LUX and SuperCDMS put a strong constraint on the mixing parameter γ . As can be seen in the left panel of figure 2, LUX excludes all favored (m X , m φ ) regions with γ 10 −9 except for m X 7 GeV. The SuperCDMS limit is weaker, but it can exclude SIDM models with m X > 3 GeV. Remarkably, benchmark points A, B and C are ruled out by LUX for γ = 10 −9 , while benchmark point D can not be excluded by LUX because of the small DM mass (see below). It can however be excluded by SuperCDMS for γ ∼ 10 −8 due to its lower energy threshold and its lighter target compared to xenon. It is remarkable that the LUX and SuperCDMS constraints on γ are much stronger than those from beam dump experiments [54], which exclude for instance   dominated by decays into e + e − [18]. If the mediator decays before weak freeze-out, ∼ 1 s, then we can be assured that the predicted light element abundances will not be modified from the standard scenario. From figure 2 (right) we see that direct detection experiments are probing the SIDM parameter space which is relevant to the thermal history of the early universe.
To further narrow down the allowed parameter range for γ 10 −10 , a detailed study of the impact of decays on Big Bang Nucleosynthesis (BBN) is required. It has been shown that BBN provides a model-independent constraint on γ between 10 −14 and 10 −12 , due to the destruction of 3 He and 2 H by the electromagnetic showers resulting from φ decays, with the conservative assumption that mediator particles are produced in the early universe by inverse decays only [55]. This limit needs to be reevaluated for the case where the dark photons were thermally populated in the thermal bath of the hidden sector [56,57].
Since the range 10 −12 γ 10 −10 will be tested by direct detection experiments soon, let us consider BBN constraints on γ in this range in some more detail. If the photons in the shower can scatter off background photons and produce electron-position pairs, then the distribution of photons in the shower is rapidly pushed to energies below ∼ m 2 e /(22T γ ) = 0.012(MeV/T γ ) MeV [58,59], as long as the mediator mass is sufficiently larger than this threshold. As an example, this implies that the bulk of energy in the shower created by decays at T γ = 100 keV (τ φ ∼ 100 s) will be shifted rapidly to energies below the 2 H binding energy, 2.2 MeV. There could, however, be an effect on the light element abundances because the presence of the thermalized mediator (or its decay products) modifies the expansion rate during BBN. The magnitude of this effect depends on the temperature of the hidden -8 -

JCAP10(2015)055
sector, which could be very different from that of the visible sector since kinetic mixing with the low γ values under consideration could not have kept the two sectors in thermal equilibrium [60]. Based on these arguments, there should be regions of SIDM parameter space with 10 −12 γ 10 −10 which are consistent with all the BBN constraints and accessible to direct detection experiments.
It is also suggestive to compare direct detection bounds with other astrophysical constraints. For mediator masses in the range of interest here, m φ 100 MeV, the parameter space with γ 10 −10 is strongly disfavored by supernova cooling arguments (similar to those for axions) [61][62][63]. This constraint becomes ineffective when the coupling is large enough to trap the mediator particles within the supernova core, which happens above γ 10 −7 . Interestingly, the lower limit from the supernova cooling constraint will be superseded by direct detection experiments soon. In fact, for m X ∼ 100 GeV, the current LUX bound, γ 10 −10 , is already comparable to the supernova cooling limit. We now focus our attention back on the bounds in the (m X , m φ ) plane. To get a better understanding of the bounds, it is useful to indicatively write the total average rate in a certain recoil energy interval [E R1 , E R2 ] as (3.1) Here we neglected all numerical constants and normalization factors, apart from γ , since we are only interested in the functional dependence of the rate on the model parameters, and we made explicit the dependence of v min on m X . We also neglected the nuclear form factor as it has no bearing on the following discussion. The presence of multiple isotopes in the detector and the signal smearing due to the finite detector resolution should not affect our argument, either. (2.10)), this also implies a reduction of the velocity integral for small v min . In the bottom panels of figure 3 we show again the velocity integral but now in E R -space, η 0 (v min (E R , m X )), for our four benchmark values of the DM mass, m X = 5, 10, 100, and 1000 GeV. For each value of the DM mass, the solid line corresponds to our standard set of SHM parameters while the dashed line is for the larger values reported above. The variables nuclear recoil energy E R and minimum speed v min can be used interchangeably, and one can switch from one to the other via the m X -dependent change of variables eq. (2.4). The plots in v min -space and in E R -space in figure 3 show the same thing in different manners, and we think it is useful to present them together.
Due to the correspondence between E R and v min , one can think of the rate in eq.
10 100 values, and in E R -space (vertical dashed gray lines, bottom). It is clear that, while heavy DM samples small v min values, light DM samples the high-velocity tail of the distribution: this is expected, since light DM particles need higher speeds to impart a certain recoil energy to the nucleus, compared to heavier DM particles. For small enough m X , the experiment stops being sensitive to DM scattering events because there are no DM particles in the halo that are energetic enough to scatter a target nucleus above the experimental threshold. In other words, the v min range probed by the experiment lies completely above the maximum possible DM speed in Earth's reference frame, v max , and the only way of gaining sensitivity -10 -

JCAP10(2015)055
to such light DM particles is lowering the experimental energy threshold. This is the reason why a 5 GeV DM particle as in our benchmark model D lies outside of LUX's reach (with the 3 keV threshold), while being within the reach of SuperCDMS which employs a lighter target and has a lower energy threshold (see the left panel of figure 2).
In the opposite regime of heavy DM, or more quantitatively m X m T , the minimum speed in eq. (2.4) becomes independent of m X and therefore the v min range probed by the experiment does not change when increasing the DM mass. Thus, an experiment like LUX probes the same range in v min -space for a 1000 GeV DM particle as for heavier DM particles, and the same happens with SuperCDMS. For this reason, the rate eq. (3.1) at these two experiments simply scales as R ∼ 1/m X for m X > 1000 GeV.
We are now ready to understand the dependence of the limits shown in figure 2 on the model parameters. In the simplest case where an experimental result is employed to only set an upper limit R limit on the predicted rate, a bound on the model parameters is set by requiring R limit > R(m X , m φ , γ ). Recalling what we just said about the velocity integral η 0 (v min ), let us first notice that, for γ and m φ fixed, the rate in eq. (3.1) approaches zero for m X → ∞, and vanishes for all DM masses small enough that v min (E R , m X ) > v max for all values of E R probed by the experiment. For intermediate values of the DM mass, the rate is non-zero and it has therefore a global maximum, which is reached in the massless mediator limit (effectively attained in the long-range regime m 2 φ 2m T E R ). Therefore, there exist a value of the mixing parameter γ below which no bounds can be set, because for this value the rate is smaller than R limit in all of the parameter space.
Above this critical γ , bounds can be determined in the following way. Starting from a point in parameter space where m X is very large and m φ = 0 (bottom-right portion of the left panel of figure 2), we can decrease m X keeping the mixing parameter fixed until we obtain R = R limit . Models with DM mass larger than that so determined are viable, while those with similar but smaller m X predict a larger rate and are therefore ruled out by the experiment. For fixed γ , the limit line in parameter space can be determined by progressively 'turning on' m φ ; the limit line raises vertical in the (m X , m φ ) plane until m φ grows of the same order of magnitude as the momentum transfer 2m T E R , when the rate becomes sensitive to the mediator mass and a raise in m φ can be compensated by decreasing m X , so that the predicted rate eq. (3.1) remains equal to the limit rate. We then enter the contact-interaction regime, where m φ 2m T E R and, if the DM mass is still large enough that η 0 (v min (E R , m X )) changes slowly with decreasing m X in the E R range probed by the experiment, the limit line scales in the parameter space as 2 γ /(m X m 4 φ ) = constant, as can be seen from eq. (3.1). If we keep following the limit line for smaller and smaller m X values, eventually the DM mass becomes so small that the velocity integral vanishes quickly because all scattering events occur below the experimental threshold. Such a reduction in the rate can be compensated by a decrease in m φ , if we keep γ fixed, until we enter again the longrange regime where the rate is no longer sensitive to m φ and the limit line drops vertically in the (m X , m φ ) plane. For smaller values of the mixing parameter, however, the condition R = R limit can be obtained even before reaching the kinematical limit of the experiment.
We can now also easily understand the effect of changing the halo parameters or the experimental threshold on the parameter space limits. Figure 3 shows that increasing the typical DM speed makes η 0 (v min ) larger at large v min and smaller at low v min values, i.e. the differential rate increases at high energies and decreases at low energies. Consequently, the total rate in eq. (3.1) will be smaller for heavy DM, which probes small values of v min , while it will be larger for light DM. For fixed γ and m φ 2m T E R , this change can be contrasted -11 -

JCAP10(2015)055
by a decrease in m X , since η 0 (v min (E R , m X ))/m X increases (decreases) with the DM mass for light (heavy) enough DM at fixed E R , so to maintain the limit condition R = R limit . Therefore, a larger DM average speed will cause the limits in figure 2 to move slightly to smaller values of m X . Decreasing the experimental energy threshold has the general effect of increasing the rate, thus the limits on the parameter space will become more stringent. The bounds on the (m X , m φ ) plane will therefore move to larger m X for heavy DM and to smaller m X for light DM.

Identifying SIDM with LUX and SuperCDMS
To see whether DM direct detection experiments are able to distinguish SIDM candidates from usual WIMPs if a positive signal is detected, we perform a careful study of the SIDM signal spectrum in SuperCDMS and LUX. Given the current constraints on the mixing parameter as discussed above, we take γ = 10 −10 for the SIDM models as a baseline. We note that the BBN constraints discussed earlier may require γ 10 −12 and therefore current and next generation experiments should be able to test much of the viable parameter space of this simple SIDM setup. For each of our four SIDM benchmark models, we also consider two additional cases, for comparison: one is with the mediator mass to be three times the value of m φ in the benchmark model; the other is a WIMP model, where a spin-independent (SI) contact interaction between DM and nuclei is assumed. The SI cross section can be taken to be that in eq. (2.2) in the limit m φ q, and with arbitrary values of the coupling constants. We normalize the coupling strength of the two additional models such that they have the same total number of events (proportional to the integral of the measured rate) as the benchmark model. The model with three times larger mediator mass is meant to provide a comparison point between the benchmark SIDM model (exhibiting long-range dynamics due to the light mediator) and the contact interaction of the SI model.
The quantity that is ultimately measured by the experiments is the differential rate in primary scintillation light/photoelectrons for LUX, dR/dS1 in eq. (2.12), or in detected energy for SuperCDMS, dR/dE in eq. (2.13). Therefore, it makes sense to compare these rates for our different DM models, rather than the theoretical recoil rate dR/dE R ≡ T dR T /dE R . To understand how the experiment-dependent effects, as resolution and efficiency, affect the spectrum when the rate in S1 and E is computed, we show in figure 4 the theoretical spectrum (top) and how it changes after the detector resolution has been taken into account (bottom), corresponding to neglecting the experimental efficiency in the measured rate. We choose our benchmark model B for this illustration, with m X = 100 GeV and m φ = 15 MeV (solid red line); the model with three times the mediator mass, m φ = 45 MeV (dashed purple line) is plotted along, together with the SI contact interaction (dotted green line). While theoretical and measured spectra are almost identical for SuperCDMS (before considering ) due to its high resolution, the effect of adding the detector resolution is significant for LUX, in particular in the low E R region. This is due to the involved process of conversion of the nuclear recoil energy into a signal in noble liquid detectors (and possibly due to the conservative 3 keV cut in the integral in eq. (2.12)); our treatment, while lacking a detailed description of such process, already captures some of the modifications occurring to the spectrum. Also plotted in the bottom panels of figure 4 is the S1 range of LUX and the detected energy range of SuperCDMS (dashed vertical lines). For LUX, we also include on the top horizontal axis the E R scale corresponding in average to the S1 values on the bottom axis, defined by E R = ν −1 (S1) (see footnote 2).  . All curves are normalized to yield the same number of measured events when the measured rate is considered. Top: theoretical scattering rate as a function of the recoil energy E R . Bottom: theoretical scattering rate integrated with the resolution function, or in other words, the measured rate prior to including the effect of the experimental efficiency , as a function of the detected signal (S1 for LUX and E for SuperCDMS). Also plotted is the range in detected signal probed by the experiments (dashed vertical lines). For LUX, we also provide on the top axis the average recoil energy E R corresponding to the detected signal S1 in photoelectrons.
From figure 4 it can be noted that the SIDM model has a spectrum that is peaked towards low nuclear recoil energy, compared to the SI contact interaction model, as expected from the long-range nature of the interaction. For its contact nature, instead, the SI rate is virtually energy independent, if not for the drop at high energies due to the nuclear form factor, which parameterizes the loss of coherence experienced when the projectile (in this case the DM particle) stops 'seeing' the nucleus as a whole and starts resolving its internal structure. As anticipated, the model with three times the mediator mass as the SIDM model has a spectrum in between these two cases, and all three spectra are very similar at high energy, when the mediator mass becomes negligible. Figure 5 shows the measured event spectrum at LUX (left) and SuperCDMS (right) for all SIDM benchmark models, each compared to the model with three times the mediator mass and to the SI model with same DM mass. With respect to figure 4, here we include both detector resolution and efficiency. The dashed vertical lines enclose the range in detected signal probed by the experiment. The three curves in each plot are normalized to have the  . Measured rates at LUX (left) and SuperCDMS (right), for different DM masses. For each of our benchmark SIDM models (solid red line), a model with three times the mediator mass (dashed purple line) and a SI model with contact interaction (dotted green line) are also considered. The spectra are normalized to have the same area within the signal range, enclosed by the two vertical dashed lines. For LUX, we also provide on the top axis the average recoil energy E R corresponding to the detected signal S1 in photoelectrons. Notice that a 5 GeV DM particle is below threshold for LUX, and therefore the measured rate is zero.

JCAP10(2015)055
same area within the signal range (same number of measured events), exactly as in figure 4. The experimental efficiency has the general effect of reducing the overall rate; moreover, it is smaller at low energies, as could be guessed from a comparison with figure 4. SuperCDMS's efficiency is taken to vanish below the 1.6 keV threshold, hence the rate is zero below that energy. LUX's efficiency instead goes to zero smoothly at small values of S1. The high energy behavior of the spectrum is dictated by the form factor, which depends solely on the target nucleus and its recoil energy, for m X = 1 TeV and 100 GeV, and in fact the exponential tails of the differential rate are similar in these two cases. The velocity integral, which is almost constant within the LUX and SuperCDMS signal region for these DM masses, moves to lower energies when decreasing the DM mass, and for m X as small as 5 or 10 GeV it provides the main source of suppression of the rate at high energies (see figure 3). For this reason the high energy tail of the differential rate moves to lower energies when decreasing the DM mass, while the low energy tail is fixed by the efficiency function. For light enough DM, the velocity integral vanishes at energies below the experimental threshold and therefore there is no detectable signal, unless the energy threshold is lowered. This is what happens with a 5 GeV DM particle at LUX, as can be seen in figure 3, for which the measured rate is zero.
The lighter mediators move the recoil spectrum to lower energies as expected from eq. (2.2). A more quantitative understanding of the transition between the long-range regime and the contact regime can be obtained by looking at table 1. For a 1 TeV DM particle as in our benchmark model A, the typical momentum transfer is 127 MeV for LUX and 74 MeV for SuperCDMS. Therefore, the DM-nucleus scattering for benchmark model A, with a 3 MeV mediator, is deeply in the long-range regime m φ q for both experiments and the corresponding recoil spectrum is clearly different from the spectrum of the SI model with contact interaction. The same can be concluded for the model with a 9 MeV mediator. Scattering for the benchmark model B, with m X = 100 GeV and m φ = 15 MeV, again occurs in the long-range regime, while for a 45 MeV mediator we have m φ ≈ q and therefore the scattering occurs in between the pure long-range and the contact regimes. The m φ ≈ q regime is also relevant for scattering in SuperCDMS for the benchmark model C, with a 10 GeV DM particle and m φ = 20 MeV, while a 60 MeV mediator is much closer to the contact interaction model. In LUX this model is difficult to distinguish from WIMP models with contact interaction, because the scattering occurs close to the energy threshold of the detector. In fact, because of the small DM mass, the high energy tail of the spectrum (dominated by the exponential fall-off of the velocity integral η 0 at large v min ) gets very close in energy to the low energy tail (dominated by the efficiency function ). Neither η 0 nor depend on the mediator mass, and the only difference between the various models can be seen at the peak of the recoil spectrum, where the (q 2 + m 2 φ ) −2 dependence of the cross section has a tiny effect. The same behavior can be seen for benchmark model D (m X = 5 GeV) in SuperCDMS.
We have argued that long-range and contact interactions can be distinguished by the different spectrum at direct detection experiments, since the recoil spectrum is flatter at low energies for contact interactions. However, since only energies above a certain threshold are accessible at the experiments, a light WIMP could fake the steeper SIDM spectrum. In fact, as discussed above, the flat part of the spectrum for a light WIMP lies below threshold for both LUX and SuperCDMS, so that these experiments can only observe the tail due to the exponential drop of the velocity integral. This is the reason why the assumption of a heavy mediator favors lighter DM in the analysis of direct detection data, while assuming a longrange interaction points toward heavier DM, as noticed in ref. [23]. From a visual inspection of the recoil spectra we ascertain that 10-20 GeV WIMPs can mimic the SIDM spectra, depending on the target material but almost independently of the DM and mediator mass in the SIDM model. To assess the ability of the experiments to distinguish the two signals, we show in figure 6 the spectrum of a 20 GeV DM particle with contact interaction and compare it with the spectrum from our benchmark model B (m X = 100 GeV, m φ = 15 MeV), for both LUX (left) and SuperCDMS (right). It is evident from the upper panels that the WIMP and SIDM spectra at LUX are very similar. However, the two spectra could be distinguished by an experiment employing a different target, such as SuperCDMS. Furthermore, the two signals have very different modulated spectra, as shown in the bottom panels of figure 6. The modulated spectrum of the WIMP model is almost entirely positive above threshold, as is typical for light DM particles, while the SIDM spectrum is partially negative because of the large DM mass. The shape of the modulated spectrum is explained in detail in the next section for a sodium iodide detector like DAMA. Since the modulated rate is about 30-60 times smaller than the unmodulated rate, we expect that hundreds of events will be needed in order to distinguish SIDM from WIMPs using modulation signals.
Our results for this section can be summarized as follows. Despite the limitations imposed by detector resolution and efficiency, direct detection experiments such as LUX and SuperCDMS can potentially distinguish SIDM and WIMP recoil spectra, provided the DM mass is heavy enough that the scattering does not occur too close to the experimental threshold. While light WIMPs may fake a SIDM signal at a direct detection experiment, the degeneracy between the two spectra could be lifted either by a second experiment employing a different target, or by observing the modulated part of the spectrum. While LUX has a larger sensitivity to heavy DM and a much higher exposure, SuperCDMS is more sensitive to low-mass DM particles because of its lighter target and low energy threshold.

Modulation signal
Together with the average rate at LUX and SuperCDMS, we also study the prospects to distinguish SIDM from WIMPs using the annually modulated part of the rate. While the modulated rate is usually much smaller than the average rate and therefore more difficult to detect, it is a cleaner signature of DM detection because the modulation phase has to be consistent with Earth's motion through the DM distribution.
The DAMA Collaboration has claimed a 9.3σ evidence for modulation in their data [51], although this result is in strong tension with experiments measuring no signal. Here we take DAMA as an example to study the annual modulation signal. Our results will apply straightforwardly to similar future experiments such as KIMS-NaI, ANAIS, DM-Ice17, and SABRE (see e.g. [52] and references therein). Our analysis follows that in the previous section, with the modulated rate given by eq. (3.1) with η 0 substituted by η 1 . Analogous to figure 3, we plot in figure 7 the functions η 1 (v min ) in v min -space (top) and η 1 (v min (E R , m X )) in E R -space (bottom), for the two target nuclei used in the experiment, sodium (left) and iodine (right). The solid lines are for our fiducial SHM parameters, while the dashed lines are for the larger values reported in the previous section. The horizontal black line indicates the zero of η 1 ; positive values indicate that the modulation phase t 0 in eq. (2.11) corresponds to the time of maximum of the signal at v min 200 km/s, while negative values indicate that t 0 becomes the time of minimum rate for v min 200 km/s. The experimental range in detected energy, E ∈ [2 keV, 20 keV], is mapped onto a different range in E R for each target because of the different quenching factors for sodium and iodine, Q Na = 0.3 and Q I = 0.09, respectively (the quenching factor being the average fraction of recoil energy that is recorded as detected energy). The corresponding v min ranges are displayed as colored bands for each choice of DM mass in our four benchmark models in the top panels, and the E R range is delimited by the two vertical dashed lines in the bottom panels. For m X = 10 GeV and m X = 5 GeV, the signal window continues to the right of the plotted range for sodium, and lies completely outside of the plotted range for iodine. From the figure it is apparent that scattering off sodium is marginally sensitive to a 5 GeV DM, while it is fully sensitive to m X = 10, 100, and 1000 GeV. On the other hand, because of its larger mass, scattering off iodine is completely below threshold for light DM, m X = 5 and 10 GeV.
As mentioned in section 2, our form of the modulation signal relies on the DM distribution being isotropic in its rest frame. DM substructure, as e.g. dark disks and streams, will in general modify the modulation amplitude and its phase. We are also neglecting the anisotropy induced by the Sun's gravitational potential, which modifies the phase of the modulation below ∼ 200 km/s [40,41].
In figure 8 we plot the modulated spectrum for DAMA, in analogy with the unmodulated spectrum for LUX and SuperCDMS in figure 5. The modulated spectrum is shown for each of our four benchmark models (solid red lines), together with the additional model with three times the mediator mass (dashed purple lines) and the WIMP model with SI contact interaction (dotted green lines). The three spectra are normalized such that the total modulated rate in the signal window of the experiment is the same (same number of events after subtracting the unmodulated rate). The plots show that it may be possible to distinguish SIDM from WIMPs with the modulated component of the rate, especially if the experimental sensitivity extends to low energies.

Conclusions
DM self-interactions are a promising avenue to solve the small-scale puzzles in galaxy formation. If the self-interations are modeled as a Yukawa potential, the mediator mass required to generate large enough self-interactions to affect the internal structure of galaxies should be generically below 100 MeV. In order to avoid overclosing the universe, the mediator must be unstable. We have assumed that the mediator decays to lighter SM particles, which opens up the possibility of searching for SIDM particles in direct detection experiments.
In this paper we have studied SIDM direct detection in detail. While we have focused on SIDM benchmarks, our work is applicable to all DM models with long-range interactions with nuclei. We have shown that DM direct detection experiments are remarkably sensitive to the SIDM parameter space, even if the coupling between the two sectors is extremely feeble. For example, LUX has excluded all the favored SIDM region with DM heavier than 7 GeV if the kinetic mixing parameter is larger than 10 −9 . Models with DM particle masses down to 3 GeV can be excluded by SuperCDMS for values of the kinetic mixing parameter an order of magnitude larger.  In contrast to the point-like interaction in usual WIMP-nucleus scattering, SIDM interacts with nuclei through a long-range force, which leads to a measurably different signal spectrum. When the mediator mass is comparable to the momentum transfer of nuclear recoils, the SIDM spectrum is peaked more towards low recoil energy compared to usual WIMPs, and can be potentially tested by direct detection experiments such as LUX, Super-CDMS, and sodium iodide experiments like DAMA. Our analysis shows that we may be able to distinguish long-range (SIDM) and contact interaction spectra with detections in experiments with different targets or by measuring both the time-average and the modulated rate.