Direct detection and complementary constraints for sub-GeV dark matter

: Traditional direct searches for dark matter, looking for nuclear recoils in deep underground detectors, are challenged by an almost complete loss of sensitivity for light dark matter particles. Consequently, there is a signi(cid:12)cant e(cid:11)ort in the community to devise new methods and experiments to overcome these di(cid:14)culties, constantly pushing the limits of the lowest dark matter mass that can be probed this way. From a model-building perspective, the scattering of sub-GeV dark matter on nucleons essentially must proceed via new light mediator particles, given that collider searches place extremely stringent bounds on contact-type interactions. Here we present an updated compilation of relevant limits for the case of a scalar mediator, including a new estimate of the near-future sensitivity of the NA62 experiment as well as a detailed evaluation of the model-speci(cid:12)c limits from Big Bang nucleosynthesis. We also derive updated and more general limits on DM particles upscattered by cosmic rays, applicable to arbitrary energy- and momentum dependences of the scattering cross section. Finally we stress that dark matter self-interactions, when evaluated beyond the common s -wave approximation, place stringent limits independently of the dark matter production mechanism. These are, for the relevant parameter space, generically comparable to those that apply in the commonly studied freeze-out case. We conclude that the combination of existing (or expected) constraints from accelerators and astrophysics, combined with cosmological requirements, puts robust limits on the maximally possible nuclear scattering rate. In most regions of parameter space these are at least competitive with the best projected limits from currently planned direct detection experiments.


Introduction
So far no unambiguous signal for new physics at the electroweak scale has been identified at the Large Hadron Collider (LHC) [1,2], despite seemingly intriguing theoretical arguments that have been brought forward why the appearance of new physics should be expected at these energies (with low-scale supersymmetry being the most popular example, see e.g. [3]). In consequence, while some natural islands remain [4,5], the experimental focus in the search for physics beyond the standard model of particle physics (SM) presently undergoes a substantial broadening in scope, both concerning energy scales and theoretical frameworks for such searches [6]. At the intensity frontier, in particular, there is a plethora of both ongoing and planned activities that aim to explore new physics in the sub-GeV -1 -JHEP03(2020)118 evolution of the dark sector in section 5. In that section we also derive bounds from BBN and DM self-interactions that apply to the specific scenario considered here, and mention further astrophysical bounds. In section 6 we then combine the various constraints, and compare them to (projected) bounds from direct DM detection experiments. Finally, in section 7 we discuss our results and conclude.
2 Models for light dark matter with portal couplings Light dark sector particles are required to have small couplings to SM states in order to be allowed phenomenologically and therefore naturally correspond to fields that are singlets under the SM gauge interactions. They may then directly couple to the SM via the well known portal interactions [33], i.e. gauge-invariant and renormalisable operators involving SM and dark sector fields. If the DM particle χ is stable and fermionic, as assumed in this work, no direct renormalisable interaction is available and an additional particle X mediating the interactions with the SM is required. 1 In recent years mediator searches at colliders together with complementary constraints from direct detection have therefore received a large amount of interest, both for searches at the LHC [56][57][58][59][60][61][62] and at low energy colliders [63][64][65]. When comparing the sensitivities of collider searches with direct detection experiments it is important to take into account the large difference in energy scale between the centre of mass energy at colliders and the typical momentum transfer in nuclear recoils. In particular the relative sensitivity of direct detection experiments is significantly increased for light mediators, implying that while scenarios with heavy mediators are strongly constrained by collider searches, those constraints are significantly weakened for light mediators. Another appealing feature of light mediators, adding predictivity, is that DM can be produced within the standard thermal freeze-out paradigm: 2 For sufficiently large couplings the dark sector will thermalise in the early universe and the DM relic abundance is set via annihilations into two mediators, χχ → XX, but also via annihilations into SM fermions via an s-channel mediator, χχ → X →f f , if the dark sector is not fully decoupled at the time of freeze-out. Two particularly interesting and often studied options are vector mediators kinetically mixed with the SM hypercharge gauge boson or scalar mediators with Higgs mixing.

Vector mediators
Let us start with a brief discussion of the vector mediator case. In the simplest scenario the field content consists of only a dark matter fermion charged under a dark U(1) X with kinetic mixing (see e.g. [34,39,66]). For light mediators the coupling structure will basically be that of a photon, so that X predominantly decays to charged SM fermions such as electron 1 For scalar DM, on the other hand, there is a direct (Higgs) portal term [53], constituting the most minimal DM model that is phenomenologically viable; for a recent status update see ref. [54]. Another portal term exists for a new heavy neutral lepton mixing with the SM model neutrinos; such a particle (often called 'sterile neutrino') is not stable (decaying e.g. into three SM neutrinos), but can be sufficiently long-lived to constitute DM [55]. We do not consider these options here. 2 It has been noted that for heavy mediators DM overproduction can only be avoided in rather special corners of parameter space [60][61][62].

JHEP03(2020)118
positron pairs. An important observation is that DM annihilations proceed via s-wave for both channels discussed above. If DM was ever in thermal contact with the SM (not necessarily through the kinetic mixing) such that the dark sector temperature is not much smaller than the photon temperature, there are strong constraints from Cosmic Microwave Background (CMB) observations, ruling out DM masses m χ 10 GeV [67]. In fact these bounds extend to significantly higher DM masses for mediators parametrically lighter than DM due to the Sommerfeld enhancement of the annihilation cross section [68].
There are a number of ways to evade this CMB limit, but they do involve some nonminimal component in the DM model. For instance DM may be asymmetric with only a sub-leading symmetric component such that residual annihilations during CMB are sufficiently suppressed [69]. For consistency such a setup will however require the existence of an additional dark sector state to compensate the charge of the DM (reminiscent of electrons and protons). Another possibility would be to introduce a scalar whose vacuum expectation value (vev) generates small Majorana mass terms for the DM fermion, resulting in two dark matter states with slightly different mass, coupled off-diagonally to the vector boson (this is often referred to as inelastic DM) [70]. If the heavier state decays before the time of the CMB, s-wave annihilations χ 1 χ 2 → X →f f are no longer possible and constraints are evaded. A third possibility would be to couple the vector mediator to a light hidden sector state such that the decays of X are invisible [71], in which case the CMB bounds can also be evaded. Finally, if the abundance is set via freeze-in [72] rather than freeze-out, the annihilation cross section may be sufficiently small to be in accord with observations. While all these options are viable and possess an interesting phenomenology, we wish to concentrate on a minimal setup in the current study. As discussed below, a model for light DM which still survives in its simplest form is that of a scalar mediator with Higgs mixing.

Scalar mediators with Higgs mixing
In contrast to the case of a vector mediator, DM annihilations proceed via p-wave for a scalar mediator and the setup is correspondingly much less constrained by residual annihilations during CMB times. If the dark sector was in thermal contact with the SM heat bath at even earlier times, however, dark sector masses are typically still required to be larger than m χ 10 MeV in order not to spoil the agreement between predicted and observed primordial abundances of light nuclei (we will study the relevant limits in detail below).
Let us consider a new real scalar S that mixes with the SM Higgs and further couples to a new Dirac fermion χ that can play the role of the DM particle (see e.g. ref. [45] and references therein), Here m χ is the mass of the DM fermion, H is the Higgs doublet of the SM and V (S, H) is the scalar potential. The terms involving the singlet scalar can be written as with V (S) = ξ s S + 1/2µ 2 s S 2 + 1/3A s S 3 + 1/4λ s S 4 . Without loss of generality the field S can be shifted such that it does not obtain a vev, implying ξ s = A hs v 2 /2 (where the Higgs vev is given by 246.2 GeV). After electroweak symmetry breaking the singlet S mixes with the physical component of H such that the singlet S naturally acquires a coupling to all SM fermions while the Higgs h acquires a coupling to χ, with mixing angle (2.4) The usual Higgs quartic coupling λ h is fixed in the SM via the observed Higgs mass and we are interested in the parameter region m h m S . In our convention where S does not acquire a vev the mixing angle is therefore approximately fixed by A hs . While the mixing angle clearly has a very large impact on most of the experimental observables, it does not fully specify the phenomenology of the scalar sector. For example, the decay width of the SM-like Higgs boson into two light singlets is determined by the S 2 H † H coupling, When we evaluate constraints e.g. from the Higgs signal strength we will assume that λ hs 0 to be conservative. Similarly we do not rely on λ hs for the thermalisation of the SM with the dark sector, yielding conservative limits from BBN. We also assume the trilinear coupling A s to be small, so that the 3 → 2 annihilation rate of singlet scalars is negligible and no phase of 'cannibalism' [73] occurs after freeze-out, again leading to conservative bounds.
For the calculation of DM-nucleus scattering rates we will also need the effective Yukawa coupling between a nucleon and the scalar mediator, Here the constants f q,G correspond to the quark and gluon content of the nucleon. It is well known that the couplings to protons and neutrons are very similar for Higgs exchange with g n ≈ g p ≈ 1.16 · 10 −3 sin θ, using state-of-the-art values for the f q [74].
3 Constraints from direct dark matter searches

Conventional light dark matter detection
Direct detection experiments probe the elastic scattering cross section σ SI χN between DM particles χ and nuclei N (since we only consider scalar mediators, we restrict our discussion to spin-independent scattering) at finite (spatial) momentum transfer

JHEP03(2020)118
where T N is the nuclear recoil energy. For better comparison, however, these results are typically reported in terms of the inferred cross section per nucleon, σ SI , at zero momentum transfer. An assumption that is often adopted for the sake of this translation is that of isospin-conserving couplings, which is almost perfectly satisfied for a scalar with Yukawalike coupling structure. This leads to the familiar coherent enhancement of Here µ χN and µ χp are the reduced masses of the DM/nucleus and DM/nucleon system, respectively, and A is the atomic mass number of the nucleus N . Compared to this, the cross section at finite momentum transfer is suppressed by a nuclear form factor G N , This form factor is conventionally computed as the Fourier transform of the nuclear density profile, i.e. under the assumption that the scattering on the nucleons does not induce an additional momentum dependence. For an interaction mediated by a scalar S it is straightforward to calculate the nonrelativistic scattering cross section as [45] 4) where g N denotes the coupling between S and the nucleus, i.e. g N = Ag p if isospin is conserved. In the case of Higgs mixing, using eqs. (2.6), (3.2) and (3.4), this translates to the DM-proton scattering cross section For a heavy mediator, this expression can directly be compared to standard limits on σ SI because scattering on nucleons is essentially momentum-independent. If the mediator is light compared to the typical momentum transfer, however, the cross section probed in the detector is smaller than expected from eq. (3.3) and limits have hence to be re-evaluated taking into account all the relevant experimental information. An approximate -but still reasonably accurate -way of taking into account the momentum suppression consists in simply rescaling (see, e.g., ref. [75]) whereσ χN is the limit reported under the assumption of a constant scattering cross section -in terms of σ SI as given in eq. (3.2) -and Q 2 ref is an experiment-specific reference momentum transfer (see also ref. [76] for a discussion of how to explore light mediators with direct detection experiments).
In figure 1 we summarise the most stringent (projected) direct detection constraints at low DM masses, along with the value of Q 2 ref that we use for the corresponding experiment. The latter was either estimated by using eq. (3.1) for the minimal recoil energy adopted in the respective analysis, or by directly fitting to data provided by the experiment (for PandaX-II [25]). We note that carefully modelling inelastic scattering processes, resulting in the emission of a photon or an atomic electron, in principle allows to improve sensitivities in the few 100 MeV range [83][84][85]. There is also a number of proposed direct detection experiments, and ideas, that would probe even smaller cross sections in the mass range shown in figure 1, but the status of those is presently less certain (for a recent compilation, see refs. [35,48]).

Cosmic ray-accelerated dark matter
The right panel of figure 1 clearly illustrates the exponential loss of sensitivity of conventional direct detection experiments to sub-GeV DM, reflecting the fact that non-relativistic DM particles with such small masses do not carry enough momentum to allow for nuclear recoils above the experimental threshold. As recently pointed out, however, there is a small yet inevitable component of relativistic DM that alleviates this limitation [26]: 3 if DM can elastically scatter with nuclei, then also the well-established population of high-energy cosmic rays will scatter on DM, thus accelerating them from essentially at rest (in the galactic frame) to GeV energies and beyond -in principle for arbitrarily small DM masses.
In order to handle scattering via light mediators we extend the formalism developed in ref. [26] to allow for arbitrary relativistic scattering amplitudes (rather than only a constant σ χN as assumed there). As the derivation follows the same steps as in ref. [26], we only briefly state our results here and refer to that reference for further details (see also JHEP03(2020)118 ref. [94]). The flux of cosmic-ray accelerated DM (CRDM) before a potential attenuation in the Earth or the atmosphere is given by Here, ρ local χ and Φ LIS CR are the local interstellar DM density and the cosmic-ray flux, respectively, and T min CR is the minimal kinetic cosmic-ray energy needed to accelerate DM to kinetic energy T χ ; we take into account elastic scattering of cosmic-ray nuclei N = {p, 4 He} with DM, including in each case the same dipole form factor suppression as in ref. [26]. 4 D eff ∼ 8 kpc, finally, is an effective distance out to which we assume that the source density of CRDM is roughly the same as it is locally (which, for a standard DM distribution, corresponds to a sphere of about 10 kpc diameter). The scattering rate of relativistic CRDM particles in underground detectors is then determined as where the scattering cross section dσ χN /dT N must be evaluated for the actual DM energy T z χ at the detector's depth z (which is lower than the initial DM energy T χ due to soil absorption [95][96][97][98]), and T χ (T z,min χ ) denotes the minimal initial CRDM energy that is needed to induce a nuclear recoil of energy T N (again taking into account a potential attenuation of the flux due to the propagation of DM through the Earth and atmosphere). In order to relate T z χ to the initial DM energy T χ = T z=0 χ , we numerically solve the energy loss equation where T max N denotes the maximal recoil energy T N of nucleus N , for a given DM energy T z χ , and we sum over the 11 most abundant elements in Earth's crust. It is worth stressing that the momentum transfer in a direct detection experiment is given by eq. (3.1) also in the relativistic case. In particular, the form factor in the nuclear scattering cross section does not depend on the energy of the incoming DM particles, only on the relatively small range of Q 2 that falls inside the experimental target region. This makes it straightforward to translate direct detection limits reported in the literature for heavy DM, assuming the standard DM halo profile and velocity distribution, to a maximal count rate in the analysis window of recoil energies and in turn to limits resulting from the CRDM component discussed here [26]. The updated routines for the computation of the resulting CRDM flux and underground scattering rates have been implemented in DarkSUSY [99], which we also use to calculate the resulting limits from a corresponding re-interpretation of Xenon-1T [24] results. Direct detection constraints on dark matter accelerated by cosmic rays for fixed mediator masses. Cross sections below the lower boundaries lead to recoil rates too small to be detectable, while cross sections above the upper confining boundaries prevent the dark matter particles to reach the detector, due to efficient scattering in the overburden. As a rough indication of how large cross sections are in principle possible, we also show in each case the parameter range where the couplings are well inside the perturbative regime (for a more detailed treatment, see ref. [100]). Right panel. Same, for fixed mediator to DM mass ratios.
In order to do so, we still need the full relativistic scattering cross section of DM with nuclei, mediated by a scalar particle. For fermionic nuclei we find where σ SI,NR χN is the scattering cross section in the highly non-relativistic limit, as stated in eq. (3.4), s = E 2 CM and G N (Q 2 ) is the conventional nuclear form factor. While the non-relativistic result is of course recovered for Q 2 → 0 and s → (m χ + m N ) 2 , this cross section is actually enhanced for Q 2 m 2 χ when compared to the standard estimate given in eq. (3.6). This is particularly relevant both for very light DM (m 2 χ Q 2 ref ) and the production of the CRDM component stated in eq. (3.7), for which the momentum transfer is typically much larger than expected in underground experiments.
In figure 2 we show the resulting limits from Xenon-1T on light DM. An important feature of a constant scattering cross section is that these constraints (almost) flatten for very small DM masses [26]. Compared to that, as expected from the above discussion (see also ref. [94]), we observe a significant strengthening of our constraints at fixed mediator masses. However the figure also clearly demonstrates that for light mediator masses the production of the CRDM component becomes suppressed by the mediator momentum; when considering only mediators that are lighter than the DM particle, in particular, the resulting constraints become less and less stringent. Also the behaviour of the maximal cross section (due to soil absorption) is rather instructive, as it falls into two clearly distinguishable regimes: i) for heavy (GeV-scale and above) mediators the upper boundary essentially follows that of the constant cross section case [26], roughly rescaled by an additional m 2 χ dependence (for small m χ ) with the same origin as discussed above for the lower -9 -

JHEP03(2020)118
boundary; ii) for lighter mediator masses, the momentum suppression starts to become relevant, strongly favouring scattering events in the overburden with small momentum transfers -which in turn leads to a significantly increased penetration depth, and hence weaker constraints. Let us, finally, stress that the limits presented in figure 2 in principle apply to any model with scalar mediators, i.e. they are not restricted to the specific structure of the DM-nucleon coupling given in eq. (2.6).

Constraints from particle physics experiments
Let us now turn to constraints on the scalar portal model from particle physics experiments. In the following we concentrate mostly on the case m S m χ so that the annihilation channelχχ → SS is kinematically allowed in the early universe. The reason is that for m S m χ only direct annihilations into SM states via an s-channel scalar singlet are allowed,χχ → S → SM (see section 5.1 for a more detailed discussion). The corresponding annihilation rate, however, is typically constrained to be too small to allow for the observed DM relic abundance (see e.g. [45]), making this case less appealing. Note that m S m χ naturally implies that the singlet scalar S can only decay to SM states that can potentially be observed in detectors ('visible decays'). Depending on the mixing angle θ, however, the lifetime of S can be so long that the decay happens only outside of the detector and the signature is therefore identical to an invisibly decaying scalar. While we mainly concentrate on this case, we will also briefly comment on the case m S m χ .
An important property of the inherited Yukawa-like coupling structure is that the production of S may well proceed via one of the larger Yukawa couplings, while its decay is typically controlled by smaller couplings because only light states are kinematically accessible. In particular, flavour changing transitions induced at the loop level are typically very relevant (see, e.g., [101]) and lead to the production via rare meson decays such as B → KS and K → πS, which are strongly constrained by a variety of experiments [15,[102][103][104]. Constraints on light scalars as well as projected sensitivities have been evaluated frequently in the literature, with a recent compendium of limits shown e.g. in ref. [105] (see also ref. [106], pointing out that such a light scalar could even drive cosmological inflation). In addition invisible decays of the SM Higgs into DM, h →χχ can give relevant constraints on the same product of couplings, g χ · sin θ, that is relevant for direct detection. In the following we briefly summarise the limits that we use in our analysis.

Invisible Higgs decay and signal strength
Data on Higgs bosons created at the LHC in principle constrain the SM Higgs mixing angle θ in two ways. First, invisible Higgs decays are constrained as BR(h → inv.) < 0.19 (95% C.L.) [107]. Second, the observed Higgs signal strength where σ h is the Higgs boson production cross section and N exp,SM h is the number of observed and expected Higgs bosons, respectively, is constrained to be µ > 0.89 (95% C.L.) [108]. In -10 -JHEP03(2020)118 our model the latter constraint is more stringent because the Higgs boson production cross section can only be reduced compared to the SM case, thus implying BR(h → inv.) < 0.11.
Specifically there are three effects that lead to a reduction of the signal strength: 1. Reduction of production cross section and decay widths for h, due to mixing.
3. Decays into two scalars, which depletes the branching ratio in the remaining channels.
In our case the ratio of the production cross sections is simply given by and for the branching ratios we have Here Γ 0 ≈ 4.1 MeV is the total SM Higgs width (without mixing), is the partial decay width for invisible decays and Γ SS is the Higgs boson decay width into two scalars (see eq. (2.5) and related discussion). Here we conservatively assume λ hs to be negligibly small, and hence set Γ SS ≈ 0. Taken together, the limit resulting from the predicted Higgs signal strength is thus given by which for m χ m h implies This limit will soon be improved with data from the 13 TeV run (see e.g. [109]). For the high luminosity phase of the LHC the direct bound on the invisible branching ratio will become more constraining and we use as an estimate of the future bound [110], thus strengthening the bound in eq. (4.6) by a factor of about 0.21. 5 JHEP03(2020)118

Limits on sin θ from beam dumps and colliders
As already mentioned, singlet scalars S can be efficiently produced via the decay of heavy mesons which in turn are copiously produced at the collision energies available at SPS and LHC based intensity frontier experiments, see e.g. ref. [101]. 6 For scalars too heavy to be produced in B meson decays this production mechanism is not available and the most constraining limit comes from LEP [105].
For scalar masses kinematically accessible to B → KS but too large to allow for K → πS, the strongest constraints are typically obtained from experiments where both the scalars can be created and their decay products can be detected. For experiments of this type the number of detected events thus scales with sin 4 θ (at the lower bound of sensitivity). For our analysis we use the results from LHCb [112] as well as a reinterpretation [105] of the past beam-dump experiment CHARM [113]. Values of sin θ just below the current sensitivity will be tested by LHCb in the high luminosity phase and we estimate the corresponding sensitivity (see appendix A for details). For smaller values of sin θ the SHiP experiment [15,17,114] and MATHUSLA [19,21] have the best sensitivities m S 5 GeV [6]. Both experiments aim at working in the background-free regime (see refs. [16,17] for detailed simulations for SHiP and refs. [20,21,115,116] for MATHUSLA). In reality, it is very difficult however to completely exclude the possibility of residual background events to take place in the detector. In case of SHiP such events can be distinguished from the signal events by making use of the spectrometer, mass reconstruction and particle identification. These options are not available in the case of MATHUSLA, implying that it is not straight-forward to compare its reported formal sensitivity (based on 2.3 events in the detector) to the one from SHiP. For this work we will therefore concentrate on the expected bounds from SHiP.
For smaller scalar masses, m S < m K −m π ≈ 350 MeV, experiments that search for rare kaon decays are typically more sensitive. The reason is that, unlike for particles heavier than kaons, one can perform precision measurements of the final state pion energy, on an event by event basis. The number of confirmed signal events thus no longer depends on the detection of the scalar decay products and therefore scales as sin 2 θ, i.e. is much less suppressed than in the case of heavier scalars. Both experiments E949 [117] and NA62 [7] search for rare K + → π + + MET decays and give bounds on the scalar mixing through the process K + → π + S. As this is independent of the decay of S, scalars with arbitrarily small masses can be constrained. In addition to the limit from E949 we estimate the sensitivity of NA62 during LHC Run 3 (see appendix B for details); we note that the resulting sensitivity is very similar to the sensitivity of the proposed KLEVER experiment shown in ref. [6]. used for this study. 7 For comparison, we also show existing limits from astrophysics; those will be discussed in more detail in section 5.4. While we mainly concentrate on the case m S < m χ as discussed above, we will also consider parameter regions in which m S > 2m χ and therefore invisible decays of the scalar naturally occur. In this case not all collider limits shown in figure 3 directly apply. To be specific, for m S = 0.1 GeV we will use the limits from E949 and NA62 as shown in figure 3 while for m S = 1 GeV the most stringent bound comes from the BaBar measurement of BR(B + → K +ν ν) < 1.6 · 10 −5 [118]. Making use of the partial decay width B → KS (see e.g. [63]), this translates into sin θ 6 · 10 −3 for m S 4 GeV.

Cosmological evolution of the dark sector
In this section, we describe the full thermal evolution of the dark sector particles, χ and S, which can be qualitatively divided into five, partially overlapping stages.

T > T dec
At high temperatures, the dark and the visible sector can be in chemical equilibrium due to the processes χχ ↔ ff , S ↔ ff and SS ↔ ff . In that case both sectors also share the same temperature, through efficient scattering of the involved particles, so 7 Nominally, there is a small gap in projected sensitivity at around mS ≈ 1 GeV and sin θ ≈ 5 · 10 −5 between the future exclusion power of the HL LHCb and the upper range of validity of the ShiP limits. This window however, will most likely be closed by i) slightly stronger (upper) limits of FASER2 [6] compared to SHiP and ii) the fact that in addition to B + → K + µ + µ − the channel B + → K + ππ will also be analysed by LHCb. The corresponding limit is expected to be more stringent than our estimate around mS ∼ 1 GeV owing to the fact that the branching ratio of S into pions is strongly enhanced compared to the branching into muons in this mass range, see e.g. [105]. When presenting our final results for the projected sensitivity of future experiments in this mass range, we will thus just use the lower sensitivity bound of SHiP, sin θ ∼ 10 −6 . is simply unity. For very small values of the mixing angle θ, however, the total interaction rate Γ DS↔SM between the two sectors is never large enough to bring them into thermal contact. Adding additional high-scale interactions to our model Lagrangian (2.1), on the other hand, would still allow to achieve chemical equilibrium for very large temperatures, without affecting the low-energy phenomenology. In this work, we will consider both of these possibilities and separately indicate the relevant parts of parameter space in our results.
T < T dec . At some temperature T dec the dark sector decouples from the visible sector.
To an approximation sufficient for our purpose, this happens when the total interaction rate equals the Hubble rate, This relation, in other words, allows to determine T dec as a function of our model parameters (m χ , m S , g χ and θ). In practice we compute Γ DS↔SM only for the three numberchanging reactions that keep up chemical equilibrium as mentioned in the previous paragraph (T > T dec ); elastic scattering Sf ↔ Sf will enforce kinetic equilibrium to be maintained slightly longer -but this is a small effect given that both scattering partners are relativistic around T dec . After decoupling the scalar mediators still retain a thermal distribution with temperature T S as long as they are relativistic (while non-relativistic scalars start to build up a chemical potential, even if a large quartic coupling λ S can keep them in local thermal equilibrium). Moreover, for sufficiently large dark couplings g χ , they are also kept in thermal equilibrium with the DM particles. Taken together, this leads to separate entropy conservation in the dark and visible sectors, and hence a temperature ratio that changes with the respective number of effective entropy degrees of freedom as It is worth stressing that this equation only holds as long as i) the scalars are still relativistic and ii) entropy is actually conserved, i.e. before S has decayed.
T > T fo . Above a certain temperature T fo , the DM particles will typically be in chemical equilibrium with the mediators and/or the SM heat bath. The former is achieved via the annihilation process χχ ↔ SS, while the latter is only relevant for T > T dec and happens dominantly through the s-channel process χχ ↔ ff . As the temperature approaches T fo , the DM number density freezes out and thereby sets the relic abundance of χ in the usual way. We numerically calculate the relic abundance with DarkSUSY [99] including the Sommerfeld enhancement of the annihilation rate, by modifying the implemented dark sector mediator model (vdSIDM) such as to fully include the temperature evolution of -14 -JHEP03(2020)118 ξ(T ) as specified in eq. (5.3); while the t/u-channel annihilation processes are identical to that model, we update the s-channel annihilation rate to the one relevant for our case, where Γ S represents the total width of S, and Γ S ( √ s) the width to SM particles assuming that S had a mass of √ s rather than m S . Assuming that DM is entirely produced via thermal freeze-out thus essentially fixes the dark coupling g χ as a function of the other model parameters. For our final results we will both use this assumption and demonstrate how the constraints on the model are affected if this assumption is relaxed (thus allowing for alternative DM production mechanisms). If m S m χ only s-channel annihilation is kinematically accessible; obtaining the correct DM abundance via thermal freeze-out would then require larger values of sin θ than compatible with the constraints presented in section 4.2 [45]. For 0.1 m S /m χ 1, on the other hand, the freeze-out process would involve two species that are no longer in chemical equilibrium, and where eq. (5.3) no longer applies. As such a situation would require a dedicated analysis, we will in the following leave this small part of the parameter space unexplored.
T < T fo . After the freeze-out of the dark matter particle, the mediator simply acts as an additional contribution to the energy density -until it decays to SM particles. Both stages have an impact on BBN, as discussed below. The corresponding lifetime of the scalar depends on the available decay channels and we adopt the results from ref. [119] in the following.
Let us, finally, stress again that, depending on the values of masses and couplings, DM freeze-out can in principle happen both before (T fo > T dec ) and after (T fo < T dec ) the decoupling of the two sectors. In this work we will restrict ourselves to light mediators with m S ≤ 0.1 m χ when discussing thermal DM production (see discussion above). In this case, taking into account the constraints on sin θ that result from direct searches for S, it turns out that we are always in the domain of T fo < T dec . Ultimately, this is a consequence of the Yukawa structure of the dark sector coupling, and the fact that we restrict our analysis to light DM.

Big Bang nucleosynthesis
In this section, we calculate BBN constraints for the scalar portal model using the formalism that was developed in [120][121][122][123], carefully taking into account the cosmological evolution of the dark sector as described in section 5.1. Specifically, in order to use the modelindependent constraints from ref. [120], we have to evaluate the values of τ S , T fo and ξ(T fo ) for every combination of m S , m χ and sin θ, thereby fixing g χ by the requirement of reproducing the observed relic abundance as described above. We then confront the predicted abundances of light nuclei in each parameter point to the following set of observed (Here Y p denotes as usual the primordial mass fraction of 4 He; we find that the requirement of obtaining the correct nuclear abundance ratio for 3 He/D leads to weaker limits than the above two constraints for all parts of parameter space.) Theoretical uncertainties associated to the nuclear rates entering our calculation are taken into account by running AlterBBN v1.4 [125,126] in three different modes corresponding to high, low and central values for the relevant rates. We then derive 95% C.L. bounds by combining the observational and theoretical errors as described in more detail in refs. [120,121].
In figure 4 we present the constraints from BBN as a function of m S and sin θ for fixed mass ratios m χ /m S (panels on the left) and as a function of m χ and sin θ for fixed mediator masses m S (panels on the right). The solid black line indicates the overall limit. In addition we also show which nuclear abundance causes an exclusion in which part of parameter space. In the pink (blue) region the 4 He abundance is too large (too small) compared to the observationally inferred values, while the constraints due to an underand overproduction of deuterium are shown in grey and purple, respectively. In addition, we show the lifetime of S as a function of sin θ for reference (green dashed lines); the fact that we identify excluded regions with τ S < 1 s implies that the often adopted 'conservative BBN limit' of τ S = 1 s is not always conservative.
It can be seen in the figure that the limits depend significantly on the value of m S as this quantity determines the possible decay channels of the mediator and therefore the lifetime, while the dependence on the dark matter mass is more indirect (but still not negligible) via its impact on the temperature ratio ξ. More concretely, we can distinguish the following regimes: • For m S > 2m µ the limit on sin θ is rather weak due to the small mediator lifetime above the muon threshold. For m S > 2m π also hadronic decays would become relevant for very small values of sin θ, see e.g. ref. [119]. 8 Overall we find that for values of sin θ relevant to this study, the mediator already decays before the onset of BBN for m S > 2m µ , therefore not causing any observable consequences for the range of direct detection cross sections we consider.
• For 2m e < m S < 2m µ the scalar can decay before, during or after BBN, depending on the value of sin θ. In this region of parameter space, the presence of the dark sector influences BBN via two different effects: (i) an increase of the Hubble rate due to the extra energy density of the dark sector and (ii) entropy injection into the SM heat bath due to scalar decays into electromagnetic ration. Both affect the synthesis of light elements as discussed in detail in ref. [120] and are fully taken into account in our evaluation. For lifetimes τ S 10 4 s photodisintegration also becomes relevant,

JHEP03(2020)118
but does not exclude any additional regions of parameter space (at least in case DM is produced via freeze-out).
• For m S < 2m e , the scalar S can decay only into photons, which leads to a drastically increased lifetime. Consequently, for comparably small values of sin θ, the mediator outlives the creation of the light elements, thereby acting as an additional relativistic degree of freedom, whose presence can be robustly excluded by current BBN constraints (even stronger constraints in the case of such late decays arise from photodisintegration of light elements [120,127,128] and the CMB [129]). For very large values of sin θ, on the other hand, S again decays during BBN -but since this case is strongly excluded by other considerations, cf. section 4, we do not perform a detailed study of BBN limits in this regime.
As discussed in section 5.1, for sufficiently small values of the mixing angle θ, the two sectors will never thermalise via the Higgs portal; this is indicated by the dashed black line in figure 4. The conservative BBN limits which do not make any additional assumptions about early universe cosmology therefore naïvely end at this line. 9 While the overall limit looks very similar for a given scalar mass m S , the thermalisation line is rather sensitive to the DM mass m χ which directly translates into rather different conservative BBN limits for the different cases. For the calculation of the BBN limits below this line we assume ξ(T → ∞) = 1, i.e. that both sectors were thermalised at very large temperatures by some additional processes that are not covered by the model Lagrangian in eq. (2.1). If the two sector never thermalised, on the other hand, the bound would depend on the initial temperature ratio ξ (or ratio of energy densities). For ξ(T → ∞) = 0 only the freeze-in contribution would remain. Nevertheless, even in this case stringent bounds from photodisintegration and the CMB are expected for sizeable regions of parameter space. To indicate this additional uncertainty the BBN limits below the thermalisation line have a lighter shading.

Dark matter self-interactions
The exchange of a scalar particle generates an attractive Yukawa potential between two DM particles, resulting in a self-interaction rate that strongly depends on the couplings g χ , the DM and mediator masses m χ and m S , and the relative velocity v of the scattering DM particles. For the range of parameters we are interested in here, in particular, these interactions typically show a characteristic resonant structure, resulting in a large enhancement or suppression of the momentum transfer cross section σ T when varying, e.g., the dark coupling g χ . In this regime, analytic expressions are available that result from restricting the analysis to s-wave scattering and approximating the Yukawa potential by a Hulthén potential [43]. While these expressions result in a reasonable estimate for height 9 The bounds may be significantly stronger taking into account an irreducible contribution from freezein production of either the DM particle or the mediator. In fact even a small abundance of mediators is constrained if the lifetime is such that photodisintegration is relevant. A detailed exploration of this parameter region is left for future work. Also, as stated before, we assume λ hs 0. Sizeable values for this quartic inter-sector coupling could further shift the region of thermalisation to smaller values of sin θ.  and as a function of m χ for fixed masses m S (right panels). Below the dashed black line, the Higgs portal by itself is insufficient to ever thermalise the dark sector with the SM. In addition to the overall limit (solid black line), we also separately show the regions of parameter space which are excluded due to D underproduction (grey), D overproduction (purple) and/or 4 He underproduction (blue), 4 He overproduction (pink). For m S 0.1m χ the thermal evolution is not fully captured by our calculation and thus the limits are only approximate, as indicated by the hatch pattern.

JHEP03(2020)118
and location of resonances in σ T , we find that they significantly underestimate the numerical value of σ T in the vicinity of anti-resonances (see also [130]). In our analysis, we thus always solve the underlying Schrödinger equation for the full Yukawa potential numerically, including also higher partial waves. We do so by following the treatment in ref. [131], thus also correctly taking into account the indistinguishability of DM particles in the definition of the momentum transfer cross section.
In the cosmological concordance model, DM is successfully described as a collision-less fluid. In fact, astrophysical observations from dwarf galaxy to cluster scales stringently limit how much the properties of the putative DM particles can deviate from this simple picture (for a review, see ref. [132]). Here we adopt as a fiducial maximal value for the allowed effective momentum transfer cross section at a relative DM velocity of v rel χ = 1000 km/s. This corresponds roughly to the constraint inferred from the observation of colliding galaxy clusters [132][133][134], which are highly DMdominated systems and hence often argued to be less prone to modelling uncertainties of the baryonic component. Another advantage is that these observations more directly constrain the DM self-interaction rate, while the widely used translation of individual halo properties -like their core size -to bounds on σ T is subject to a non-negligible intrinsic modelling uncertainty (for a recent discussion, see ref. [135]). 10 We note that astrophysical observables do not depend on σ T alone, but may also depend on the frequency of scatterings [139]. For the case of a light mediator as studied here, this can be substantially different from a contact-like (isotropic) DM scattering, which is usually assumed in N -body simulations (see however refs. [140][141][142] for first studies including angular dependent and frequent scatterings). On the other hand, bounds on the self-interaction rate have been reported that are significantly stronger than the fiducial maximal value(s) of σ T that we adopt in our analysis [138,[143][144][145]. Overall, taking into account the above mentioned uncertainties, we expect that eq. (5.7) leads to realistic constraints on the dark coupling g χ .
Due to the characteristic resonant structure of σ T , the inversion of these limits to constraints on g χ is in general not unique. For given values of m χ and m S , in particular, there is always a maximal value g min χ such that eq. (5.7) is satisfied for all values of g χ < g min χ . In the resonant regime, however, it is possible to hit anti-resonances, and thus to decrease the cross section by increasing the coupling beyond g min χ . In other words, there can be further -sometimes only very narrow -parameter ranges with g χ ∈ (g min χ , g max χ ) that satisfy eq. (5.7). Here, g max χ denotes the maximal value for which the self-interaction constraint can in principle be satisfied, due to the appearance of anti-resonances in the scattering amplitude. Requiring g χ < g max χ is thus the most conservative way of implementing the self-interaction constraint, but it neglects the fact that many values of g χ < g max χ are ac- 10 Observations of dwarf galaxies, and their translation to limits on σT , are prone to much larger uncertainties [136][137][138]. On the other hand, for light mediators the effective self-interaction rate is considerably stronger in these systems than in galaxy clusters. Adopting for example σT /mχ < 100 cm 2 /g for relative DM velocities of v rel χ = 30 km/s, as a relatively conservative fiducial constraint, would lead to stronger constraints than eq. (5.7) only for mS 10 −3 mχ.
-19 -JHEP03(2020)118 tually excluded; requiring g χ < g min χ is more aggressive, but in some sense more generic (because it applies even outside the range of couplings where anti-resonances appear). To reflect this situation, we will in the following show results for both sets of constraints independently. We note that numerically it is straight-forward to determine g min χ and g max χ because the anti-resonances are much less pronounced than what the analytic approximation for s-wave scattering [43] would suggest.

Further astrophysical and cosmological bounds
Weakly coupled light particles can be copiously produced in the interior of stars or in the hot core of a supernova (SN) via their interactions with electrons or nucleons. For sufficiently weak couplings these particles escape the celestial body without further interactions and therefore constitute a new energy loss mechanism which is constrained by observations. We take the resulting limits on sin θ from red giants (RG) and horizontal branch stars (HB) from ref. [146]. The bound from SN 1987A which extends to larger masses because of the correspondingly larger core temperature we take from ref. [105], noting that this is an order of magnitude estimate with inherently large uncertainties. The lower boundary of this limit is determined by estimating the additional energy loss due to escaping scalars which would shorten the observed neutrino pulse. For m S < 2m χ the SN limit does not extend to arbitrarily large couplings due to efficient trapping of light scalars inside the SN; for larger scalar masses, on the other hand, there is no upper boundary of the limits because the scalar will decay invisibly to DM particles that escape the SN without interacting. These constraints are shown as grey areas in the sin θ − m S plane in figure 3, where we have used a hatched filling style for the SN bounds (assuming m S < 2m χ ) to stress their intrinsic uncertainty.
We remind that the usually very strong CMB bounds on annihilating light DM can be evaded in this model because the annihilation proceeds via a p-wave and is hence strongly velocity-suppressed (also, there are no remaining light degrees of freedom that would change the expansion rate at these times). While the elastic scattering of DM with SM particles is enhanced for light mediators S, the coupling of S to photons is still not sufficient to prolong kinetic decoupling until times where an appreciable cutoff in the matter power spectrum, and hence a potential imprint in e.g. Ly-α data, would be expected (see ref. [147] for a more detailed discussion).

Results
As motivated in the introduction, we now want to combine the various constraints discussed in the previous sections in the 'direct detection plane', i.e. as bounds on σ SI as a function of the DM mass. For a given value of the scalar mass m S (and fixed m χ ) only the invisible Higgs decay constrains the same combination of parameters (sin θ · g χ ) that enters the expression for σ SI , cf. eq. (3.5). In all other cases, we thus need to combine two types of observations to constrain these parameters individually. As the particle physics constraints discussed in section 4, but also the astrophysical constraints from section 5.4,

JHEP03(2020)118
are essentially insensitive to g χ , a handle on the dark coupling has to be provided by cosmology. Concretely, we will distinguish three versions of cosmological constraints (roughly ordered by decreasingly strong underlying assumptions): • Cosmo 1 ('thermal production'). The present dark matter abundance can be fully explained by the production of χ particles via freeze-out in the early universe, as detailed in section 5.1, which requires dark and visible sector to have been in thermal equilibrium at some point. No further interactions than those specified in eq. (2.1) are assumed. 11 • Cosmo 2 ('generic self-interactions'). No connection between the dark coupling g χ and the DM abundance is assumed, allowing for other DM production mechanisms. 'Generic' constraints on DM self-interactions are adopted, as detailed in section 5.3, i.e. we assume that g χ does not lie close to an anti-resonance in the elastic scattering cross section.
• Cosmo 3 ('conservative self-interactions'). As Cosmo 2, but implementing the most conservative constraints from DM self-interactions; larger values of g χ would thus violate eq. (5.7) even when finely tuned to lie on an anti-resonance.
BBN constraints are tightly coupled to the assumed thermal history, so for those we will always assume thermally produced DM ('Cosmo 1'). Finally, we always adopt a hard 'perturbative unitarity limit' of g χ < √ 4π in case the respective cosmological constraint would be weaker.
In figure 5 we show our results for selected mass ratios m S /m χ < 2 where invisible decays of the scalar are kinematically forbidden. In each case, the left column displays current limits while the right column displays projected limits. Conventional direct detection limits (rescaled from figure 1) are shown as grey areas. Current limits from cosmic-ray upscattering (figure 2) are shown in light blue; for the projected limits we take the sensitivity of DARWIN [79], based on the assumption that the recoil threshold can be lowered to 1 keV. Limits from invisible Higgs decay, cf. eqs. (4.5) and (4.7), are shown in green. In red, we combine the particle physics limits shown in figure 3, while the yellow lines show a combination of the astrophysical limits discussed in section 5.4. 12 The various ways of implementing cosmological limits are indicated with dotted lines ('Cosmo 1'), dashed lines ('Cosmo 2') and solid lines ('Cosmo 3'), respectively. Given the difficulties in accurately computing the thermal evolution of the dark sector for m S 0.1 m χ , see the discussion in 11 For mS mχ the relic density actually depends on the same combination of couplings (sin θ · gχ) as direct detection rates; as noted before, it is not possible to obtain the correct relic density and at the same time satisfy existing bounds on θ in this case. 12 From the shape of these limits, the potentially controversial part that derives from SN bounds is clearly discernible. We note that only in the case of 'Cosmo 1', which fixes gχ by the requirement of obtaining the correct relic density, the range in θ excluded in figure 3 (for a given value of mS) translates to a correspondingly excluded range of σSI. If instead there is only an upper limit on gχ, as in the case of 'Cosmo 2' and 'Cosmo 3', the direct detection cross section σSI ∝ g 2 χ sin 2 θ remains essentially unconstrained by this bound (in other words, while gχ still cannot be chosen so small that sin θ > 1, for a given value of σSI, this only results in a limit too weak to be visible in the figure).
-21 -JHEP03(2020)118 section 5.1, we do not display 'Cosmo 1' limits in this regime. For the case of BBN limits (shaded blue areas), we also indicate (as in figure 4) the parameter region where additional high-scale interactions would be required to thermalise the dark and visible sectors in the very early universe; BBN limits that rest on this additional assumption are plotted with a hatched filling style. (Note that, compared to figure 4, BBN limits appear to have a stronger dependence on m χ here; this is exclusively because the relic density constraint fixes g χ as a function of m χ .) The figure nicely illustrates the complementarity of the different approaches to test models with light mediators. If DM is thermally produced, in particular, current bounds already reduce the remaining parameter space for sub-GeV DM to a relatively small region of mediator masses above a few MeV, and mass ratios 0.01 m S /m χ 0.1 (see also [148] for a discussion of BBN limits in a similar context). Here it is worth commenting that BBN limits far below the thermalisation line essentially just constrain the assumed highscale temperature ratio between the two sectors; in this sense they simply exclude this additional assumption and are completely independent of the specific model Lagrangian stated in eq. (2.1). On the other hand the robust bounds may significantly extend below the thermalisation line taking into account the irreducible contribution from freeze-in production.
Even if no assumptions about the thermal history and production of DM is made, on the other hand, the combination of the limits displayed in figure 3 with those stemming from DM self-interactions translate to highly competitive constraints on σ SI . For mediator masses close to the DM mass, in particular, those constraints already fully cover the expected reach of upcoming direct detection experiments. Interestingly, we arrive at this conclusion independently of which set of SIDM constraints we implement ('Cosmo 2' or 'Cosmo 3'); let us stress, however, that the limits presented in figure 5 indeed strongly depend on fully solving the Schrödinger equation to obtain the self-interaction cross section σ T in the resonant regime, rather than following standard practice and simply adopting analytic results for s-wave scattering. From the perspective of future direct detection experiments, the most interesting parameter range to be probed -fully orthogonal to what can be tested by particle physics experiments -is the sub-GeV DM range combined with scalar masses significantly lighter than DM (but heavier than about 0.2 MeV, where astrophysical limits start to dominate).
In order to add a slightly different angle to the above discussion, we show in figure 6 the same constraints for selected fixed scalar masses m S instead. This includes kinematical situations with m S > 2m χ where the scalar can decay very efficiently to two DM particles, i.e. through an invisible channel. As discussed in section 4.2, the particle physics constraints hence need to be adapted correspondingly, and we thus only keep those limits shown in figure 3 that are still relevant in this situation (and add that from BaBar [118] for the case of m S = 1 GeV). 13 For invisible decays, furthermore, there is also no upper boundary to the area excluded by energy loss arguments in supernovae (as in figure 3). This implies that for JHEP03(2020)118  Figure 5. Current (left) and projected (right) limits on the elastic scattering cross section with nucleons in the zero-momentum transfer limit, for fixed scalar to DM mass ratios m S /m χ that do not allow invisible decays of S. For astrophysical and particle physics limits combined with cosmological limits, dotted lines assume thermal DM production via freeze-out ('Cosmo 1'), dashed lines instead implement generic DM self-interaction constraints ('Cosmo 2') while solid lines result from tuning g χ such as to resonantly suppress the DM self-scattering rate ('Cosmo 3'). See text for further details. 2m χ < m S 0.2 GeV, unlike the situation in figure 5 for visible decays, the combination of SN bounds and SIDM constraints indeed does combine to a very competitive bound on σ SI (though, as discussed in section 5.4, SN limits should be interpreted with care).
While the limits from cosmic-ray upscattered DM now become more visible, it is clear that they are never competitive to other limits in this model. Also limits from invisible Higgs decay, while significantly more stringent, are rarely strong enough to be competitive; this would change only with a dedicated Higgs factory like the ILC. In general, one can say that astrophysical, particle physics and direct detection limits probe the parameter space from rather orthogonal directions. While astrophysical constraints are most relevant for small DM (or, rather, mediator) masses, direct detection experiments place the strongest limits for large DM masses. The m χ -dependence of constraints on σ SI stemming from particle physics, on the other hand, is somewhat weaker. Consequently, particle physics (combined with cosmological input) tends to place the most relevant constraints on the model at intermediate DM masses (for the sub-GeV range that we consider here), and the most promising avenue for direct DM searches appears to lie in lowering the detection threshold, even slightly, in a way that compromises the overall sensitivity as little as possible (this, in other words, will test more of the so far unprobed parameter space than focussing on very low thresholds at the expense of overall sensitivity).

Discussion and conclusions
In this work we have considered the prospects of future direct detection experiments to test uncharted parameter space for light (sub-GeV) dark matter. It is natural in this context to expect additional light particles mediating the interactions between dark matter and the target nuclei in order to achieve a sufficiently large scattering cross section. To alleviate the strong cosmological bounds from the CMB we have concentrated on a scenario in which dark matter couples via a scalar mediator (with coupling g χ ) such that dark matter annihilations proceed via p-wave and are therefore strongly suppressed at the time of the CMB. This allows the dark matter relic abundance to be set via thermal freeze-out, although other production mechanisms are possible, and our bounds also apply to more general cases. We assume that couplings to Standard Model states are induced by the well-known Higgs portal with mixing angle θ.
The DM scattering cross section off nuclei is then proportional to the product of couplings, g 2 χ · sin 2 θ. To map out the available parameter space we evaluate and compile the relevant limits both on sin θ and on g χ from current and near future particle physics experiments, BBN, astrophysics and cosmological considerations. We also show limits on light DM particles upscattered by cosmic rays, which turn out never to be competitive in the model considered here. In our analysis we paid special attention to cosmological constraints which, while requiring certain assumptions, cannot be avoided altogether in a given model. Indeed, they provide quite in general the missing link to translate a variety of constraints on portal models to constraints on the scattering cross section relevant for direct dark matter searches.

JHEP03(2020)118
In our analysis we update and carefully extend previous recent work similar in spirit [45,46,48,49]. Concretely, we re-scale direct detection limits according to the full Q 2 dependence in section 3, and present a genuinely new treatment of cosmic-ray up-scattered DM for such a case -which, as we stress, is applicable also to other scenarios that invoke Q 2 -dependent scattering. We present updated invisible Higgs decay constraints on this specific model (in section 4.1), and add genuinely new estimates of the future LHCb and NA62 sensitivities (see appendix) to our compilation of up-to-date particle physics constraints and projected sensitivities in section 4.2. A significant refinement compared to what is typically done in the literature is further that we consistently implement the cosmological evolution of this scenario in full detail (section 5.1), which we use both for precision calculations of the relic density and a careful evaluation of BBN constraints (going well beyond the standard procedure of simply translating model-independent constraints to the model at hand). We also point out that the self-interaction cross section (section 5.3) has to take into account the (in)distinguishability of external particles and needs to be evaluated beyond the s-wave approximation to avoid the appearance of artificially deep anti-resonances and hence overly weak constraints; correcting for this, we instead find that DM self-interactions generically lead to limits comparable to those resulting from the assumption of DM thermally produced via the freeze-out mechanism.
Overall we find that, almost independently of the dark matter production mechanism, strong bounds on the maximally possible nuclear scattering rate exist for large regions of parameter space. Nevertheless, some regions remain safe from the combination of existing (or expected) constraints from accelerators, astrophysics and cosmology, motivating the development and construction of future direct detection experiments which could explore these regions. This not only requires low thresholds for the recoil energies, but at the same time sensitivities better than what is presently achievable at dark matter masses around 1 GeV.  Figure 7. Lifetime of the scalar particle S as a function of its mass, thereby fixing sin θ to the lower bound of the current LHCb sensitivity as shown in figure 3 (blue line). The black dashed lines correspond to lifetimes of 1 ps and 10 ps, thus indicating the borders between the prompt region (τ S < 1 ps), the intermediate region (1 ps < τ S < 10 ps) and the large displacement region (τ S > 10 ps). See text for further details.
for a dataset with collision energy √ s = 7 and 8 TeV and integrated luminosity L 0 = 3 fb −1 . For this analysis, the parameter space of the scalar was divided into (i) a prompt region with the lifetime of the scalar τ S < 1 ps, (ii) an intermediate region with 1 pc < τ S < 10 pc and (iii) a large displacement region for τ S > 10 pc. Background events were expected in the first two regions, while the last region was background free. In figure 7 we show τ S , fixing sin θ to the lower bound of the current sensitivity of the LHCb experiment. We conclude that no background is expected for m S < 3.7 GeV (which is the region of interest to us), while for higher masses we need to consider a non-zero background contribution.
To estimate the sensitivity of a similar analysis in the high-luminosity (HL) phase of the LHC, we assume the total integrated luminosity of LHCb to be L HL = 300 fb −1 and the centre-of-mass energy to be √ s = 13 TeV. The corresponding increase of the number of produced B mesons in the direction of LHCb can be estimated as R = L HL · σ 13 (pp → B + + X) L 0 · σ 8 (pp → B + + X) ≈ 162.2 , (A.1) where σ 13/8 (pp → B + + X) is the production cross section of B + mesons which fly into direction of the LHCb detector for energies 13 and 8 TeV respectively. We estimated these cross sections using FONNL (Fixed Order + Next-to-Leading Logarithms) -a model for calculating the single inclusive heavy quark production cross section, see [149][150][151][152] for details. For the region in which background events are expected, we assume for simplicity that the number of background events also increases by the factor R. We estimate the future sensitivity as For the case of large displacements, while no background events are expected, we need to take into account the probability of the scalar to decay inside the region where displaced -27 -JHEP03(2020)118 vertices can be observed, l min ≤ l decay ≤ l max with l min = 3 mm and l max = 0.6 m, [112]. This probability can be written as P decay (θ) = e −l min /l decay (θ) − e −lmax/l decay (θ) , (A. 3) where l decay = cγ S τ S is decay length of the scalar in the lab frame and γ S is the corresponding Lorentz factor. We estimate the average Lorentz factor of the scalar to S be (see appendix C in [114]) where γ S,rest is the Lorentz factor of S in the rest frame of the B-meson. The average energy of the B-mesons in the direction of LHCb we take from FONNL, E B ≈ 80 GeV for both centre-of-mass energies. Taking everything together, we estimate the future sensitivity in this regime as θ 2 future (m S )P decay (θ future (m S )) = 1 R θ 2 current (m S )P decay (θ current (m S )) .

B Estimate of future NA62 sensitivity
In this appendix we briefly describe how we estimate the sensitivity of the NA62 experiment with respect to light scalars produced in K + → π + S. 14 One of the main physics goals of NA62 is the measurement of the rare decay K + → π + νν, allowing for a direct determination of the V td CKM matrix element [7]. The observed final state is a π + plus missing momentum. If the scalar S is sufficiently long-lived to decay outside of the detector it would contribute to the same final state and can therefore be constrained with this search mode. A crucial difference between the expected signal from K + → π + S compared to the SM process K + → π + νν is the distribution of the 'invisible mass', which in the case of decays into S peaks at the mass m S while for the SM process (as well as other SM backgrounds which contribute to this final state) the distribution is rather flat, 15 see e.g. [155]. The number of kaons expected during LHC Run 3 [6] is estimated to be N K 10 13 , which (scaling up the results from [155]) corresponds to about 35 events in the signal region with a rather flat distribution in the missing mass.
To compare this to the expected signal from K + → π + S we have to take into account the corresponding branching ratio as well as the total selection efficiency for this process, which in general will depend on m S . The expected number of events is then N obs S = N K · BR(K + → π + S) · . (B.1) For our analysis we approximate the total efficiency as = 0.3 [156]. The relevant branching ratio is given by (see e.g. [104]) BR(K + → π + S) 1.85 · 10 −3 sin 2 θ 1 − (m S + m π ) 2 m 2 14 For a sensitivity estimate of NA62 to light scalars with different coupling structure see e.g. [153]. 15 Due to large backgrounds from K + → π + π 0 the mass range 130 MeV < mS < 150 MeV should be excluded from the analysis [154].

JHEP03(2020)118
Taking into account that the experimental resolution of the missing mass is about 1/35 of the signal region, we expect about 1 background event from SM processes after all cuts. The 95% CL upper limit on the scalar would then correspond to ∼ 5 events, which is what has been required for our result shown in figure 3.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.