The Double-Dark Portal

In most models of the dark sector, dark matter is charged under some new symmetry to make it stable. We explore the possibility that not just dark matter, but also the force carrier connecting it to the visible sector is charged under this symmetry. This dark mediator then acts as a Double-Dark Portal. We realize this setup in the \emph{dark mediator Dark matter} model (dmDM), featuring a fermionic DM candidate $\chi$ with Yukawa couplings to light scalars $\phi_i$. The scalars couple to SM quarks via the operator $\bar q q \phi_i^* \phi_j/\Lambda_{ij}$. This can lead to large direct detection signals via the $2\rightarrow3$ process $\chi N \rightarrow \chi N \phi$ if one of the scalars has mass $ \lesssim 10$ keV. For dark matter Yukawa couplings $y_\chi \sim 10^{-3} - 10^{-2}$, dmDM features a thermal relic dark matter candidate while also implementing the SIDM scenario for ameliorating inconsistencies between dwarf galaxy simulations and observations. We undertake the first systematic survey of constraints on light scalars coupled to the SM via the above operator. The strongest constraints are derived from a detailed examination of the light mediator's effects on stellar astrophysics. LHC experiments and cosmological considerations also yield important bounds. Observations of neutron star cooling exclude the minimal model with one dark mediator, but a scenario with two dark mediators remains viable and can give strong direct detection signals. We explore the direct detection consequences of this scenario and find that a heavy $\mathcal{O}(100)$ GeV dmDM candidate fakes different $\mathcal{O}(10)$ GeV WIMPs at different experiments. Large regions of dmDM parameter space are accessible above the irreducible neutrino background.


Introduction
The existence of dark matter (DM) is firmly established by a myriad of astrophysical and cosmological observations [1]. Nevertheless, the exact characteristics of dark matter particles remain almost completely mysterious. Weakly Interacting Massive Particles (WIMPs) are the most popular DM candidate since they arise in supersymmetry and can naturally occur with the correct relic abundance [2], but many other scenarios are possible.
Direct detection via DM-nucleus scattering [3] has made tremendous strides, with experiments like LUX [4], Super-CDMS [5] and XENON100 [6] achieving sensitivities to WIMPnucleon scattering cross sections of σ SI n ∼ 10 −45 cm 2 for a O(100 GeV) WIMP. There have also been several anomalies in the O(10 GeV) mass range [5,[7][8][9] that seem to conflict with each other, as well as with various exclusion bounds by the above experiments when assuming a WIMP-like scattering. It remains possible that some or all of these hints will be explained by something other than dark matter, especially given how challenging these measurements and their background suppression is in that mass range. Even so, past and current anomalies naturally stimulate a great deal of work by the theory community in an attempt to reconcile conflicting experimental results. The myriad of plausible models demonstrates the necessity to explore as many different dark matter scenarios as possible, lest a crucial signal be overlooked.
In most models of the dark sector, dark matter is charged under some new symmetry to make it stable. However, in light of the complex structure of the Standard Model (SM) there is no particularly strong reason to assume the dark sector to be so simple. We explore the possibility that not just dark matter, but also the force carrier connecting it to the visible sector is charged under this symmetry. This dark mediator then acts as a Double-Dark Portal.
In [10] we introduced a model to realize this scenario: Dark Mediator Dark Matter (dmDM). It features a fermionic dark matter candidate χ with Yukawa couplings to one or more light scalars φ i . These scalars carry dark charge and can only couple to the SM in pairs, realized as a nonrenormalizable coupling to quarks,qqφφ/Λ. For sufficiently light φ this can lead to direct detection via a 2 → 3 nuclear scattering process, shown in Fig. 1.
Bounds from direct detection experiments are usually analyzed assuming a contact operator interactionχχqq/Λ 2 . The shape of the resulting nuclear recoil spectrum is entirely determined by the nuclear form factor and dark matter velocity distribution. Many past models feature different nuclear recoil spectra. Examples include the introduction of a mass splitting [11][12][13]; considering matrix elements |M| 2 with additional velocity-or momentum transfer suppressions (for a complete list see e.g. [14]), especially at low DM masses close to a GeV [15]; light scalar or 'dark photon' mediators (see e.g. [13,16]) which give large enhancements at low nuclear recoil; various forms of composite dark matter [17][18][19][20][21] which may introduce additional form factors; and DM-nucleus scattering with intermediate bound states [22] which enhances scattering in a narrow range of DM velocities. Notably missing from this list are alternative process topologies for DM-nucleus scattering. This omission is remedied by the dmDM scenario, which generates a functionally unique recoil suppression and overall cross section dependence on DM and nucleus mass. Direct detection constraints on dmDM are explored in this paper in detail, and we show that a ∼ 100 GeV dmDM candidate fakes different O(10 GeV) standard WIMPs at different experiments.
Dark Mediator Dark Matter has important consequences outside of direct detection. Coupling dark matter to a light scalar can ameliorate inconsistencies between simulations and observations of dwarf galaxies [23][24][25] while being compatible with a thermal relic. Perhaps more drastic however is the unique pair-wise coupling of light scalars to SM quarks.
We conduct the first systematic survey to constrain operators of the formqqφ i φ * j /Λ ij where φ i is a very light scalar, checking a large variety of cosmological, astrophysical and collider bounds. The heaviest stable dark mediator has to be lighter than ∼ eV to avoid overclosing the universe. This makes emission during direct detection plausible. The most stringent bounds on its coupling come from observations of neutron star cooling, which require Λ 10 8 TeV for a single dark mediator. However, all constraints are easily circumvented in a model with two mediators, which can generate a strong direct detection signal. The constraints we derive are important outside of the dmDM context as well, applying to any light scalars with the above coupling to the SM.
The pairwise dark mediator coupling to quarks is not gauge invariant above the electroweak breaking scale, necessitating a UV completion. We present one such possibility featuring dark vector quarks, leading to discoverable TeV scale LHC signatures. This paper is organized as follows. In Section 2 we define the dark mediator Dark Matter model and outline how dmDM could be realized in a UV-complete theory with its own set of LHC signatures. Section 3 summarizes bounds on the dark matter Yukawa coupling to dark mediators. In Section 4 we derive stellar astrophysics bounds on dark mediators coupled to SM quarks, which give the most powerful constraints on our scenario. Cosmology and LHC experiments also yield important bounds, which are discussed in Section 5. A realistic model of dmDM, which avoids all constraints, is defined in Section 6. Section 7 reviews the direct detection phenomenology of dmDM, and we conclude in Section 8. Some technical details and additional calculations are presented in the Appendices.

Dark Mediator Dark Matter
In this section we define the Dark Mediator Dark Matter model and discuss a possible UVcompletion involving heavy vector-like quarks that could be discoverable at the LHC.

Model Definition
Given its apparently long lifetime, most models of DM include some symmetry under which the DM candidate is charged to make it stable. An interesting possibility is that not only the DM candidate, but also the mediator connecting it to the visible sector is charged under this dark symmetry. Such a 'dark mediator' φ could only couple to the SM fields in pairs, at leading order. There are several possibilities for writing down a dark-mediator model. However, if the mediator couples via additional derivatives or through loops, direct detection is suppressed below observable levels. This limits the choice of dark mediator couplings to the simple construction introduced in [10], which we repeat here.
Consider real or complex SM singlet scalars φ i coupled to quarks, along with Yukawa couplings to a Dirac fermion DM χ. The relevant terms in the effective Lagrangian are where . . . stands for φ, χ mass terms, as well as the rest of the dark sector, which may be more complicated than this minimal setup. This interaction structure can be enforced by a Z 4 symmetry. The first two terms dictate the dark sector's interaction with the SM, while the quartics are only important in the early universe (see Section 5). 1 The leading order process for DM-nucleus scattering is χN →χN φ if m φ O(10 keV). However, an elastic scattering χN → χN is always present at loop-level since it satisfies all possible symmetries, see Fig. 1. This low-energy 2 → 2 loop process is equivalent to the φ φ u χ v ψ Q 1 ψ Q 2 ψ u 1 ψ u 2 Q χ c Figure 2. The 2 → 3 direct detection scattering process within the UV completion of dmDM. When treating Higgs vev as a mass insertion, the propagator of heavy Dirac quark is dominated by the chirality-flipping piece, , at low energy. This gives the suppression scale in eq. (2.7). operator y 2 (for n φ = 1) in the massless φ limit, where q = √ 2 m N E r is the momentum transfer in the scattering. 2 Effectively, this is identical to a standard WIMP with aχχN N contact operator, but with an additional 1/E r suppression in the cross section. This gives a similar phenomenology as a light mediator being exchanged at tree-level with derivative coupling.
The main new features of this model for direct detection in Section 7 are captured by the minimal case with a single mediator n φ = 1. However, the actual number of dark mediators is important for interpreting indirect constraints in Sections 3, 4 and 5. It also affects the relative importance of the two nuclear scattering processes. When n φ = 1, the 2 → 3 process will dominate direct detection for Yukawa coupling y χ below some threshold as long as m φ keV. If n φ = 2, however, the dominant scalar-DM coupling could beqqφ 1 φ * 2 /Λ 12 . In that case, the 2 → 2 operator above is ∝ y φ 1 χ y φ 2 χ and can be suppressed without reducing the 2 → 3 rate by taking y φ 1 χ y φ 2 χ . Both processes will be considered for direct detection in Section 7. The effect of strong differences between proton and neutron coupling to DM have been explored by [26]. To concentrate on the kinematics we shall therefore assume the operator qqφφ * /Λ is flavor-blind in the quark mass basis.
We point out that depending on the UV completion of the model, a leptonic coupling via¯ φφ * is also possible. We do not consider it here, since direct detection would be very difficult, but indirect constraints, in particular from white dwarf cooling, could be sensitive to such a scenario.

A possible UV-completion
Above the electroweak symmetry breaking scale theqqφφ * /Λ operator is realized asQ L Hq R φφ * /M 2 . This is suggestive of a particular UV completion involving heavy vector-like fermions coupling to φ and SM quarks via Yukawa couplings. The minimal particle content to realize dmDM is therefore a light scalar mediator φ, heavy vector-like quarks ψ Q, q in the same gauge representations as the SM Q L , u R , d R respectively, and a Dirac fermion dark matter candidate χ. Their charges are shown in Table 1. The Lagrangian 3 contains Yukawa couplings where the index 1, 2 represents the chirality component of Dirac fermion ψ's. ψ Q,u,d have Dirac masses The DM mass and its coupling to φ are given by Note that in this limit, the process has an IR pole similar to tree-level t-channel exchange, hence the q −1 dependence. 3 We show the n φ = 1 complex scalar case, generalization to more or real dark mediators are trivial.
We assume all the couplings are flavor universal and M Q = M q , y Q = y q for simplicity. The direct detection 2 → 3 scattering process is shown in Fig. 2. When the momentum transfer through heavy quarks is much smaller than M Q , we can integrate out the lower part of the diagram to generate the dimension 6 operator Below the scale of electroweak symmetry breaking this becomes the operator of eq. (2.1) with where v = 246 GeV is the SM Higgs VEV. As we will show, M Q could easily be TeV scale, allowing for discovery of these heavy vector-like quarks at the LHC. As long as the LHC with √ s = 8 TeV has not produced them on shell they are not trivially excluded despite being new colored states that couple to the Higgs. Since they do not receive their mass primarily from the Higgs vev, their contribution to the hγγ loop coupling is strongly suppressed. As we discuss in Section 5, the collider constraints on additional vector-like quark generations can be satisfied for M Q TeV. The quark Yukawa couplings do receive a flavor-universal correction which may lead to the light quark Yukawa couplings being tuned to the order of 0.1%, but like the origin of the light scalar φ we put these naturalness issues aside to concentrate on the phenomenology of dmDM.

Constraining the DM Yukawa Coupling
The dark matter Yukawa coupling y χ χ c χφ can be constrained by various astrophysical and cosmological observations, the most important of which we summarize here. For simplicity these bounds are formulated for n φ = 1, but can also be applied directly to n φ > 1 scenarios if one Yukawa coupling dominates.
The dark matter relic density Ω CDM = 0.1196 ± 0.0031 has been accurately measured by the Planck Satellite [1]. Under the assumptions of a simple thermal relic this fixes y χ to a specific value (which depends on m χ ). The lowest-order annihilation cross section for the process χχ → φφ * is assuming no sizable φ 3 couplings. Performing the standard WIMP freeze-out calculation [27] we find that the φφ * ↔χχ process freezes out at the usual T ∼ m χ /20. Requiring that Ω χ = Ω CDM gives  Figure 3. Bounds on the Yukawa coupling y χ χ c χφ for n φ = 1. (These bounds can also be applied directly to n φ > 1 scenarios if one Yukawa coupling dominates) Magenta: required value of y χ for χ to be a thermal relic. Cyan and Orange: upper bounds on y χ from bullet cluster and ellipticity observations. The green shaded region implements the SIDM solution to the core-cusp and too-bigto-fail problems of dwarf galaxies [23][24][25], while the pink region can modify the halo of milky way size galaxies. See text for details. Black curve: 2 → 3 dominated direct detection requires y χ to lie below this curve if n φ = 1, see Section 7.3. This is generically very small, of order 0.01 for ∼ 10 GeV DM, and is compared to the other y χ bounds in Fig. 3 (magenta line). We emphasize that this constraint will be shifted if χ is nonthermally produced. Although DM interaction is mediated by light scalars, the Sommerfeld enhancement, which is proportional to [91] is negligible due to the small Yukawa coupling y 2 χ , as well as the relatively large velocity v 0.3 during freeze-out.
An upper bound on the dark matter self-interaction may be obtained from observations of the Bullet Cluster and galactic ellipticities. This was done by the authors of [28] for a massless mediator. We can apply those bounds directly to our model as long as m φ is much smaller than the momentum transfer of a characteristic DM-DM collision (q MeV for m χ GeV). The bullet cluster bound is considered quite reliable, but concerns have been raised about the ellipticity bound, the strength of which may have been overestimated [29]. Both upper bounds are shown in Fig. 3 (cyan and orange lines). Rather than merely requiring the light mediator to not spoil well-understood aspects of galaxy formation and interaction, one could go one step further and use the dark matter selfinteraction to address existing inconsistencies between prediction and observation. Current N -body simulations of Cold Dark Matter halos predict an overabundance of dwarf spheroidals, as well as dwarf galaxy halos that are more cusped than observed. These inconsistencies are called the too-big-to-fail and core-cusp problems. It has been shown that the disagreement between simulations and observation can be ameliorated by introducing a sizable dark matter self-interaction, dubbed the Selft Interacting Dark Matter (SIDM) scenario [23][24][25].
The presence of a light scalar in the m φ MeV mass range allows dmDM to act as a realization of SIDM. To derive the preferred range of y χ we follow the procedure in [25].
The small ratio between the potential energy of φ mediation and the kinetic energy of DM in galactic halos, 2α χ m φ /(m χ v 2 ) 1, shows that DM self-interaction should be described in the classical limit. The transfer cross section for DM scattering, is just the total cross section weighted by fractional longitudinal momentum transfer. A value of could reconcile the inconsistencies between N -body simulations and observations. The required coupling depends on the ambient dark matter velocity, which is ∼ 30 km/s for dwarf galaxies and ∼ 300 km/s in larger milky way size galaxies. Fig. 3 shows the preferred bands of y χ to achieve the cross section eq. (3.6) in these two systems. In this plot, m φ = MeV, but the change for m φ = eV is not substantial. 4 As we can see, the dmDM model with a thermal relic DM does provide a potential solution to the core-cups and too-big-to-fail problem of dwarf galaxies. Finally, as we discuss in Section 7.3, there is an upper bound on y χ for the 2 → 3 process to dominate direct detection when n φ = 1. If y χ is larger, direct detection proceeds via the 2 → 2 loop process. This is shown for the LUX experiment as the black line in Fig. 3. (The corresponding upper bound for other experiments is somewhat weaker.) Note that this boundary between the two direct detection regimes is arbitrarily shifted for n φ = 2.
In summary, Fig. 3 shows both the preferred values of y χ for a thermal relic and to resolve inconsistencies between observations and simulations for dwarf galaxies and the milky way; it also shows the upper bounds on y χ to satisfy bullet cluster and self-interaction bounds, and to ensure 2 → 3 dominated direct detection. Roughly speaking, the most relevant values of y χ are ∼ 10 −3 − 10 −2 . is produced in the early universe, as well as stellar cores and high energy colliders. In this section we compute m φ -dependent bounds on Λ from stellar astrophysics. The light scalar φ is produced in stellar cores if m φ T . This can affect the length of the neutrino burst in supernovae explosions, radiative heat transfer and energy loss in the sun, and the cooling of stellar relics. We assume n φ = 1, but the constraints are easily applied to the more general case.

Constraining the Dark Mediator φ through Stellar Astrophysics
The derivation of these bounds differs from the corresponding calculations for axions, since light scalars couple more strongly at low energy due to the scaling of the operator eq. (4.1). In the regime where respective bound can be set, φ fully thermalizes in the sun and white dwarfs. By far the strongest constraints are obtained from observations of neutron star cooling: Λ 10 8 TeV for m φ 100 keV, which excludes this scenario for direct detection completely. However, in Section 6 we construct n φ = 2 scenario with one eV and one MeV dark mediator that evades all constraints while allowing for sizable direct detection signals.
It is useful to keep in mind the range of Λ relevant for direct detection. As discussed in Section 3, the preferred range for the dominant DM Yukawa coupling is y χ ∼ 10 −3 − 10 −2 . Direct detection bounds on dmDM were computed in [10] and are reviewed in Section 7. For dmDM to be detectable above the irreducible neutrino background, Λ 10 4 TeV in the relevant dark mediator coupling to quarks.

φ interaction and production cross sections
Computing stellar astrophysics and cosmological bounds requires an understanding of the φ-nucleus scattering cross sections at sub-GeV energies. This is easily computed analytically using standard methods for DM scattering and is shown in Fig. 4 for Λ = 10 TeV. For illustration we also compare these cross sections to some relevant SM scattering processes, νN → νN and Compton scattering. Note the different energy scaling of these cross sections, with σ(φN → φN ) being independent of energy for E φ 100 MeV.
At tree-level, φ only couples hadronically. Therefore, the most relevant production processes for φ in stellar cores are Again we are only concerned with sub-GeV energy scales. We can model the first two processes, shown in Fig. 4, in MadGraph5 by treating the proton as a fundamental fermion and multiplying the cross section by a quark-nucleon matrix element factor, see eq. (7.3). The Helm form factor eq. (7.4) is also included for nuclei. The one-pion exchange approximation was employed for the first process [31], and the obtained cross section should be seen as an O(1) estimate. The first process can occur off any nucleus, with N = p, He shown in Fig. 4 (the cross sections for N = He, C, O are nearly identical), which is relevant in the Sun and white dwarfs. The second process proceeds identically for protons and neutrons and is relevant in neutron stars, with additional subtleties due to neutron degeneracy discussed in Section 4.5.  Figure 4. Low-energy scattering and production cross sections for computing bounds on the dmDM model, compared to some relevant SM processes. Solid lines: coherent scattering of φ off a stationary nucleus via the operatorqqφφ * /Λ. Dashed lines: coherent scattering of a neutrino off a stationary nucleus via Z-exchange (see also [30]). Long-Dash-Dotted lines: Compton scattering of a photon off a stationary electron or proton. Long-Dashed lines: pp → ppφφ * , γp → pφφ * and γHe → He φφ * where one initial proton is stationary. The blue band represents a naive dimensional analysis estimate eq. (4.4) of γγ → φφ * (or the reverse annihilation process φφ * → γγ), taking B 2 = 1 − 100. In all cross sections involving theqqφφ * /Λ operator we used Λ = 10 TeV.
The photon annihilation process γγ → φφ is difficult to calculate due to unknown form factors connecting quarks to hadronic QCD resonances. A rough estimate of the amplitude can be obtained by treating it as a loop process mediated by constituent quarks. The same approach is used to calculate the photon meson couplings, for example in [32]. With the correct power of electric charge and one mass insertion for the correct chirality, the size of operator |φ| 2 F µν F µν is approximated as where B is the form factor between the φ and constituent quarks, and m u,d 263 MeV [32] is the mass of the constituent quarks within the NJL model. The resulting cross section is The blue band in Fig. 4 is a very rough estimate with B 2 = 1 − 100. At our level of precision we also take this to be the cross section for the reverse annihilation process φφ → γγ.

Supernovae
Like massless axions, production and emission of φ's can lead to rapid energy loss during a supernova explosion. This can be constrained by measuring the duration of the associated neutrino burst. There are two allowed regimes [39]. The φ are trapped in the stellar medium if they couple more strongly to the SM than neutrinos. In that case they do not affect the neutrino burst. Alternatively, if the SM copuling is 5 orders of magnitude weaker, φ production is too negligible to affect the supernova. Rescaling σ φN →φN ∝ Λ −2 at E φ ∼ 10 MeV from Fig. 4, we see that the former constraint is satisfied for Λ 10 6 TeV. Therefore, supernova roughly supply the bound Λ 10 11 TeV or Λ 10 6 TeV (4.5) on Λ involving scalars with a mass of m φ 10 MeV, the temperature of a supernova explosion.

Solar Energy Loss and Radiative Heat Transfer
A stellar core at some temperature T GeV can be seen as a fixed target experiment in which slow-moving nuclei are bombarded by photons as well as relativistic electrons and, in the case of dmDM, φ scalars. The most relevant production processes for φ are shown in eq. (4.2), with cross sections as a function of energy illustrated in Fig. 4. The φ production rate per second per volume via a process X 1 X 2 → φφ * + SM particles, with cross section σ φprod and parent particle number densities n X i , is (4.6) (We assume φ is so light that it is always relativistic.) On the other hand, the mean free path for φ before it scatters off nuclei in the star is for N i = {p, He} in the case of the Sun, with additional heavier elements in white dwarfs. Estimating energy loss due to φ emission from the star is greatly simplified if we can make four assumptions: (a) The effect of φ production is small enough so as to not significantly influence the evolution of the star, allowing us to treat it as a background source of φ's. 5 (b) φ particles are produced predominantly at the center of the star.
(c) L φ is short enough that φ scatters many times and thermalizes before leaving the star.
(d) There is negligible φ annihilation in the star.
If these conditions are satisfied, the created φ particles diffuse outwards from the center until they reach a layer of low enough density so that the surface of the star is within ∼ one scattering length, at which point they escape. Each φ carries away energy is the temperature of the 'layer of last scattering'. 6 In the absence of annihilation processes, the equilibrium rate for φ emission is equal to the total rate of φ production, which together with T escape φ gives the total energy loss from φ emission. We now perform this computation for the case of the Sun. Radial density, temperature and mass fraction profiles for the standard solar model can be found in basic astrophysics textbooks and are reproduced in Appendix A for reference. The radius of the sun is about R sun ≈ 3.85 × 10 26 cm, while the central density and temperature are ρ sun (0) ≈ 150 g cm −3 and T sun (0) ≈ 1.5 × 10 7 K ≈ 1.3 keV. The corresponding nucleus number densities are of order 10 25 cm −3 , while the density of photons obeying a Bose-Einstein distribution is n γ = 2ζ(3) Fig. 4 and eq. (4.6) it is clear that γN → N φφ * is the dominant production process for T ∼ keV.
The φ creation rate per volume as a function of distance R from the sun's center is where n γ andσ are evaluated at temperature T sun (R), andσ is the thermally averaged cross section for a Maxwell-Bolzmann distribution (excellent approximation of Bose-Einstein in the sun) of photons hitting a stationary nucleus:  The resulting r create φ (R) ∝ Λ −2 is shown in Fig. 5. About 90% of φ production takes place within 0.2 solar radii, validating assumption (b) above. The total rate of φ creation in the entire sun is (4.10) For our purposes here, define R core = 0.2R star . Since most of the φ creation takes place within that radius, is satisfied up to a factor of two. We next compute the φ mean free path L φ (R) ∝ Λ 2 via eq. (4.7) using similarly averaged scattering cross sections. This is shown in Fig. 6 for Λ = 10 TeV. For this benchmark value assumption (c) is certainly satisfied. The 'layer of last scattering' is situated at R ≈ R escape where This allows us to define the approximate temperature of the escaping φ's as Both R escape φ and T escape φ are shown as functions of Λ in Fig. 7. Assumption (c) holds for Λ 100 TeV. On the other hand, our calculations become unreliable around Λ ∼ 1 TeV since we then become sensitive to details of the sun's surface structure.
We can now estimate the fraction of the star's luminosity in the form of φ emission, making use of the (yet to be verified) assumption (d), which gives at equilibrium: (4.14) Therefore, the power of φ emission is The sun's measured power output is P sun ≈ 3.85×10 26 Watts. The ratio P φ /P sun as a function of Λ is shown in Fig. 8. The φ contribution becomes negligible 7 for However, as we will see below, this does not constitute the strongest bound obtained from the sun. 7 Compare to neutrino emission Pν /Psun ≈ 2%. [40] We still need to verify that assumption (d) holds. Evaluating the rate of φ annihilation in the sun requires us to solve for the equilibrium number density n φ (R). We can construct the associated diffusion equation with the information assembled here, but numerically solving it is beyond the scope of this work. However, we can make a ball-park estimate of the total equilibrium φ population by noting that the time taken for a single φ to escape is dominated by the time taken to diffuse from the dense core: Fig. 6). This means N φ , the equilibrium total number of φ's in the sun, is approximately given by solving Assuming all these φ's live in the core, the corresponding number denisty is Consulting Fig. 4 and comparing with number densities n γ ∼ 10 22 cm −3 and n N ∼ 10 25 cm −3 in the core, it is clear that the φ annihilation rate is completely negligible compared to the creation rate in eq. (4.8).
To make sure the sun is not disturbed by φ production we also have to ensure that radiative heat transfer, which dominates the core and radiative zone, is relatively unaffected. The radiative heat transfer due to φ should be compared to the photon heat flux [41]: Substituting eq. (4.20), n γ (T sun (0)), L φ (0) as well as L γ ∼ 10 −2 cm (see Fig. 20 in Appendix A), we obtain the following heat transfer ratio in the core: which is significantly stronger than the bound from energy loss due to φ emission.
Two final remarks are in order. Firstly, as mentioned above, Fig. 7 shows that assumption (c) of trapped and thermalized φ's starts breaking down when Λ 100 TeV. In that case φ no longer contributes to radiative heat transfer, while the lost power due to φ emission is roughly P φ ∼ T core R create φ ∼ (10 −5 P sun )(100 TeV/Λ) 2 . The free-streaming regime in the sun therefore sets no constraints on Λ. Secondly, we also point out that there is a sub-population of protons and photons with E ∼ 10 MeV produced by fusion reactions in the sun, but the total rate of fusion reactions R fusion ≈ 3.6 × 10 38 s −1 is many orders of magnitude too low for this subpopulation to affect our estimates.

White Dwarf Cooling
White dwarfs (WD) represent the evolutionary endpoint of stars up to several solar masses. They are supported by electron degeneracy pressure, which largely decouples their hydrostatic and thermal properties and results in a strong relationship between their mass and radius. Since white dwarfs do not support fusion processes in their cores they simply cool down after they are formed, with observable luminosities ranging from 0.5 to ∼ 10 −4 L sun , corresponding to core temperatures of around 10 to 0.1 keV (about 10 8 and 10 6 K) [40].
Their relative simplicity makes white dwarfs suitable for constraining new physics with light particles (see e.g. [40,42]). Unlike the Sun, where we have a single well-studied star to compare predictions to, white dwarf cooling is constrained by the White Dwarf Luminosity Function (WDLF), which is the number of observed WDs at different luminosities, see Fig. 10.
For reasonable assumptions about the star formation rate, the shape of this WDLF curve is given entirely by the WD cooling rate [40].
The large central density of white dwarfs ρ W D ∼ 10 6 g cm −3 means φ's can be copiously produced, but also thermalize completely before diffusing out of the star. This makes their emission somewhat sensitive to details of WD stellar structure, unlike e.g. free-streaming axions. Comprehensively studying the constraints on the 1 Λ qqφφ * operator set by WD cooling would therefore require modeling a representative WD population, which is beyond the scope of this work.
Fortunately, the WD population in our stellar neighborhood is strongly peaked around 0.5 − 0.7 solar masses [40,43]. This means we can obtain a preliminary estimate of the bound on Λ by studying a single star in this representative mass range.
Our benchmark dwarf (about 0.5 solar masses) started its life as a roughly one solar mass main sequence star that was evolved forward in time using the stellar evolution code MESA [44]. 8 Most of the observational data in the WDLF is for luminosities 0.1L sun , which corresponds to a bolometric magnitude M bol > 7. Photon cooling, well-described by Mestel's Law [45], dominates for such cool white dwarfs. We therefore compute φ-cooling in our benchmark dwarf for L WD < 0.1L sun .
Since the degenerate electron gas in WD cores is an excellent conductor of heat, radiative heat transfer is unimportant. We therefore only compute the total power loss due to φ emission, in an identical manner to the previous subsection. As we will see, assumptions (a) -(d) are satisfied throughout as long as Λ is large enough. Radial profiles of density, composition and temperature produced by MESA for our benchmark dwarf are shown in Appendix A.
The φ creation rate per volume is shown in Fig. 5 for Λ = 10 TeV when the white dwarf has 0.1 solar luminosity. Due to the similar temperature but larger density, it is 5 orders of magnitude higher than in the sun. The pp → ppφφ process is still strongly temperaturesuppressed. Figs. 6 and 7 show that φ's do not escape until they are very close to the white dwarf surface. The resulting power loss is shown in Fig. 8, and φ emissivities are compared to known photon and neutrino emissivities in Fig. 9. To a reasonable approximation, (4.25) Requiring φ emission to represent a subdominant 10% fraction of the total WD luminosity requires Λ 40 TeV, but as it turns out the actual bound on Λ from the white dwarf luminosity function is significantly less constraining. We now compute this bound following the procedure in [40].
The white dwarf looses internal energy U with time due to emission of photons, neutrinos and (in our case) φ's, so that dU/dt = −(L γ +L ν +L φ ), where L γ is the total photon luminosity of the star. Assuming a constant star formation rate B, the number density of white dwarfs in a given magnitude interval is proportional to the time it takes to cool through that interval,  Figure 9. Photon (short-dashed) and neutrino (long-dashed) emissivities for typical white dwarfs with ρ core = 10 6 g cm −3 as a function of core temperature [40,45]. Solid lines indicate φ emissivities obtained for our benchmark dwarf.
For a white dwarf, the photon emissivity can be be computed using Kramer's opacity law: (4.27) Bolometric magnitude gives the photon luminosity relative to the sun, log 10 (L γ /L sun ) = (4.74 − M bol )/2.5. Therefore T ∝ 10 −4M bol /35 and we get where we have absorbed details of the star formation rate and white dwarf heat capacity into the constant C. (When comparing to observational data it is conventional to normalize this constant to the data point with the smallest uncertainty.) For pure photon cooling, this reduces to the well-known Mestel's cooling law [45], shown as the black dashed line in Fig. 10. This already gives a reasonable fit, but full simulations (thin green curve) are needed to account for the observed data in detail. For the range of core temperature we consider, neutrino cooling can be neglected. Using the φ emissivities of our simulated benchmark dwarf in eq. (4.28), we obtain the modifications to Mestel's cooling law shown in Fig. 10 for different values of Λ. Given the crudeness of our estimate a full fit to the data is not appropriate. However, we can estimate the sensitivity of a full stellar simulation to φ cooling by the size of the deviation from Mestel's law. Given the scale of astrophysical uncertainty in the WDLF (illustrated by the gray band in Fig. 10), a reasonable rough bound on the allowed modification to standard white dwarf cooling is Λ 10 TeV. (4.29) The effect of φ cooling is more pronounced for young, hot white dwarfs (smaller M bol ). A more thorough study, including full simulation of φ cooling throughout the life of the white  Figure 10. Green and magenta datapoints: measured white dwarf luminosity function from [46] and [47], showing the number density of observed white dwarfs per bolometric magnitude interval per pc 3 . The gray band indicates how much the WDLF of [46] would change by varying the assumed scale height of the galactic disk between 200 and 350pc. Green line: WDLF from a full simulation, assuming constant star formation rate, taken from [42]. Black dashed line: Mestel's cooling law (pure photon cooling). Colored solid lines: modification of Mestel's law due to additional φ cooling for Λ = 1, 5 and 10 TeV. All cooling curves have been shifted to pass through the datapoint with the smallest uncertainty.
dwarf, might therefore give a somewhat more stringent bound on Λ. However, as we see in the next section, a much stronger constraint is supplied by neutron star cooling. Finally, one might worry about φ being produced in electron collisions or plasmon decay via its loop-induced coupling to the Z-boson, see eq. (B.8). However, this coupling is ∼ 10 −3 (10 TeV/Λ) smaller than the equivalent tree-level electroweak coupling. According to the discussion in [42], φ emission from the plasmon decay and electron Bremsstrahlung is therefore < ∼ 10 −8 (10 TeV/Λ) 2 erg s −1 g −1 at T ∼ 4 × 10 7 K, which is much lower than the nuclear production discussed above.

Neutron Star Cooling
Neutron stars are the evolutionary endpoint for heavy stars that do not collapse to a black hole. They are supported by neutron degeneracy pressure and constitute the densest form of matter in the universe. This introduces many subtleties into their cooling processes, which are not yet fully understood even in the Standard Model (see e.g. [48][49][50][51][52][53][54][55][56][57]). However, neutron stars are such powerful "φ-factories" in dmDM that we can still set very strong constraints despite these uncertainties.
Neutron stars are born in hot supernovae explosions with T ∼ 10 11 K 10 MeV but quickly cool down and enter the neutrino cooling phase when their internal temperature reaches about T ∼ 10 9 K ∼ 100 keV (see e.g. [48] for a review). Neutrino cooling dominates for ∼ 10 5 years, after which photon cooling takes over. For a given equation of state, the mass of the neutron star fixes both the radius and density profile. The radius is about 10km, while the central density is ρ ∼ 2 − 10 × ρ 0 , where ρ 0 ≈ 2.8 × 10 14 g cm −3 is the density of nuclear matter at saturation.
The neutron star core extends to about 1km below the surface and is divided into an inner core (ρ 2ρ 0 ) and an outer core (ρ 2ρ 0 ). (Light neutron stars with M 1.3M sun do not have an inner core.) The characteristics of the outer core are well constrained by nuclear theory and laboratory data, while the inner core is much less well understood, with hypotheses for its composition ranging from normal nuclear matter to hyperions, pion or kaon condensates, or a pure quark fluid (called 'strange quark matter' due to the presence of s quarks). However, the recent observation of a 2 solar mass neutron star [58] is in conflict with all core hypotheses other than normal nuclear matter, which provides the only equation of state 'stiff' enough to support such large masses. We shall therefore only consider neutron stars with nucleon cores.
The neutron star is surrounded by an outer crust of thickness ∼ few 100 m, consisting of a non-degenerate neutron gas with characteristic density of order ρ N ∼ 4 × 10 11 g cm −3 . During the neutrino cooling phase the outer crust acts as a heat blanket, thermally insulating the neutron star interior against radiative losses into space. For a nonmagnetic iron envelope the surface temperature of the star can be related to the interior temperature by [59][60][61] T surface = (0.87 × 10 6 K) g s 10 14 cm/s 2 where g s = GM/R 2 star is the surface gravity. The inner crust has a thickness of ∼ 1km and forms the transition between the heat blanket and the core. The thermal conductivity of nuclear matter is so high that the interior below the blanket is close to isothermal.
Late-time cooling is constrained by ∼ 20 observations of neutron stars for which both surface temperature and age could be determined, see green data points in Fig. 12. The mass of an individual neutron star, which is not known in the dataset, determines the cooling curve T surface (M ; t). Different cooling models can be excluded by the requirement that the observed data points fall into the range of allowed cooling curves. For non-superconducting neutron stars with non-magnetic iron heat blankets in the Standard Model this range is illustrated with the two blue dashed lines in Fig. 12 [49]. Accretion of light elements in the crust and the presence of strong magnetic fields at the surface would increase the thermal conductivity of the outer neutron star layers, increasing T surface by a factor of a few during the neutrino cooling phase. Furthermore, the core may be in different phases of neutron and/or proton superfluidity, which can affect the surface temperature by at least a similar factor. Nevertheless, quite stringent constraints on φ-cooling can be obtained from light, slow-cooling neutron stars.
It is necessary to understand how the standard range of allowed cooling curves changes when φ emission is included. We will therefore estimate first the emissivity φ (T core ) and then the cooling curves T surface (t) for a light, slow-cooling neutron star and a heavy, fast-cooling neutron star, which bounds the range of allowed cooling curves. The relevant parameters of our benchmark stars are given in Table 2.
M/M sun R (km) R core (km) ρ core /ρ 0 1.1 13 11 2 2.0 11 10 10 Table 2. Parameters of light and heavy neutron stars to determine the range of allowed cooling curves in dmDM. ρ core is the central density. Adapted from [48], which assumed nucleon cores.
We model the neutron star core as a sphere of constant temperature and density. Assume for the moment that the annihilation process φφ * → γγ can be ignored, and that φ is freestreaming in both the core and the crust. In that case, the φ emissivity is given simply by (4.31) φ creation proceeds via the process nn → nnφφ * . Here we have to take into account Pauliblocking: since the neutrons in the core are strongly degenerate, only the subpopulation living on the Fermi surface can participate in reactions, and furthermore the phase space of reactions is suppressed since neutrons cannot scatter into the occupied Fermi volume. The neutron Fermi Energy E F = 2 (3π 2 n n ) 2/3 /(2m n ) is 95 MeV (280 MeV) for ρ = 2ρ 0 (10ρ 0 ). The fraction of neutrons on the Fermi surface is roughly T /E F , so we define the number density of 'available neutrons' (with kinetic energy ≈ E F ) as σ F nn→nnφφ * is much smaller than σ nn→nnφφ * from Fig. 4 due to Pauli Blocking: two neutrons with kinetic energy ∼ E F interact softly to produce two φ's with energy T so that their final energy is still on the Fermi surface. We can roughly estimate this phase space suppression ζ(E F , T ) using MadGraph, shown in Table 3. This gives where σ 0 prod = B prod × 7 × 10 −9 pb and B prod = 0.3 − 3 is a parameter we vary to account for the uncertainty of this estimate. Interestingly the cross section is constant up to a factor of ≈ 2 for E F in the range of 95 to 280 MeV, so we absorb this ρ core dependence into the uncertainty.
DefiningẼ F ≈ 95 MeV and n 0 n ≈ 3.3 × 10 38 cm −3 to be the Fermi energy and neutron number density when ρ core = 2ρ 0 , and specifying the actual number density via the dimen-  Table 3. Cross section of nn → nnφφ * for two neutrons both with kinetic energy E F , computed in MadGraph in the one-pion exchange approximation [31]. The third column gives the phase space suppression when requiring both final state neutrons to have kinetic energy in the range (E F −T, E F + T ). All quantities are understood to be ∼ estimates.  Figure 11. φ emissivity φ (T core ) for different Λ in the light neutron star defined in Table 2. The allowed range for each Λ comes from the production cross section uncertainty in eq. (4.34). Also shown is slow neutrino emission (dashed) via the modified Urca process, which occurs in all neutron star cores [62], and the effective emissivity from photon emission [59][60][61] (long-dashed). sionless ratioñ n = n n /n 0 n , we obtain This is shown in Fig. 11 for the light neutron star as a function of core temperature, and compared to the effective emissivity from neutrino and photon emission. Requiring that φcooling be subdominant to standard cooling mechanisms in the light neutron star sets the strong constraint Λ 10 8 TeV. The constraint derived from the heavy neutron star is much weaker, since the powerful direct Urca neutrino emission process is active when the central density is ρ core 2ρ 0 [63]. We have checked that for Λ 10 4 TeV, the equilibrium φ density in the neutron star is indeed small enough to render the annihilation process φφ * → γγ irrelevant. Furthermore, φ becomes free-streaming in the crust (core) when Λ 5000 TeV (500 TeV). This validates the assumptions of our estimate, and allows us to circumvent the subtleties of φ-scattering inside the neutron star core and crust (see [64] for some of the involved issues).  Figure 12. Green data points: surface temperature and age of observed neutron stars [49]. Blue dotted lines: cooling curves for heavy and light non-superconducting neutron stars with non-magnetic iron heat blankets [49]. Solid black lines: our corresponding estimate of these cooling curves using eq. (4.35) and eq. (4.36). Orange contours: estimate of the light neutron star cooling curve with φ emission for Λ = 10 6 , 10 7 , 10 8 TeV and B prod = 0.3. (For Λ > 10 6 TeV, the cooling curve for the heavy NS does not change perceptibly.) In all our estimates we multiplied T surface (t) by 0.6 (0.2 units on vertical axis) to bring them into better agreement with the full calculation by [49].
We can explicitly demonstrate the effect of φ emission on neutron star cooling. Following [65], a reasonable estimate of the cooling curve can be obtained by solving the differential equation where the specific heat for a gas of non-interacting fermions is and the Fermi momenta are p N F = (340 MeV)(2ñ n ) 1/3 and p n,e F = (60 MeV)(2ñ n ) 2/3 . The surface temperature is then approximately given by eq. (4.30). The resulting cooling curves for the heavy and light neutron star are shown in Fig. 12. (Since we are interested in the effect of introducing φ-cooling compared to the standard scenario, we multiply all our T surface (t) by 0.6 to bring our estimates into better agreement with complete cooling calculations. This corresponds to a uniform downward shift of 0.2 units on the vertical axis of Fig. 12.) The heavy neutron star cooling curve is not visibly affected for Λ 10 6 TeV. To avoid altering light neutron star cooling curves by much more than the plausible size of the effects of surface accretion, magnetic fields, and the likely presence of a superfluid component in the core [48], requires Λ 10 8 TeV. (4.38) This confirms our earlier estimate of the constraint. Light neutron stars therefore supply a very strong bound on Λ in dmDM for n φ = 1. However, as we shall see in Section 6, this constraint is easily circumvented when n φ = 2.

Other Constraints on φ
While stellar astrophysics provides the most impressive constraints on the operatorqqφφ * /Λ, cosmology and LHC searches bound the dmDM parameter space in complementary directions. We find that all cosmological constraints are satisfied as long as the only stable dark mediator is lighter than ∼ eV. LHC searches provide a constraint of Λ few TeV that does not depend strongly on n φ or m φ . Flavor physics bound could restrict the allowed SM flavor structure of the coupling in eq. (4.1), but we avoid those constraints by making the operator diagonal in the SM quark mass basis.
Dark mediators can also be probed, in principle, using fixed target experiments, precision electroweak measurements or indirect detection of dark matter annihilation. However, as discuss in Appendix B, these measurements yield no meaningful constraints.

LHC Searches
The LHC is sensitive not just to the effective coupling in eq. (4.1) but also to the UV completion of dmDM. We therefore analyze constraints in terms of the dark vector quark model of Section 2.2.
The di-jet + MET search by CMS [33] is sensitive to on-shell production of two heavy vector-like quarks via the process pp → ψ QψQ → φφ * jj. The constraint is straightforward to apply in our model, since the vector quarks are produced by gauge interactions. The resulting bound is M Q > 1.5 TeV. (5.1)

The 20/fb CMS mono-jet search [34] is sensitive to the processes
by doing a counting experiment in different missing energy bins. We simulated the dmDM signal expectation in MG5+Pythia 6.420+PGS4 [35,36] and validated our simulations by reproducing the CMS jZ(νν) background prediction with an overall scaling factor of K = 1.4. The same scaling factor was also applied to the dmDM signal. The resulting 95% CL lower bound on theqqφφ * /Λ operator depends on whether the intermediate dark vector quark is produced on-shell: Λ eff 2 (6.6) TeV for M Q = 4 (1.5) TeV.
where Λ eff = Λ for n φ = 1. For n φ > 1, since the total signal production cross section is given as the sum of all the φ i φ * i cross sections.

Cosmological Constraints
The mass of the dark mediator should be smaller than about an MeV to allow for sizable direct detection of χ. This means φ can be thermally produced in the early universe even after χ freezes out. Such a stable light degree of freedom can overclose the universe and affect Big Bang Nucleosynthesis (BBN) as well as structure formation. In this section, we discuss the thermal history of φ and derive the relevant cosmological constraints at each step.

Thermal φ production
The relic density of a light φ is given by [27] Ω where g φ = 2 is the number of degrees of freedom (d.o.f.) associated with a complex scalar and g * S is the total number of d.o.f. at the temperature at which φ decouples from the thermal bath. Given the possible size of g * S , it is clear that dark mediators with sizable couplings to the SM will overclose the universe unless the heaviest stable scalar has a mass of m φ eV. This is effectively massless for the purpose of computing all other bounds in this section, which we shall assume from now on.
Assessing the effect of φ on BBN requires knowing its freeze-out temperature more precisely. For values of Λ relevant to direct detection, the hadronic coupling to the SM bath keeps φ in thermal equilibrium at least until pions decay at T ∼ 100 MeV. After pion decay, the process φφ * ↔ γγ maintains thermal equilibrium until the time taken for two photons to annihilate exceeds the hubble time, i.e.
Substituting eq. (4.4) and the smallest possible g * ≈ 3 to slightly underestimate the freeze-out temperature, we obtain For Λ TeV, φ will decouple from the SM bath before BBN.

Big Bang Nucleosynthesis
The presence of φ during the BBN epoch (T around 10 to 0.1 MeV) can affect the generation of light elements in two ways. First, even though φ-nucleon scattering does not change the relative number of neutrons and protons, the presence of an additional light degree of freedom speeds up the expansion of the universe and makes the neutron-proton ratio freeze out at a larger value. This leads to an over production of 4 He, an effect that can be constrained by measuring the effective number of neutrino flavors, N eff during BBN. Current observation gives N eff = 3.3 ± 0.6 [1] at 95% CL from Plank+WMAP+HighL CMB observations. Since φ is relativistic it will contribute to an effective number of light neutrino flavors, to the SM value of N eff = N ν = 3. Assuming all φ i are in thermal contact with the SM bath during BBN, this restricts n φ < 2 (1) if φ is real (complex). 9 Note that this constraint is weaker if φ i decouples earlier.
When light elements are formed around 0.1 MeV, φ can dissociate the nuclei if it gives a recoil energy larger than the nuclear binding energy. However, due to the lightness of φ, the energy that φ can give to a nuclei is rather small. For example, in the 2 H rest-frame, the maximum recoil energy of the 2 H nucleus being hit by a φ is E max R = 2 E 2 φ /m2 H . This is only larger than the 2 H binding energy of 2.2 MeV if E φ > 47 MeV, which is much higher than the φ temperature at the same time. The effect on nuclear number densities is negligible.

Structure Formation
During the structure formation era (around a temperature of 10 eV), the scattering length between φ and 4 He was about 3 × 10 4 Mpc for Λ = 10 TeV, so we can treat φ as a collisionless particle. φ therefore generates a Landau damping to the primordial density fluctuations, with a free-streaming length that can be estimated as [27] λ F S, φ 20 Mpc m φ 10 eV −1 .

(5.9)
This is close to free-streaming neutrinos with λ F S, ν 20 Mpc mν 30 eV −1 , and φ should satisfy similar constraints as a thermally produced sterile neutrino, with cold dark matter still dominating relic density. As discussed in Section 5.2.1, this latter requirement of a sub-dominant hot dark matter φ component requires m φ eV. The scenario is then similar to the case studied by [37]. The existence of sterile neutrinos at sub-eV scale can relax the tension between Planck result and the local measurements of galaxy clusters on matter perturbation and the expansion rate of the universe. Similar conclusions apply to scalars, meaning sub-eV scale φ's are compatible with structure formation bounds.

Dark Acoustic Oscillations
When DM particles χ couple to a bath of nearly massless φ scalars, we expect the DM-φ system to give rise to dark acoustic oscillations (DAO), similar to the acoustic oscillations of baryons. The temperature and polarization spectra obtained from CMB data strongly constrain this effect, which translates to an upper bound on χ-φ scattering.
In dmDM, the only tree-level χ-φ interaction that is not suppressed by m χ is the process χφ → χ c φφ. This is mediated by a t-channel φ and is generated both by the DM yukawa coupling to φ and the quartic coupling λφ 4 in eq. (2.1). The transverse cross section of this process is only suppressed by the energy transfer m 2 χ v 4 χ and decouples at a very late time for light φ. Therefore, for scalars with mass below about 10 eV (the temperature of structure formation), CMB data sets stringent upper bounds on the coupling combination λy χ to ensure DAO do not generate a sizable effect.
Although the coupling between four of the lightest scalars has no direct implication for direct detection signals, it is still useful to understand this constraint for completeness. A detailed analysis using CMB data is beyond the scope of this work (for an example of a analysis, see [38]). However, we can estimate a conservative bound on y χ λ by requiring the scattering to decouple before structure formation (T ∼ 10 eV).
There are two ways for the χφ → χ c φφ process to freeze out before structure formation.
1. The lightest scalar could have a mass above 10 eV. To avoid overclosure, eq. (5.5) then requires φ to freeze out when g * S 10 2 , i.e. before the electroweak scale. This can be the case if Λ is very large. Indeed, for n φ = 1, this is required by neutron star cooling, see eq. (4.38). However, such a scenario would be sterile with respect to direct detection.

2.
In the next section we will define a n φ = 2 model which avoids neutron star bounds while allowing for direct detection. In that case, the lightest scalar must have a mass below ∼ eV to avoid overclosure. The light scalar therefore remains in thermal contact with dark matter, and remains relativistic during structure formation, which translates to a strong constraint on its quartic coupling.
The extent to which DM motion is influenced by φ scattering is given by the transverse cross section for χφ → χ c φφ, which we can estimate as The logarithmic factor comes from the Coulomb potential of the long-range φ interaction, and the additional phase space suppression of emitting an additional scalar is approximated by a factor of (16π 2 ) −1 . Assuming the scattering rate to be smaller than Hubble before 10 eV, the upper bound on the couplings translates to  for benchmark parameters m χ = 10 GeV, m φ L = 1 eV, and DM with kinetic energy ∼ 10 eV.
A more detailed study may relax the rather conservative bound, but a sizable tuning with λ ∼ 10 −8 is expected for y L χ ∼ 10 −4 . However, this bound has no bearing on direct detection.

A Realistic dmDM Scenario for Direct Detection
A summary of our derived constraints on light dark mediators and their coupling to the SM, formulated for n φ ≥ 1, is shown in Table 4. To place these in context, recall from Section 3 that the dominant Yukawa coupling should be y χ ∼ 10 −3 − 10 −2 . Furthermore, as we review in Section 7, direct detection of dmDM in the 2 → 3 regime is feasible if Λ ij 10 4 TeV with m φ i keV and m φ j MeV, so that one scalar can be emitted while the other acts as a light mediator.
With this in mind it is clear that any n φ = 1 scenario of dmDM with realistic direct detection prospects is completely excluded neutron star bounds. In fact, the minimum value of Λ required by neutron star cooling truncates the length of the supernova neutrino burst, so the actual lower bound on Λ becomes 10 11 TeV.
However, there is a very simple n φ = 2 scenario which behaves almost identically to the minimal n φ = 1 model for purposes of direct detection, yet is not excluded by any of the bounds in Table 4.
Consider a dmDM setup like eq. (2.1) with two light dark mediators, real scalars φ L and φ H having masses m φ L eV and m φ H ∼ MeV. We also add a quartic coupling to allow φ H to decay into φ L : Other quartic couplings are omitted for simplicity 10 . When λ > 10 −9 , φ H → φ L φ L φ L is instantaneous when the temperature drops below one MeV, leaving φ L with a similar relic density to eq. (5.5). Now let Λ LL > 10 8 TeV to comply with neutron star bounds, while Λ HH , Λ LH < 10 6 TeV avoids supernova bounds by trapping both φ L and φ H in the stellar medium. In that case, all the bounds in Table 4 are satisfied. Importantly, Λ LH , which can be relatively small, now controls direct detection. This can give a large rate for the processχN →χN φ L . Since φ H is much lighter than the typical momentum exchange of 10 MeV for ambient DM scattering off nuclei, the nuclear recoil spectrum is nearly identical to the n φ = 1 case with m φ ∼ eV.
A small Λ LH will generate an effective Λ LL coupling through a loop of constituent quarks and φ H . The size of this effective operator is where m q ≈ 263 MeV is the constituent quark mass. The neutron star bound on Λ LL then translates to a bound on Λ LH : Λ LH 10 TeV (6.3) which is the bound we adopt when discussing direct detection in the next section. Finally, the presence of two Yukawa couplings to dark matter gives additional freedom. A very modest hierarchy y H χ /y L χ 10 would suppress the χN → χN loop process. This makes it possible for y H χ to be large enough for a thermal relic χ and ameliorate the inconsistencies between dwarf galaxy simulations and observation, all while being in the 2 → 3 regime of direct detection (see Fig. 3).
This scenario can be realized in the UV completion of Section 2.2 by assuming hierarchical Yukawa couplings between dark mediators, dark vector quarks and different chiralities of the SM quarks. 10 The quartic coupling for φ 4 L would have to obey the constraint from dark acoustic oscillations, see eq. (5.11).

Direct Detection of dmDM
In this section we outline in detail our computation of nuclear recoil spectra and direct detection constraints on dmDM, first summarized in [10]. We work with the minimal n φ = 1 scenario with effectively massless φ for simplicity, with the understanding that this phenomenology can be replicated by the unexcluded n φ = 2 scenario defined by eq. (6.1). The discussion of the previous two sections derived constraints on the Yukawa couplings between dark mediators and dark matter, and the coupling between dark mediators and SM quarks. Direct detection is sensitive to a combination of the two. We predict the dmDM signal at XENON100 [6], LUX [4], CDMS-Si [5] and CDMSlite [66] and demonstrate that large regions of the direct detection plane are not yet excluded.
It is instructive to compare the dmDM interaction with nuclei to the contact operator since it is the standard choice for showing constraints from different direct detection experiments in the same (m χ , σ n SI )-plane. Referring to the above interaction model as the "standard-WIMP", we find that O(100 GeV) dmDM will fake a different lighter O(10 GeV) WIMP at different experiments. This is due to energy loss from the outgoing φ, which leads to an underestimate of the DM energy when assuming the above contact operator. We study this interesting phenomenon by first examining at the parton level cross section, understanding the parametric dependence of the recoil spectrum, and then produce the full experimental recoil prediction including form factors and the velocity distribution.

Differential cross section calculation
We are interested in the differential cross section for a dark matter particle hitting a stationary nucleus which then recoils with kinetic energy E r . This is given by dσ bare N /dE r is the 'parton-level' differential cross section evaluated for the process qχ → qχφ or qχ → qχ with the substitution of m q → m N . This is because ambient dark matter is extremely non-relativistic with velocites of order a few 100 km/s, interacting with the entire nucleus coherently. dσ bare N /dE r is easily evaluated analytically for the 2 → 2 loop process using eq. (2.2), reproducing the result of a standard WIMP with an additional suppression at high momentum transfer. For 2 → 3 scattering we adopt a Monte-Carlo approach 11 by defining a MadGraph5 [35] model containing the DM Yukawa coupling and theN N φφ * /Λ effective operator using FeynRules1.4 [67]. We discuss the resulting spectrum below.  Figure 13. Examples of nuclear recoil spectra with dmDM at 'parton-level' (without nuclear/nucleus form factors and coherent scattering enhancement) for different m N , m χ and a given incoming energy E χ . The blue datapoints are given by the MG5 simulation, and the red curve are the analytical approximation of the spectrum in Eq. eq. (7.5).
The factor (ΣB) 2 is a quark-nucleon form-factor to convert the amplitude from quark-to nucleon-level by taking into account the values of quark currents inside the proton or neutron (see [68,69] for a review). Since the momentum transfer q 2 is much less than the QCD confinement scale we can take this form factor to be constant. The relevant case for dmDM is the scalar operator N |m qq q|N = f N q m N , which is interpreted as the contribution of quark q to the nucleon mass m N . Importantly, the contribution of all sea quarks is additive, giving a large matrix element enhancement. f N q < 1, since each sea quark contributes more than its bare mass to the proton mass, and can be computed from lattice techniques. This gives matrix elements N |qq|N ≡ B N q , where Going from nucleon-to nucleus-level, the cross section is enhanced by A 2 (assuming equal φ coupling to protons and neutrons) and must be convolved with the Helm Form Factor [70,71]. This is just the Fourier transform of the radial nuclear density distribution, where j 1 is a Bessel Function, q = √ 2m N E r , s = 1 fm, r 0 = √ r 2 − 5s 2 and r = 1.2A 1/3 . It is instructive to compare dσ bare N /dE r for the 2 → 3 scenario to the simple WIMP case generated by the contact operator eq. (7.1). We examine the case of massless emitted φ. m φ ∼ few keV could be interesting to introduce shape features into the recoil spectrum but is cosmologically disfavored, see Table 4. As shown in Fig. 13, the recoil spectrum of dmDM can be well described by the function is the maximum allowed nuclear recoil energy for a given incoming DM velocity, same as for the standard WIMP. The above approximation holds for massive dark mediators as well, provided the intermediate t-channel φ has mass MeV and the emitted φ has mass keV. Eq. 7.5 can be decomposed into the phase space part of χ N → χ N φ scattering via a contact interaction, times the propagator of the light mediator φ. The phase space can be approximated by which vanishes when E R reaches its maximum value, or when E R → 0, since a relativistic φ itself cannot compensate both the energy and momentum of a non-relativistic DM particle. The non-relativistic scattering also requires the spatial momentum exchange to be much larger than the kinetic energy, which makes the spatial momentum of the relativistic φ negligible in energy-momentum conservation. Because of this, the propagator of φ in the dmDM scattering is (2m N E R ) −2 , which is dominated by the spatial momentum exchange between N and χ and gives the spectrum in eq. (7.5). The contact operator eq. (7.1) produces a flat parton-level nuclear recoil spectrum for E r < E max r . On the other hand, eq. (7.5) features a suppression at large recoil. The functional form of this recoil suppression is different than for 2 → 2 scattering with light mediators and/or derivative couplings. Furthermore, the scaling of total cross section with m N , m χ is unique. This necessitates a full re-interpretation of all direct detection bounds to understand how a heavy dmDM candidate fakes different light WIMPs at different detectors. We expect the recoil suppression to increase the sensitivity advantage enjoyed by low-threshold Xenon detectors over CDMS.

Nuclear Recoil Spectra
To compute the expected nuclear recoil spectrum at a direct detection experiment, the differential scattering cross section must be convolved with the dark matter speed distribution in the earth frame,  Figure 14. The fraction of the dmDM (solid) and WIMP DM (dashed) direct detection cross section above experimental threshold (blue for CDMS II Si E r > 7 keV, black for LUX S1 > 2) as a fraction of the total cross section. m φ < keV. S1 light collection efficiency is taken into account but signal selection cuts have not been applied. The speed distribution is given in Appendix C. In our Monte Carlo calculation for 2 → 3 scattering, we simulate χN → χ c N φ for different m N , m χ , m φ and incoming DM velocities v in bins of 20 km/s to build up a table of the various required dσ N dEr and perform this convolution numerically. For verification, we applied our pipeline to WIMP-nucleus scattering, reproducing the expected analytical results.
The spectrum of nuclear recoil events that occurred in the detector must be translated to actual experimental observables. This involves folding in efficiencies, as well as converting the nuclear recoil signal to a scintillation light signal in the case of liquid Xenon detectors. These details are also given in Appendix C.  Figure 17. The measured nuclear recoil spectrum produced by a dmDM candidate with mass m χ = m 2→3 is very similar to that of a WIMP with mass m 2→2 < m 2→3 , interacting with nuclei via the contact operator eq. (7.1). m 2→2 (m 2→3 ) is shown for XENON100 (S1 > 3 with 6% light gathering efficiency, dashed red line), LUX (S1 > 2 with 14% light gathering efficiency, dash-dotted black line), CDMS II Silicon (E r > 7 keV, solid blue line), and CDMSlite (Germanium, Er > 0.2 keV, dotted purple line) before selection cuts.
The detection efficiency for 2 → 3 scattering in dmDM is about 100 − 1000 times smaller compared to the standard WIMP (and also 2 → 2 scattering in dmDM), see Fig. 14. This is expected, given the additional E r -suppression. In the next subsection we will take care to understand the parameter regions where 2 → 3 scattering dominates over the 2 → 2 process in dmDM Fig. 15 shows some 2 → 3 nuclear recoil spectra at CDMS II Si before taking detection efficiency into account. dmDM is compared to WIMPs for different DM masses, and the principal experimental feature of our model is apparent: a ∼ 50 GeV dmDM candidate looks like a ∼ 10 GeV WIMP. Fig. 16 shows different S1 spectra at LUX, where a ∼ 50 GeV dmDM candidate looks more like a ∼ 20 GeV WIMP. This mass remapping compared to the standard contact operator interpretation is shown for different experiments in Fig. 17. This dependence of recoil suppression on the detector and DM parameters is unique to dmDM, and could be added to other DM models by including the emission of a light particle..

Direct Detection Constraints
We compute direct detection bounds on dmDM in two ways. The first is by remapping the bounds provided by the respective experimental collaborations using the remapping of dmDM to standard WIMP parameters [10], part of which is shown in Fig. 17. These results are then reproduced, for verification, by using a full modified maximum likelihood analysis [72] for each experiment. The resulting bounds in the direct detection plane for dominant 2 → 3 and 2 → 2 scattering in dmDM are shown in Fig. 18 and Fig. 19. To provide a lower boundary on the relevant parameter space we indicate where in the direct detection plane the dmDM signal gets drowned out by the irreducible neutrino background [73].
For the 2 → 3 and 2 → 2 scattering regimes, direct detection probes y χ /Λ and y 2 χ /Λ respectively. The neutron star cooling bound eq. (6.3) for the n φ = 2 model Λ 10 TeV and the bounds on dark matter Yukawa coupling y χ can be combined to be shown in the direct detection planes of Figs. 18 and 19. The assumption of a thermal relic then sets bounds which supersede the liquid Xenon experiments for m χ 10 GeV.
For n φ = 1, the 2 → 2 loop process in Fig. 1 dominates if y χ 10 −3 . This is indicated, together with the neutron star bound, by the dashed orange line in Fig. 18. However, in the n φ = 2 model of Section 6, the vertical axis of Figs. 18 and 19 is (y H χ /Λ) 2 and (y H χ y L χ /Λ) 2 respectively, so this orange line can be moved arbitrarily upwards. This means the n φ = 2 model can realize 2 → 3 dominated direct detection while being consistent with a thermal relic, as well as the SIDM solution to the inconsistencies between dwarf galaxy simulations and observation.

Conclusion
Previous theoretical investigations have shown that direct detection can proceed very differently from the standard WIMP scenario. Investigating all possibilities is of vital importance. For one, the current list of experimental anomalies naively conflicting with other collaborations' bounds motivates the search for alternative interpretations of the data. Another general reason for achieving 'full theoretical coverage' is the looming irreducible neutrino background [73] that direct detection could become sensitive to in about a decade. Hitting this neutrino floor without a clear dark matter signal is an undesirable scenario, but being left without alternative options to explore would be an even more dire situation.
Dark Mediator Dark Matter is the first example of a slightly non-minimal dark sector where the mediators connecting dark matter to the Standard Model are themselves charged under the same symmetry that makes dark matter stable. Phenomenologically, this closes a long-standing gap in the list of investigated scenarios by realizing 2 → 3 nuclear scattering at direct detection experiments.  Figure 18. Direct detection bounds on the 2 → 3 regime of n φ = 1 dmDM. The vertical axis is proportional to σ χN →χN φ , and is understood to be (y H χ /Λ) 2 for the n φ = 2 model of Section 6. Solid lines: 90% CL bounds by XENON100 (red), LUX (black) and CDMSlite (purple), as well as the best-fit regions by CDMS II Si (blue, green). The large-dashed black line indicates where the dmDM signal starts being drowned out by the irreducible neutrino background [73]. Small-dashed magenta line: y χ = y relic χ (m χ ) and Λ = 10 TeV. Need to be below this line for a thermal relic to be compatible with the neutron star cooling bound eq. (6.3). Lower dotted orange line: for n φ = 1, below this line y χ is small enough to ensure the 2 → 3 process dominates direct detection while also satisfying the neutron star cooling bound. This line can be arbitrarily moved when n φ = 2.
We carry out the first systematic exploration of light scalar mediators coupling to the SM quarks via operators of the formqqφφ * /Λ. Their existence and coupling can be strongly constrained by cosmological bounds, LHC direct searches and stellar astrophysics, see Table 4. Neutron star cooling excludes detectable dmDM scenarios with a single dark mediator completely, but an n φ = 2 scenario can easily evade all bounds while giving identical direct detection phenomenology.
The presence of a light mediator and additional particle emission means that the nuclear recoil spectrum of dmDM at direct detection experiments is strongly peaked towards the origin. The functional form of this recoil suppression and the overall cross section dependence on nucleus and DM mass is unique. As a consequence of this suppression, a ∼ 100 GeV dmDM candidate fakes different O(10 GeV) standard WIMPs at different experiments. We compute direct detection bounds on dmDM for both nuclear scattering processes, χN → χN φ and the loop suppressed χN → χN and find large regions that are not excluded but discoverable in the future. The abovementioned n φ = 2 scenario can realize 2 → 3 direct detection while being compatible with a thermal relic and the SIDM solution for the inconsistencies between dwarf galaxy simulations and observation.  Figure 19. Direct detection bounds on the 2 → 2 regime of n φ = 1 dmDM. The vertical axis is proportional to σ χN →χN , and is understood to be (y H χ y L χ /Λ) 2 for the n φ = 2 model of Section 6. Same labeling as Fig. 18.
Our model represents an interesting combination of light mediator and inelastic scattering ideas, since the latter is realized by having a light scalar φ from a direct-detection point of view. This allows us to smoothly map dmDM spectra to similar WIMP spectra, and the resulting map of dmDM parameters to WIMP parameters makes transparent how the direct detection bounds are re-interpreted (see also [10]). While dmDM does not reconcile the conflicting signals and constraints, it may point the way towards another model that does. For example, it might be interesting to explore how this new scattering process changes models with non-standard form factors or exothermic down-scattering.   [75]. The profiles for our benchmark white dwarf were produced with the MESA code [44] by Max Katz.

A Radial Profiles for Stellar Cooling Calculation
The solar energy loss and radiative heat transfer calculation in Section 4.3 makes use of standard radial profiles for temperature, density and composition of the sun, shown in Fig. 20. These can be found in basic astrophysics textbooks like [74]. The radius and power output of the sun are R sun ≈ 3.85 × 10 26 cm and P sun ≈ 3.85 × 10 26 Watts.
For the white dwarf cooling calculation in Section 4.4 we gratefully acknowledge the help of Max Katz, who simulated the evolution of an approximately one solar mass sun-like star from the main sequence to a very old white dwarf using the MESA stellar evolution code [44]. The mass of this white dwarf, ≈ 0.5 solar masses, is representative of the majority of white dwarfs in the luminosity function dataset [40,43]. The luminosity (power output in solar units) and core temperature of this dwarf over time are shown in Fig. 21, along with the relationship between core temperature and power output.
Radial density, temperature and composition profiles for this white dwarf when its power  Figure 21. Top two plots: Evolution of power output and core temperature for the one solar mass white dwarf simulated by Max Katz using MESA [44]. Bottom: relationship between power output an core temperature. output was P WD /P sun ≈ 0.1 are shown in Fig. 20. This point in the evolution of the dwarf marks the start of dominant photon cooling [41], which we compare to φ emission in Section 4.4.

B Checking other constraints on φ
Here we briefly demonstrate that fixed target experiments, precision measurements bounds, and indirect detection do not constrain dark mediators.

B.1 Fixed target experiments
The fixed target experiments MINOS (E p = 120 GeV) [76,77], T2K (E p = 30 GeV) [78][79][80], MiniBooNE (E p = 8.9 GeV) [81] and LSND (E p = 0.8 GeV) [82] bombarded graphite or beryllium targets with ∼ 10 20 − 10 23 protons-on-target (N pot ). Dark mediators can be produced in these collisions via the process pp → ppφφ, which has cross section σ produce ∼ 10 −5 to 10 −3 pb for Λ = 10 TeV, depending on E p . (This was computed in MadGraph5.) We can estimate whether these experiments are sensitive to φ production with this cross section. Hitting a target of some length and proton number density L target , n target with N POT protons produces the following number of φ's: For Λ = 10 TeV, the interaction cross section of a high-energy φ with stationary protons is σ pφ→pφ ∼ 0.1 pb over the relevant range of proton energy E p . This makes the mean free path of a high-energy φ in a typical material with densities ∼ g cm −3 about ∼ 10 10 meters, so we can assume all φ's leave the target. The chance that a single φ is detected is where (likely 1) is some efficiency factor to account for the finite φ-beam width as well as likelihood of detecting the target nuclear recoil. n detect and L detect are the number densities Figure 22. One of the penguins in the K + → π + ν ν and the corresponding K + → π + φ φ process. and lengths of the detector material, and σ detect ∼ σ pφ→pφ is the interaction cross section of φ with the target material. Therefore, the number of φ's detected by this experiment is Substituting graphite and iron densities for target and detector respectively, as well as the production cross section at E p = 150 GeV and = 1 to strongly overestimate φ detection, we get N detected  10m). Therefore, fixed target experiments have no sensitivity to φ production for Λ 10 TeV.

B.2 Precision Measurement Bounds
Since φφ * has potentially sizable coupling to the scalar quark current it could contribute to meson decays and the invisible Z-width. We show below that no meaningful constraints are derived from these processes. However, the heavy vector-like quarks in the UV-completion do contribute to the S-parameter, which bounds their Yukawa coupling to the Higgs.

B.2.1 Bound from the meson decays
Pions are the lightest QCD pseudoscalars and cannot decay to φφ * , but Kaon decay with φφ * in the final state is possible. The measurement [83] Br(K + → π + νν) exp = 17.3 +11.5 −10.5 · 10 −11 , (B.5) is quite close to the SM value [84] Br(K + → π + νν) SM = (8.5 ± 0.7) · 10 −11 . (B.6) Since K + → π + φφ gives the same signal as the neutrino process, the same measurement sets a lower bound on the suppression scale Λ. The dominant diagrams involve a W -up-quark loop which suffers from a GIM mechanism and requires up-type mass, dominated by the top. The leading order SM result is given by [85] and we can use it to estimate the dmDM matrix element: This is clear from Fig. 22 by chirality and dimensionality arguments. For Λ = 10 TeV the ratio is about 0.03, much less than the current experimental precision. Kaon decay therefore supplies no meaningful dmDM bounds.

B.2.2 Electroweak Precision Measurement
As explained in Section 2.2, theqqφφ * /Λ operator is most plausibly generated by a single generation of vector-like quarks. Since these quarks have SU (2) L × U (1) Y charge they contribute to the oblique parameters S, T, U [86]. The resulting constraints on such a Top Partner Doublet model have been computed by [87]. For M Q ≈ 1 TeV, the mass splitting between the up and down type vector-like quarks must be less than about 10 GeV, which sets a strong bound on the flavor structure. However, for flavor-diagonal couplings there are no constraints.

B.2.3 The Z φ φ coupling
The properties of the Z boson are extremely well measured. The |φ| 2Q q/Λ operator contributes to the invisible Z-width. Of all the electroweak precision constraints, this contribution gives the strongest bound on dmDM. Assuming flavor-diagonal |φ| 2 coupling to quarks in the SM mass basis, the dominant contribution to Z → φφ * comes from a top loop with a single mass insertion. This can be expressed as an effective Zφφ * coupling This is much smaller than the current precision of Γ(Z → invisible) = 499.0 ± 1.5 MeV [83], meaning electroweak precision constraints supply no meaningful bounds on dmDM.

B.3 Indirect Detection
A potential indirect detection signal may arise from the annihilation process χ c χ → qqφ * (assuming the 2 → 2 loop processχχ →qq is suppressed like in direct detection). The total annihilation cross section is found using MadGraph5 to be The Fermi-LAT collaboration [88] reports dwarf galaxy bounds 12 on (σ χχ→qq v) ann between ∼ 10 −26 and 10 −25 cm 3 s −1 for m χ between 2 GeV and 100 GeV. The dmDM annihilation cross section is many orders of magnitude smaller. Furthermore, the Fermi bounds assume 2 → 2 annihilation, resulting in a monochromatic quark spectrum E q = m χ before hadronization and decay. In dmDM the spectrum is triangularly rising towards E q = m χ , which is harder to detect. It is therefore clear that dmDM leaves no detectable signal in the gamma ray sky.

C.1 Dark matter speed distribution
The velocity of thermalized cold dark matter in the halo frame has an approximate Maxwell-Bolzmann distribution: where v 0 ≈ 220 km/s and v esc ≈ 544 km/s are the temperature of the distribution and the galactic escape speed 13 [89]. Transforming this distribution to the earth frame moving at an average velocity of v e ≈ 233 km/s through the galaxy gives the relevant speed distribution for our calculations.
e 2vve/v 2 0 − e 2 cos φmaxvve/v 2 0 , 12 For these objects there is an independent handle on the local DM densities by means of measuring the peculiar velocity, thereby significantly reducing dependence on any particular halo model compared to the galactic center. 13 We implement the galactic escape speed with a hard cutoff, which is not entirely realistic, but the effect of vesc < ∞ is very small in the earth frame so this is sufficient for our purposes.
where η is a normalization factor, and cos φ max as a function of v is given by

CDMS II Silicon
The CDMS II Silicon detectors (m N = 38 GeV) have accumulated 140.2 kg·days of expsure and use simultaneous measurement of ionization and non-equilibrium phonons to measure nuclear recoil and distinguish from electron recoil background. The recoil cutoff is 7 KeV (higher for some sub-detectors), which together with the fiducial volume and phonon timing cuts (to eliminate background) results in the WIMP-nucleon scattering efficiency curve shown in [5]. To obtain experimental predictions for recoil spectra we simply multiply each recoil bin by this efficiency, which is about 20% for E r = 10 keV and asymptotes to about 40% at 30 keV.

C.2.2 XENON100 and LUX
Xenon100 [6] and LUX [4] (m N = 131 GeV) have accumulated exposures of 7636 and 10065 kg·days and detect nuclear recoil with two experimental signals: the number of produced scintillation photons (S1) and the charge signal once the ionization travels up to the gaseous phase (S2). Nuclear recoils can be distinguished from electron recoil backgrounds using the time difference between the S1 and S2 signals, as well as their ratio S2/S1. XENON100 has a light gathering efficiency of about 6% for S1 photons. The expected number of scintillation photons for a given nuclear recoil event is [6] S1 = E r × S nr L y S ee L ef f (E r ), (C.5) where S ee = 0.58, S nr = 0.95 and L y = (2.28 ± 0.04) (photo electrons)/(keV ee ). The relative scintillation efficiency L ef f must be extracted from experiment. The fit used by the collaboration can be found in [90]. Due to limitations of the L ef f measurement, no DM signal below E r = 3 keV is included. This makes the bounds conservative. The same applies for LUX, appropriately rescaled to account for the highter 14% light gathering efficiency.
The expected S1 spectrum can then be computed from the expected nuclear recoil spectrum, dN dS1 = dE r dN dE r Poi(S1, S1 ), (C. 6) where S1 is a function of E r as per eq. (C.5). Once an S1 signal has been detected it must pass selection cuts. The probability of a WIMP signal passing these cuts is S1-dependent and asymptotes to about 0.8 for S1 5 at XENON100 and about 1 for S1 2 at LUX. The XENON100 signal region is S1 ∈ (3, 30), while for LUX it is S1 ∈ (2, 30).

C.2.3 CDMSlite
CDMSlite [66] was a light DM search using a single Super-CDMS iZIP detector operated at higher bias voltage to lower the nuclear recoil detection threshold to 0.84 keV at the cost of giving up background discrimination. The accumulated exposure is 6.18 kg·days. Instead of computing the dmDM CDMSlite signal with experimental efficiencies, we transform the collaboration's DM bounds using the dmDM→WIMP parameter map discussed in Section 7.3, having validated the method on XENON and CDMS-Si data.