Looking for the WIMP Next Door

We comprehensively study experimental constraints and prospects for a class of minimal hidden sector dark matter (DM) models, highlighting how the cosmological history of these models informs the experimental signals. We study simple"secluded"models, where the DM freezes out into unstable dark mediator states, and consider the minimal cosmic history of this dark sector, where coupling of the dark mediator to the SM was sufficient to keep the two sectors in thermal equilibrium at early times. In the well-motivated case where the dark mediators couple to the Standard Model (SM) via renormalizable interactions, the requirement of thermal equilibrium provides a minimal, UV-insensitive, and predictive cosmology for hidden sector dark matter. We call DM that freezes out of a dark radiation bath in thermal equilibrium with the SM a WIMP next door, and demonstrate that the parameter space for such WIMPs next door is sharply defined, bounded, and in large part potentially accessible. This parameter space, and the corresponding signals, depend on the leading interaction between the SM and the dark mediator; we establish it for both Higgs and vector portal interactions. In particular, there is a cosmological lower bound on the portal coupling strength necessary to thermalize the two sectors in the early universe. We determine this thermalization floor as a function of equilibration temperature for the first time. We demonstrate that direct detection experiments are currently probing this cosmological lower bound in some regions of parameter space, while indirect detection signals and terrestrial searches for the mediator cut further into the viable parameter space. We present regions of interest for both direct detection and dark mediator searches, including motivated parameter space for the direct detection of sub-GeV DM.


Introduction
The existence of some form of dark matter (DM) constituting 26% of the present-day energy budget of our universe is well-established through its gravitational imprint on baryonic matter [1]. No evidence to date indicates that DM must interact in any way beyond gravitationally. The cosmological history of DM, however, will typically require DM to have some non-gravitational interaction(s) responsible for establishing its observed relic abundance, and these interactions can leave potentially observable footprints. For instance, the cosmic coincidence of the "weakly interacting massive particle (WIMP) miracle" implies that new stable weak-scale particles with weak interactions that freeze out of the thermal Standard Model (SM) plasma in the early universe can provide a good DM candidate. The cosmic abundance of WIMPs is directly determined by their coupling to the SM, and thus this class of models makes sharp predictions for signals accessible to a variety of experiments. While the parameter space for thermal WIMPs is now acutely limited by the interplay of null results at direct and indirect detection experiments and at the Large Hadron Collider (LHC) [2], thermal DM that freezes out directly to SM particles via new beyond-the-SM (BSM) mediator(s) similarly has a cosmological abundance directly set by the strength of its interactions with the SM, and has thus driven the terrestrial DM discovery program in recent years. These models, too, are becoming increasingly challenged by the lack of signals to date [3,4].
This class of thermal relics, however, represents only a fraction of possible identities for dark matter. Hidden sector freezeout (HSFO [5][6][7][8]), where the DM relic abundance is chiefly determined by interactions internal to a thermal dark sector with little to no involvement of the SM, provides a much broader class of models. In this paper, we survey the current constraints and future discovery prospects in the simplest exemplars of hidden sector freezeout. In these simple and minimal models, DM is a thermal relic that annihilates not to SM states, but to pairs of dark mediators that subsequently decay via small couplings into the SM. We take these small couplings to be the leading interaction between the HS and the SM, and consider the well-motivated and generic case where this interaction is renormalizable.
Any theory where DM arises from an internally thermalized dark sector must also address the question: how was this dark sector populated in the early universe? The most minimal cosmological history for a dark sector is for it to interact strongly enough with the SM that the two sectors were in thermal equilibrium at early times. In this case, the existence of a thermal SM plasma in the early universe guarantees the population of the dark sector. We call DM that freezes out from a thermal dark radiation bath in thermal equilibrium with the SM a WIMP next door. Mandating this cosmological history for the dark sector imposes a lower bound on the interactions between the dark sector and the SM today, the thermalization floor. The parameter space for WIMPs next door is bounded: the DM mass must lie between ∼ 1 MeV (to preserve the successful predictions of Big Bang Nucleosynthesis (BBN) and ∼ few TeV (from perturbative unitary), while the coupling between the SM and the HS must be sufficiently strong to thermalize the dark sector with the SM prior to DM freezeout.
The aim of this paper is to establish this bounded parameter space for two minimal models of HS freezeout and systematically map out how this parameter space can be tested in indirect detection, accelerator, and direct detection through a variety of experiments spanning the cosmic, energy, and intensity frontiers. The characteristic signatures of hidden sector freezeout are largely dictated by the Lorentz quantum numbers of the DM and the mediator, together with the choice of portal operator. We focus here on dark sectors which have a leading renormalizable coupling with the SM, through either the vector portal interaction, 2 cos θ Z Dµν B µν , (1.1) or the Higgs portal interaction, (1.2) We will use two simple reference models in this work, • HSFO-VP: fermionic DM χ, annihilating to vector mediators, Z D , that couple to the SM through the vector portal; and • HSFO-HP: fermionic DM χ, annihilating to scalar mediators s that couple to the SM through the Higgs portal.
These models of HSFO can be probed via complementary methods across different experimental frontiers. Direct searches for the dark mediator are the most sensitive test at accelerator-based experiments, far outpacing more traditional collider searches for DM that rely on a missing energy signature.
Direct detection experiments can access the cosmological lower bound on the portal coupling in significant portions of the parameter space. Indirect detection remains a powerful probe, provided the DM has an appreciable s-wave annihilation cross-section, as in our minimal vector portal model. Our minimal Higgs portal model, on the other hand, freezes out through p-wave interactions, placing traditional cosmic ray signals largely out of reach. The constraints on our simple reference models provide a reasonably general guide to the physics of more complicated hidden sectors, as we discuss below. We begin with a discussion of WIMPs next door in Sec. 2, where we establish the physical parameter space of our models. In Sec. 3 we discuss different experimental avenues to test this parameter space. In Secs. 4 and 5 we show the consequences for vector and Higgs portal models respectively, and in Sec. 6 we summarize our results. Three Appendices describe details of our calculations of thermal scattering rates, Sommerfeld enhancements, and bounds from dwarf galaxies.

Parameter Space for Minimal Hidden Sector Freezeout and the WIMP Next Door
In hidden sector freezeout, DM is part of a larger dark sector that is thermally populated in the early universe. As the universe expands and cools, the relic abundance of DM is determined by the freezeout of its annihilations to a dark mediator state, χχ → φφ, with little to no involvement of SM particles [5][6][7][8]. In the simplest realizations of hidden sector freezeout, these dark mediators, φ, are cosmologically unstable, decaying into the SM through a small coupling. These decays must occur sufficiently rapidly to avoid disrupting the successful predictions of BBN, thus generally requiring τ φ 1 s, and providing a cosmological lower bound on the strength of the coupling of the mediators to the SM. When the interaction that allows φ to decay to the SM is the leading interaction between the two sectors, it will additionally control the thermalization of the dark sector and the SM in the early universe.
Requiring the SM and the dark sector to be in thermal equilibrium prior to DM freezeout is the simplest and most minimal cosmology for the origin of the dark sector. DM freezing out from a thermal dark radiation bath in equilibrium with the SM radiation bath is what we define as a WIMP next door. We will focus here on the well-motivated cases where the vector or scalar portal operators (1.1-1.2) mediate the leading interactions between sectors, and establish the observable consequences. If the portal coupling is sufficiently large to ensure that the SM and the dark sector were in thermal equilibrium for some temperature above the DM freezeout temperature, T eq > T f , the existence of the SM thermal bath is then sufficient to guarantee the population of the dark sector. If, on the other hand, the portal interactions cannot thermalize the two sectors prior to DM freezeout, then some other mechanism, such as asymmetric reheating [9][10][11], must be invoked to populate the dark thermal bath in the early universe.
When the leading interaction between sectors is renormalizable, this minimal cosmology is additionally UV-insensitive: the scattering rates controlling thermalization obey Γ ∝ T , and become more important in comparison with H ∝ T 2 /M P l as the temperature drops. Thus min (T eq ), the minimum value of consistent with thermalization at a temperature T eq , does not depend on the unknown reheating temperature of the universe (provided T RH > T eq ) or other unknown UV particle content. This cosmic origin for DM also significantly sharpens predictivity by limiting the degree to which the temperature of the dark sector can differ from the temperature of the SM.
In order to determine the thermalization floor min (T eq ), we have to distinguish between two cases: first, when the hidden sector contains (at least) one relativistic species at T eq , and second, when all species in the hidden sector are already nonrelativistic at thermalization, T eq < m/2.46 for all masses. The value of T = m/2.46 is the point where a bosonic species contribution to g * drops by a factor of 2, and is our definition for when a species transitions from relativistic to non-relativistic. In the first case, the energy in the hidden sector radiation bath is the same per degree of freedom as in the SM, and thermalization requires that inter-sector reactions are efficient enough to transfer a sizable amount of energy per SM degree of freedom. In the second case, all hidden sector species have exponentially suppressed number densities at T eq , and the energy that must be transferred from the SM to thermally populate the hidden sector is thus exponentially reduced. The resulting bounds on minimal coupling strengths are correspondingly much weaker.
We focus here on the first case where the HS has a radiation bath at T eq . In this cosmology the lower bound of the thermalization floor is typically far more stringent than the lower bound from requiring mediator decays to occur prior to BBN. We require that the two sectors thermalize at least at T eq = T f . For simplicity, we consider minimal models that consist only of a dark matter species χ and a dark mediator φ. In order to have a dark radiation bath at DM freezeout, we thus require the mediator to have m φ < 2.46 T f . When T m φ , 2 ↔ 2 scatterings (SM)φ ↔ (SM)(SM) are the dominant process responsible for equilibrating the two sectors. When T ∼ m φ , 1 ↔ 2 scatterings φ ↔ (SM)(SM) become dominant. This temperature scaling is evident in Fig. 1, where we show 1 ↔ 2 and 2 ↔ 2 scattering rates in each of our models as a function of the temperature. In the absence of mass thresholds, Γ 2↔2 ∝ T at high temperatures, while Γ 1↔2 ∝ m 2 φ /T . The SM has many mass thresholds, which makes the temperature dependence of the net scattering rates less transparent. Full details of the calculation of these scattering rates are presented in Appendix A; as discussed there, the thermalization floor that we obtain is an initial estimate, computed up to a factor of ∼ 2. The resulting new cosmological lower bound on portal couplings is shown in Fig. 2 as a function of T f , in the regime where m φ 0.1 m χ . The thermalization floor is insensitive to the mediator mass as long as 2 ↔ 2 rates dominate the scattering, a condition that holds generically (but not always) when the mediator is relativistic at the time of freezeout.
Both our minimal models can be described by four independent parameters, namely the DM mass, the mediator mass, the portal coupling , and the coupling α D between DM and the mediator. Simplified model approaches can be effective at highlighting the key physical features of classes of   Figure 1. Rates for 2 ↔ 2 (in black) vs. 2 ↔ 1 transitions (in different shades of red depending on the mediator mass) as a function of the temperature. Dotted lines correspond to the 2 ↔ 1 rates after the mediator becomes non-relativistic, T < m/2.46, with m the mediator mass. The gray region at low temperatures corresponds to the uncertain regions near Λ QCD ≡ 300 MeV, above which α s becomes large, and the leading order calculation is unreliable, and below which 3 ↔ 2 pion processes (included in the 2 ↔ 2 rate, shown in light blue) dominate the equilibration. The gray region at high temperature is near the electroweak phase transition T c = 160 GeV; see Appendix A for more details. Left: HSFO-VP. In this model the 2 ↔ 2 scattering rate is nearly linear with temperature above the chiral phase transition. After pion processes become ineffective and QED processes dominate, the scaling is nearly linear again. Right: HSFO-SP. In this model the 2 ↔ 2 rate is more sensitive to mass thresholds. It drops sharply after the chiral phase transition (kaon processes have been neglected).
DM theories [3,[12][13][14][15], and, in that spirit, our simple HSFO models can be taken as useful guides to the physics of a general WIMP next door, as we discuss further below. We emphasize, however, that our minimal HSFO models are, themselves, UV-complete and self-consistent. WIMPs next door have a sharply defined and bounded parameter space. The dark matter-dark mediator coupling, α D , is fixed by the dark matter relic abundance, while the coupling of the dark sector to the SM is bounded from below by the thermalization floor. Previous estimates of these thermalization floors (e.g. [16,17]) have considered a subset of processes and/or studied equilibration at a fixed temperature.
As for standard WIMPs, the upper limit on the mass of DM is TeV-scale, arising when the interaction governing freezeout becomes non-perturbative. The precise value of this upper bound will depend in detail on the particle content of the dark sector. For instance, for DM freezing out via annihilations to massive dark photons, the upper bound depends on the structure of U (1) symmetrybreaking in the dark sector [18]. Perturbative unitarity constraints in specific models can further tighten the upper bounds on the DM mass (e.g., [19]). We will indicate in our parameter spaces where obtaining the correct relic abundance in our simple models requires the dark matter-dark mediator coupling to become non-perturbative, α D ≥ 1. This occurs for m χ ∼ 10−150 TeV in both simplified models, where the lower end of the mass range is for small DM-dark mediator mass splitting, and the upper end is for large splitting. The Sommerfeld enhancement (discussed in Appendix B) included in our freezeout calculation heavily sculpts this range. When the Sommerfeld effect becomes very large,  The orange region is where the dark sector is in thermal equilibrium with the SM at freezeout, while the blue region has the two sectors never in thermal equilibrium. In the green region, the SM and hidden sector were in equilibrium at some higher temperature (here near the QCD phase transition), but fell out of equilibrium by T f , so that the temperatures of the two sectors may drift apart (A.11). The hatched regions are near either the chiral or electroweak phase transitions, where our calculation is less reliable. Right: Thermal coupling regions for the scalar model as a function of the freezeout temperature, T f . Colors and hatching are as in the vector model. The kink in the blue region near 4 GeV is when the effectiveness of the equilibration from the ss → ff process near 30 GeV exceeds that from top processes near 200 GeV. For more details, see Appendix A.
our numerical freezeout calculation becomes less reliable, and we will further indicate these regions in presenting our parameter space. However, as the phenomenology does not undergo qualitative changes in this m DM TeV region of parameter space, we will not discuss it in detail. Meanwhile, the number of relativistic degrees of freedom that can be present at temperatures T 2m e ∼ MeV are restricted by BBN, which mandates that 2.3 < N eff < 3.4 [20]. When the dark sector is in thermal contact with the SM at temperatures T f MeV, we must then have both DM and the mediator be nonrelativistic by T ∼ MeV. We here impose the simple requirement m DM , m med > MeV. A more careful treatment of the regions shown in green in Fig. 2 where the dark sector has departed from equilibrium with the SM prior to DM freezeout would relax these bounds slightly. A detailed treatment of this region is interesting, but beyond the scope of this paper.
Direct detection. Both our vector portal and Higgs portal models have a leading spin-independent scattering cross-section with nuclei. Unlike for traditional WIMPs, the size of this cross-section is not directly related to the dark matter annihilation cross-section: it is proportional to the square of the portal coupling and can be parametrically small. We will demonstrate that both current and proposed direct detection experiments have the sensitivity to test cosmologically interesting values of the portal coupling. Currently, the best constraints on spin-independent DM-nucleus scattering come from XENON1T [21], LUX [22,23] and PandaX-II [24] at higher masses, while CDMSlite [25] and CRESST-II [26] set the strongest limits at lower masses. 1 We show the current limits, along with projections for several future experiments [29][30][31][32][33][34], in Figure 3. In the figure, we also present the neutrino floor for both xenon and calcium tungstate (CaWO 4 ) [35,36].
Indirect detection. In contrast to direct detection, results from indirect detection searches are insensitive to the (small) portal coupling, and test the dark matter annihilation cross-section directly. There are multiple sensitive probes of dark matter annihilation in the universe. The most important for our models are the Fermi-LAT limits on dark matter annihilation in dwarf galaxies [38,39] and Planck constraints on DM annihilations near recombination [1,40]. Charged cosmic rays are another important source of information about galactic DM annihilation, but are subject to much larger systematic uncertainties arising from their propagation within the galaxy. While AMS-02 measurements of the cosmic antiproton flux [41] can potentially give more powerful constraints on hadronic annihilation channels than searches with gamma rays [42], the difficulty in accurately determining propagation parameters remains a serious hurdle. We follow [43] in considering AMS-02 positron results [44], which can place bounds on leptonic channels where searches in photons have little reach, but neglecting antiproton searches, as they constrain channels for which the far less uncertain gamma-ray searches of [38,39] have good sensitivity. Meanwhile, CMB limits are mainly sensitive to the net energy deposited in the e − -γ plasma by DM annihilations near recombination [45], and are thus robust and nearly model-independent. The HAWC experiment can place constraints on very high dark matter masses [46] in the highly Sommerfeld-enhanced regime; these constraints are currently exceeded by the CMB constraints everywhere, but may become more important as HAWC collects more data, or our understanding of the Triangulum II dwarf galaxy, which dominates HAWC's sensitivity, improves [47,48]. In principle, H.E.S.S. should have sensitivity to our DM models when m χ ∼ TeV, but they do not provide enough information to allow their results to be reliably reinterpreted. 2 Accelerator. On the collider front, there are several potential discovery avenues for hidden sector dark matter. The direct production of DM (or of an invisible mediator) in events with large missing energy is no longer the leading signal, as we will demonstrate below. Rather, the leading accelerator signal is the direct production of the dark mediator, followed by its decay back to visible SM states. Mediators can be produced through rare Kaon and B-meson decays, directly through their interaction with electrons and quarks at LEP and LHC, at lower energy colliders such as Babar, and at beam dump  . The shaded gray region represents the current exclusion from the combined results of XENON1T [21], LUX [22,23], PandaX-II [24] CDMSlite [25], and CRESST-II [26]. Also shown in dashed lines are the projected limits from argon-based DEAP-3600 as well as its proposed 50 ton-year upgrade [29]; the xenonbased experiments XENON1T [30], LUX-ZEPLIN (LZ) [31], and DARWIN [32]; the CaWO 4 -based CRESST-III and its Phase 2 upgrades [33] and the projected limits from SuperCDMS [34] for both silicon and germanium for both the interleaved Z-sensitive Ionization and Phonon (iZIP) detectors (thin, dotted) and those run in high voltage (HV) mode (dashed). The coherent neutrino scattering floor is shown for both CaWO 4 and xenon [35,36]. The floor for silicon, germanium, and argon is very similar to xenon, while more complex materials, like the CaWO 4 in CRESST and the proposed EURECA [37] experiment, can have a substantially different ν floor. and other intensity frontier experiments such as NA62. They can also be produced in exotic Higgs decays [50,51]. Precision tests of Z and Higgs couplings can also constrain the mixing between dark and visible states.
Astrophysical and cosmological constraints on dark mediators. Beyond the standard suite of DM search strategies, models with long-lived dark mediator states face several additional constraints from astrophysical and cosmological observations. As the requirement that the dark sector be thermalized with the SM places lower bounds on the coupling of the dark mediator, these constraints will largely be important for the HSFO-HP model in the sub-GeV regime where small Yukawa couplings help increase the mediator lifetime. Most constraining here are cooling in Supernova 1987A [17,52], and early universe limits on the dark scalar lifetime coming from potential disruptions of isotope abundances produced during BBN or dilutions of neutrino and/or baryon abundances [53].

Vector portal
We first consider a simple vector portal model, containing a fermionic DM, χ, and a dark photon, Z D . This type of model has been studied extensively in the literature, especially to address cosmic ray anomalies (HEAT, PAMELA, and ATIC first [8,54], and more recently the Galactic Center excess [18,51,55,56]). In the following, we define our model and establish notation. We introduce a massive dark photon Z D , the gauge boson for a new dark U (1) Z D symmetry, that interacts with the SM through kinetic mixing with SM hypercharge [57,58]. The dark photon mass could arise from the Stueckelberg mechanism [59,60] or from a dark Higgs mechanism. For the sake of minimality, we will assume a Stueckelberg origin, so that the only dark sector particles in our model are the dark vector and the dark matter. Including a dark Higgs boson could open up additional annihilation channels, such as χχ → Z D h D , which could become the leading process in the regime m DM m Z D , m h D [18]; we discuss this possibility further in Sec. 4.4. The dark vector Lagrangian is thus given by where θ W is the Weinberg angle and is the dimensionless kinetic mixing parameter. Additionally, we introduce a Dirac fermion χ with unit charge under U (1) Z D and with mass m χ to serve as DM.
Making the standard field redefinition to diagonalize the hypercharge and Z D boson kinetic terms rescales the dark coupling g D =ĝ D / 1 − 2 / cos 2 θ W , and results in the following mass matrix for the neutral gauge bosons after electroweak symmetry breaking, in the basis (A, Z 0 , Z D,0 ). Here η ≡ cos θ W √ 1− 2 / cos 2 θ W and δ ≡ m Z D0 /m Z 0 , with m Z 0 the mass of the SM Z boson before mixing. The resulting massive eigenstates are with mixing angle The massive eigenvalues are We can now compute the couplings of the SM fermions and the DM with the Z D gauge boson: where Y , t 3 are the hypercharge and isospin of the (Weyl) fermion f . The physical Z boson acquires a (vector-like) coupling to χ: Note that this coupling is -suppressed, contrary to the corresponding coupling of the Z D . In fact, if we expand the couplings in (4.6) and (4.7) to leading order in we find (4.8) Our simple model can be described by four independent free parameters, which we take to be α D , , m χ and m Z D .

Thermal Freezeout and Indirect Detection
When DM is heavier than the dark photon, it can annihilate viaχχ → Z D Z D . This is the only annihilation channel in the small limit, and in this limit the thermally averaged cross-section for this annihilation, is independent of . For sufficiently heavy DM (m χ TeV), Sommerfeld enhancement can be important during freezeout, σv = σv 0 S 0 (v) , which we implement via a Hulthén potential as  described in Appendix B. Requiring that this reaction yields the observed relic abundance as measured by Planck, Ω DM h 2 = 0.1186 ± 0.0020 [1], fixes α D for each choice of m χ , m Z D . The same annihilation cross-section governs indirect detection signals. In order to assess the constraints from indirect detection, we utilize the measurement of dwarf galaxies from the Fermi-LAT and DES collaborations [38]. The energy flux of photons from DM annihilations in an astrophysical source can be expressed as where we have specialized to Dirac DM. Here dN/dE is the number distribution of photons from a single DM annihilation, and J is the astrophysical J-factor, describing the line-of-sight density of dark matter in the direction of the source [62]. We use the 41 dwarf galaxies within the nominal sample of [38] to obtain limits on the DM annihilation cross-section using the procedure outlined in Appendix C. The corresponding upper bounds on the cross section as a function of the dark matter mass m X and of the mass ratio m Z D /m X are shown in Fig. 4; the red regions show where the thermal relic abundance is excluded. These regions appear in two distinct places: the region in the lower right is where Sommerfeld enhancements are important, while in the upper left they are not (see Appendix B for details).
The flux of positrons observed in the AMS-02 experiment [43,44] can constrain photon-poor annihilation channels. In order to set constraints, we use the limit for one-step e + e − channels from [43] and compare this to σv × Br(Z D → e + e − ). In principle, considering all Z D decay modes would slightly improve this result, but achieving this mild improvement is beyond the scope of this work. The resulting exclusion is shown in orange in the right panel of Fig. 4.
Dark matter annihilation during the era of recombination can broaden the surface of last scattering and distort the CMB anisotropies through the injection of electrons and photons into the plasma. For WIMPs annihilating with a velocity-independent cross-section, the effect of this energy injection can be accurately encapsulated by a redshift-independent efficiency parameter f eff (m DM ), which depends on the DM mass and the species of particles produced by DM annihilations [63]. Planck results together with results from ACT, SPT, and WMAP limit f eff (m χ ) σv /m χ < 14 pb c / TeV [40], allowing for robust bounds to be placed on dark matter models. The f eff values for DM annihilation to pairs of SM particles have been computed in [40,64]. Due to the rather soft dependence on m χ for all branching ratios except for photons and leptons, we use the f eff values in [40] evaluated at m χ /2 for non-leptonic channels, together with f π ± eff from [64]. For leptonic channels we use where m Z D governs the branching ratios. The resulting CMB bound is shown in purple in the right panel of Fig. 4, as a function of the dark matter mass m X and of the mass ratio m Z D /m X .

Direct Detection
Direct detection experiments are an excellent test of this minimal model, and over much of the parameter space place the most stringent constraints on the portal coupling . The spin-averaged, nonrelativistic amplitude-squared for DM to scatter off of a nucleus is mediated by the exchange of both dark vector and Z bosons, and is given by where A is the mass number of the target nucleus, F 2 (E R ) is the Helm nuclear form factor [65,66] as a function of the recoil energy E R , s W ≡ sin θ, m N is the mass of the nucleus, and f are given by with Z the atomic number, and g X,u and g X,d the couplings of the boson X with up and down quarks in (4.8). Here we have retained the momentum dependence in the propagator from Z D exchange, as it is needed to accurately describe the scattering when m Z D µ χN v χ , where the DM-nucleus reduced mass is µ χN = m χ m N /(m χ + m N ). When the scattering amplitude does not depend on the DM velocity, then the event detection rate per unit detector mass in the experiment can be expressed as [67,68] where ρ χ is the local DM density, (E R ) is an experiment-specific selection efficiency, and η(E R ) is the mean inverse speed [68] defined by for which, following the experiments, we use the expression in Ref. [66]. 3 If the amplitude is independent of the recoil energy to leading order, it is reasonable to approximateM in (4.14), allowing for the particle physics contributions to the rate to be entirely factorized from the experimental and astrophysical inputs. Experimental results are typically expressed in terms of a cross-section that has been factorized in this manner and further simplified by defining an effective per-nucleon cross-section, facilitating comparison between different experiments. For the HSFO-VP model, this DM-nucleon cross-section is where we have defined the nucleon-DM reduced mass µ χn = m χ m n /(m χ + m n ). However, for the HSFO-VP model, In order to correctly account for this important effect, we will determine the excluded cross-section via where the function R is determined separately for each experiment. Given masses for the DM and dark vector, the relic density constraint fixes α D . The latest XENON1T [21], LUX [22,23], PandaX-II [24], Super-CDMS [69], CDMSlite [70] and CRESST-II [71] searches then determine the maximum allowed value of the portal coupling . We show these upper bounds in the left panel of Fig. 5. Sensitivity is greatest at small values of m Z D /m χ , thanks to the 1/m 4 Z D behavior of the nuclear matrix element. However, the sensitivity saturates when m 2 Z D 2m N E min (the threshold energy of the experiment), and the propagator in the matrix element is dominated by the momentum. Over a sizable region where m Z D /m χ 0.01 and m χ ∼ 10 GeV, current direct detection limits are sensitive enough to exclude values of the portal coupling at and below the thermalization floor. This region is shown in blue in the left panel of Fig. 5; see Appendix A.1 for details of its determination. Future direct detection experiments will be able to test this cosmological origin for DM over a broader range of DM and mediator masses. In the right panel of Fig. 5, we show the direct detection parameter space consistent with our HSFO-VP WIMP next door (in tan). As Fig. 5 shows, this cosmological history for DM can predict spin-independent cross-sections well below the neutrino floor.

Accelerator and other mediator constraints
Direct searches for the Z D are the leading terrestrial signal of this model. As summarized in [72][73][74][75] and in Fig. 6, there are many constraints on massive vector bosons kinetically mixed with SM hypercharge. Most of these results come from searches for rare meson decays, beam dump experiments, precision electroweak tests, direct production at BaBar or the LHC, and Supernova 1987A. Fig. 6 shows the parameter space for a vector portal WIMP next door as a function of m Z D and . As shown there, Supernova 1987A uniquely probes the thermalization floor in a limited range of dark photon masses at around few ×10 −2 GeV. Furthermore, especially at low masses, terrestrial searches for dark photons bound more tightly than the direct detection constraints do alone.   [32,33], will greatly cut into this range (green line), even excluding down to the thermalization floor near 500 MeV.
However, while the net impact of these direct mediator searches is generally subdominant to the combined constraints from direct and indirect detection for the minimal HSFO-VP model, it is important to emphasize that they provide complementary information. In particular, Fig. 6 shows that any dark photon discovered in meson decays or at high-energy colliders is sufficiently strongly coupled to the SM to populate a dark radiation bath in the early universe, regardless of the identity of dark matter. Further, as we will discuss below, simple extensions to the minimal model can suppress the direct detection cross-section, thus leaving dark photon searches as the leading terrestrial test of vector portal WIMPs next door.
A summary of all constraints is shown in Fig. 7 as a function of m χ and m Z D /m χ . As before, we show the union of -independent constraints from Fermi dwarfs, AMS-02 positrons, and the CMB with the black shaded region. The shaded green region denotes where the most important bound on comes from collider, beam dump, and supernova searches. Above the pink dashed line, the mediator was non-relativistic at dark matter freezeout, and thus the floor from thermalization can be much lower. As in Fig. 4, the solid green in the lower left corner has m Z D < 2m e , causing issues with BBN, and the brown region on the right of the plot illustrates where the freezeout coupling as determined with and without Sommerfeld enhancement deviate by a factor of 2.
We do not show constraints from mono-X searches at the LHC, since they are generically weaker than the constraints coming from direct detection experiments. In particular, the most stringent mono-X constraint arises from the ATLAS and CMS mono-jet searches performed with Run II data [76,77].
Values of as small as ∼ 0.1 are only probed in a small region of parameter space at around m χ ∼ 100 GeV and m χ ∼ (0.8 − 1)m Z D .

Beyond the minimal model of dark vector interactions
While we have worked with a minimal two-species model consisting only of fermionic DM and the vector mediator, the salient features of this model are representative of the behavior of a broad class of dark sectors with a vector mediator. In this section, we briefly discuss the modifications of the dark sector phenomenology obtained by introducing new dark degrees of freedom (or altering the assumed quantum numbers of DM), and argue that our minimal two-species HSFO-VP model provides a reasonable general guide to the characteristic sizes and locations of signals for vector portal WIMPs next door.
To begin with, any additional relativistic species in the thermal plasma at freezeout will contribute to the Hubble parameter, thereby requiring a mild increase of the value of α D needed to obtain the thermal relic abundance. The DM relic abundance is proportional to α 2 D × g * S / √ g * at freezeout. Neglecting the logarithmic sensitivity of the freezeout temperature to g * S , we can thus simply estimate the effect of adding additional equilibrated dark species by rescaling α D to absorb the shift in g * S / √ g * . This is a minor quantitative effect, particularly at relatively high DM masses where g * SM (T f ) 1. An excellent motivation for introducing additional dark species is to provide a dark Higgs mechanism to generate m Z D [8,18,78]. Our model assumed a Stueckelberg mechanism for simplicity. Using a dark Higgs to generate m Z D is a generic alternative scenario, but with a dark Higgs comes additional model dependence, especially through the choice of the DM's U(1) D quantum numbers. When the dark Higgs is light, it can furnish additional annihilation modes: depending on the spectrum, both the s-wave χχ → Z D h D process and the p-wave χχ → h D h D process can contribute to freezeout [78]. (We continue to assume that the vector portal coupling dominates the dark sector's interactions with the SM.) The additional annihilation modes change the specific value of α D needed to obtain the thermal relic abundance, generically by no more than an O(1) amount. These additional annihilation modes also alter the detailed cosmic ray spectrum for indirect detection. When p-wave contributions are important, the expected annihilation cross section for CMB and galactic signals can be decreased, generically by a factor of no more than O(1). Thus introducing these additional annihilation modes generally changes indirect detection signals quantitatively but not qualitatively. On the other hand, a dark Higgs can drastically impact direct detection. A Stueckelberg mass for the dark photon requires Dirac dark matter, and thus yields unsuppressed spin-independent direct detection cross-sections. However, a dark Higgs mechanism allows the Dirac spinor to split into two Majorana mass eigenstates, the lighter of which is dark matter [79]. If the mass splitting is small so that the Majorana states are nearly degenerate (pseudo-Dirac), then the leading spin-independent cross-section is now inelastic. Inelastic scattering is significantly more challenging to observe at direct detection experiments, but some signals are still possible [80].
However, as the mass splitting increases, the dominant direct detection signals come from elastic processes. These processes can arise at tree level, from the now axial-vector coupling of the DM to Z D . At relatively high dark vector masses, the axial components of the SM-Z D couplings are sizable, giving rise to spin-dependent cross-sections. The vector SM-Z D couplings yield spinindependent cross-sections suppressed by DM velocities or nuclear recoil momenta, giving small but still potentially interesting signal rates [81]. Elastic spin-independent cross-sections are also induced at one loop [82][83][84][85][86]. The size of this contribution is thus sensitive to the UV field content of the dark sector. Finally, while the coupling of the dark Higgs to the SM is sub-leading for thermalization, the exchange of the dark and SM Higgses gives a spin-independent cross-section, and, depending on the size of the dark Higgs-SM couplings, could provide the leading direct detection signal; for further discussion of Higgs-portal direct detection, see Sec. 5.2 below.
If there are additional light dark sector species, φ, then Z D can have open decay modes within the hidden sector, and when these are active one generally expects Br(Z D → φφ) 1. Letting Z D decay can eliminate (if φ is stable) or modify (if φ is unstable [43]) cosmic ray and CMB constraints on DM annihilation, and, often for either case, subject the mediator to the generally weaker terrestrial searches for Z D → invisible [87][88][89]. As the vector portal coupling is the leading interaction between sectors, these new light states will have longer lifetimes than the dark vector, which can potentially lead to stringent cosmological constraints. In particular, as BBN does not allow for an additional radiation species with the SM temperature T SM , a model that maintains equilibrium with the SM through the vector portal has limited prospects for including stable dark radiation. As illustrated by the small green region in the left panel of Fig. 2, a vector portal hidden sector equilibrated with the SM cannot decouple from the SM at sufficiently early times to permit significant departures of the HS temperature from T SM . Interestingly, a model with a portal coupling of ∼ 2 × 10 −9 and DM that freezes out after the chiral phase transition may provide a WIMP next door that permits additional dark radiation. As this model could eliminate indirect detection signals, it would open up interesting parameter space for GeV-scale DM, together with dark radiation signals that could be observable at CMB-S4 [90]. Verifying the existence of this dark radiation window, via a more detailed calculation of the thermal production rate of dark photons from the hadronic plasma near the chiral phase transition, is an interesting topic for future studies.
Our reference model assumes fermionic DM. If DM is instead a complex scalar, the story is qualitatively unchanged: the leading annihilation cross-section is s-wave, while the leading direct detection cross-section is spin-independent and unsuppressed. Once again, the introduction of a dark Higgs would make only minor changes to the indirect detection signals, while potentially introducing sizeable and model-dependent changes to the direct detection signals.
Finally, in our HSFO-VP model and the variants above, annihilations of only one representation of the dark U (1) symmetry are important during freezeout. Introducing more states in different representations and allowing coannihilation to be important in determining the relic abundance can significantly alter the phenomenology and open up different areas of parameter space, but represents a much greater departure from the minimal model discussed here.
To summarize, our minimal model provides a good guide to the essential physics of vector portal WIMPs next door. Many possible additions to the dark sector would change signals qualitatively, by ∼ O(1) amounts, e.g., through affecting Hubble. Indirect detection signals are especially robust, as adding additional annihilation channels, etc., generically changes cosmic ray signals quantitatively but does not suppress them significantly below expectations for an s-wave thermal relic. For vector portal models there is very little scope to eliminate indirect detection signals via Z D decays to dark radiation. On the other hand, direct detection signals are especially sensitive to the origin of dark symmetry breaking. The direct detection signals for the minimal model we present are maximally predictive; in a dark sector with a dark Higgs mechanism, direct detection signals can be suppressed by model-dependent amounts.

Higgs portal
Here we define a simple reference model for a Higgs portal WIMP next door, HSFO-HP. We consider a Majorana fermion dark matter, χ, with a scalar mediator, S, that interacts with SM states through a (small) Higgs portal coupling. A useful simple model is [91] (see also [5,92,93]) where we use the usual conventions for the Higgs potential, We should also add to this Lagrangian the interaction of the Higgs with the quarks, leptons and gauge bosons of the SM. These interactions will be inherited by the dark scalar through its mixing with the SM Higgs. In the Lagrangian in (5.1-5.2), we have imposed a discrete symmetry taking S → −S, χ → iχ, thus forbidding cubic and linear terms in S as well as a Majorana mass for χ. Imposing this discrete symmetry allows us to expose the essential physics of this theory with the minimum number of parameters.
In order for the fermions to be massive, both S and H must acquire nonzero vacuum expectation values (vevs), S = v s + s 0 and H = 1 Minimizing the potential gives analytic expressions for the vevs, The dark matter gets a mass of m χ = yv s , while the scalars have a simple mass matrix, yielding the mass eigenvalues and a mixing angle defined by where the latter equality holds only for 1 6 λ s v 2 s = λv 2 h and therefore the mass of the scalar S quite different from 125 GeV. In this regime, for small , cos θ ∼ 1 and tan θ ∼ sin θ ∝ , so either sin θ or can be viewed as a measure of the strength of coupling between the SM and dark sectors.
We can express the Lagrangian in terms of v h = 246 GeV, m h = 125 GeV and the four free parameters, y, m χ , m s , and sin θ. In terms of these parameters, the most important couplings for our discussion, to leading order in sin θ and for 1 6 λ s v 2 s = λv 2 h , are: It will sometimes be useful to refer to a fine-structure-like constant, α Y ≡ y 2 4π .

Thermal Freezeout and Indirect Detection
Assuming freezeout of the Majorana dark matter is governed by interactions entirely within the hidden sector (i.e., effects proportional to can be neglected), three diagrams dominate the process χχ → ss: tand u-channel exchange of χ, and s-channel annihilation through an off-shell s. The spin-averaged amplitude, integrated over final state phase space, can be written as where β i = 1 − 4m 2 i /s, and the factor of one-half for identical final state particles is explicit. This quantity, Ξ, is related to the cross-section by, σ cm = β s Ξ/(16πsv rel ), where the relative velocity is conventionally defined as Defining R ≡ ms mχ , the thermally averaged cross-section, to leading order in v 2 rel , can be written which is p-wave suppressed. This cross-section is then corrected by the Sommerfeld enhancement (B.4) to give σv = σv rel 1 S 1 (v rel ) . Thanks to the p-wave annihilation cross-section, this model does not yield observable signals from DM annihilations in halos, nor do late annihilations χχ → ss → SM deposit noticeable amounts of energy into the CMB. Thus standard indirect detection strategies do not constrain the minimal HSFO-HP model. Dark matter density spikes surrounding super-massive black holes could potentially boost p-wave DM annihilation to observable rates, yielding a point-like gamma-ray signal [91]. Additionally, there are regions of parameter space at high m χ where m s α 2 Y m χ /4 and DM annihilations can proceed through bound state formation and decay. These processes are s-wave, and thus although they are unimportant during thermal freezeout, they can leave an imprint on the CMB [94]. Assuming that the velocity of DM at recombination satisfies v R and bound states are accessible, an analytic solution for the bound state formation rate (via monopole transitions into S-wave bound states) can be written as [94] σv , and a factor of 1/4 appears due to only a single available bound state with spin-0 due to Majorana dark matter [94]. In practice, the Gamma functions in (5.10) yield roughly a cosecant function with maxima at α Y /R = m ∈ N regularized at the singular points by a tiny imaginary contribution. This bound state decay can be used to bound p-wave annihilating dark matter models with 1 2  where we use the same treatment as in Sec. 4.1 to derive f eff (m χ ). The excluded region is plotted in purple in Fig. 8. 4

Direct Detection
As in the HSFO-VP model, the HSFO-HP model has a leading spin-independent direct detection cross-section, and direct detection therefore provides a powerful test of this model. The direct detection cross-section is mediated by exchanges of both the dark scalar and the visible Higgs boson. In the regime where m s m χ , momentum exchange rather than the dark scalar mass dominates the dark scalar propagator. In order to account for this important recoil energy dependence, we follow the procedure discussed in Sec. 4.2. For the HSFO-HP model, the spin-averaged, non-relativistic amplitude-squared for DM-nucleon scattering is where m n is the mass of the nucleon and Here the nucleon matrix elements are where heavy quark flavors have been removed in the second line [95]. In principle f However, it is important to retain the recoil energy dependence in the dark boson propagator to accurately describe scattering rates when m s µ χN v χ ; to handle this, we use the procedure discussed previously in Sec. 4.2. In the left panel of Fig. 8, we show the maximum allowed value of the mixing angle sin θ consistent with the current direct detection bounds [21][22][23][24][25][26]. The blue region shows where direct detection experiments are probing values of sin θ below the thermalization floor. Additionally, the dashed blue-gray region corresponds to where the direct detection bound on sin θ implies the dark and SM sectors have decoupled and their formerly equilibrated temperature can drift (A.11). Once again, direct detection experiments can probe all the way down to the cosmological lower bound on the portal coupling in some portions of parameter space. The direct detection parameter space consistent with Higgs portal WIMPs next door is shown in the right panel of Fig. 8, where the tan region is defined by cutting out the region where m s is within 10% of m h .

Accelerator and other mediator constraints
The leading collider and accelerator searches for the HSFO-HP model are again direct searches for the mediator, s. Many different low-energy and collider observables are sensitive to Higgs-portal coupled scalars. Additionally, there are astrophysical and cosmological constraints on s when it becomes sufficiently light and long-lived. Our results for the current experimental reach for a light mediator, s, are presented in Fig. 10. In order to establish these results we need the SM branching ratios of Branching Ratio the scalar. These branching ratios are notoriously uncertain for scalar masses in the range between 2m π and ∼ 4 GeV (see [97] for more details). We adopt the hadronic branching fraction derived in Ref. [98], supplemented with a simple extrapolation in the very uncertain 1-4 GeV region, as shown in Fig. 9. Features near 1 GeV in Fig. 10 are sensitive to this choice, which is conservative for signals sensitive to muon decays; the overall lifetime is also important for determining projections for SHiP [99], MATHUSLA [100], and CODEX-b [101].
In the rest of this subsection, we will explain in the detail the constraints shown in Fig. 10. While other experiments have constrained this parameter space, such as KTeV [102] and NA48/2 [103], these results have been surpassed by the bounds from other experiments and will not be discussed here. For high scalar masses and large portal couplings, additional constraints from perturbativity and electroweak precision tests can be important [104,105]; however, these constraints are modeldependent and not in our main regime of interest, and we do not discuss them further here.

LHC (ATLAS & CMS)
: Heavy Higgs Searches -There have been many searches looking for Higgs-like bosons at the LHC. The strongest limits in the region 130 GeV < m s < 1000 GeV come from a search at ATLAS and two at CMS in the diboson decay channels [106][107][108].

LEP (OPAL, DELPHI, ALEPH & L3): Higgs Searches -
The LEP Working Group for Higgs boson searches combined the data from the four experiments at LEP to place very tight constraints on Higgs states that are produced in association with a Z from 15-115 GeV [111]. For scalar masses below 15 GeV, the tightest constraints come from L3 [110], but the ALEPH h → {invisible} search [109] sets slightly stronger limits below the muon threshold where the scalar is detector stable.
LHCb: B → K(s → µ + µ − ) -LHCb has made many detailed measurements of the important observables within B → K ( * ) µ + µ − . Two of these have been specifically interpreted to constrain   [113] (brown and purple, respectively); K → π + X, with X invisible at E949 [114] (teal); and the CHARM experiment at the CNGS beam-dump [115] (mustard). The blue, green, and red contours show the projected sensitivity for the SHiP, CODEX-b, and MATHUSLA experiments respectively [99][100][101]116]. In the lower left corner of the plot we indicate approximate constraints from SN1987A (dashed blue) [17] and late-time entropy injection (dashed green). Right: Bounds from exotic Higgs decays for fixed m χ = 2m s . In blue we show the region probed by searches for h → ss(aa) → 4µ [117], in red the region probed by searches for h → ss(aa) → 4τ [118] and finally in green the region probed by searches for h → ss (aa) → 2µ2b [119]. The region shaded in gray produces an exotic branching ratio larger than 34%. The corresponding dashed lines are the possible bounds obtained at the HL-LHC. While constraints are shown for m χ = 2m s , all of these bounds on sin θ scale roughly proportionally to √ m χ .
E949 & E787: K + → π + + invisible -The E949 collaboration at BNL has made the most accurate measurement of K + → π + νν [114] by measuring the decays of stopped K ± . In this study, they also reinterpret the results (including the data from E787) to place 90% CL upper limits on BR(K + → π + X, X → {invisible}) for m X < 125 MeV and 150 MeV< m X < 250 MeV. The outer radius of the barrel veto is 1.45 m [120], and all scalars that decay before this radius are assumed to be vetoed. The detector itself has 2/3 × 4π solid angle coverage. It is straightforward to reliably apply these limits on invisible decays to metastable particles. The NA62 experiment [121, 122], which will solidly establish BR(K + → π + νν) at the ∼ 10% level, should be able to significantly improve this limit (likely by an order of magnitude on sin 2 θ). However, the experimental design of NA62 is rather different from that of E949: notably, E949 measured stopped kaons while NA62 measures the decays of kaons in flight, so a detailed study is required in order to reliably assess the reach of this powerful new experiment [116].
CHARM: proton beam-dump -The CHARM experiment, with a 35 m decay region located 480 m downstream from the 400 GeV CNGS proton beam-dump, performed a search for axion-like particles [115]. They did not find any. However, as various hadrons are amply produced in the p-Cu interactions, this null result can be used to constrain light scalars with long lifetimes [123]. The approximate number of scalars that would be expected to decay within the detector can be estimated as The proposed SHiP experiment could also achieve excellent sensitivity to light scalars. At SHiP, unlike at CHARM, many kaons will be stopped. However, the overall luminosity is expected to be much higher, resulting in N B = 3.2 × 10 13 , N K ± = 2.9 × 10 16 , and N K L = 1.4 × 10 15 [124]. SHiP should be sensitive to any of the scalar final states, which greatly increases the total branching ratio into visible states above the pion threshold. Moreover, the SHiP experiment would have dec = 64 m and det = 50 m [124] resulting in a much improved sensitivity. We follow Ref. [125] in producing our estimated sensitivity from SHiP.
MATHUSLA and CODEX-b -The proposed MATHUSLA [100] and CODEX-b [101] experiments can be sensitive to long-lived scalars emitted in rare meson decays produced at the LHC. To estimate the reach of both experiments, we use the distributions of B-mesons produced in Pythia 8.223 [126], decaying to scalars following Ref. [116]. The event rate is normalized to σ bb = 0.5 mb. For MATHUSLA, the detector is treated as a 200 × 200 × 20 m box located 100 m above the interaction point and 100 m in the beam direction [100]. A 95% exclusion contour is shown assuming a flat 75% detection efficiency with no appreciable background for 3 ab −1 of data. For CODEX-b, we consider a 10 m cubic detector starting from 5 m in the beam direction (z) and offset by 15 m in x, and centered 2 m below the interaction point in y.
A third recently proposed experiment, FASER [127] could have sensitivity to extremely forward, boosted scalars produced in B decays [128]. While we do not reproduce the sensitivity here, the experiment would be expected to exclude additional territory below the τ + τ − threshold in between the sensitivity of SHiP and LHCb [128].
Supernova 1987A -The observed duration of the neutrino pulse from supernova 1987A places a restriction on how much energy it could have radiated into light scalars that escape the core of the supernova [17,129]. We follow the treatment in [17] to estimate this constraint. The total power radiated into scalars per unit volume is where the pion-nucleon coupling from low-energy scattering data is f πnn ≈ 1.0 [130], T SN is the core temperature, the nucleon Fermi momentum is given by p F = 3π 2 ρ SN 2mn 1 3 in terms of the core density ρ SN , G(x) is a function that contributes an O (1) factor (see [129] for the full expression), and to account for finite mass effects. A scalar only contributes to the energy loss if it escapes the core. The probability for a produced scalar to escape is where the mean free path of a scalar, λ s , can be approximated using the principle of detailed balance in terms of the equilibrium abundance of s at the core temperature, T SN , giving λ s ≈ ρ s,eq (m s , T SN )/P s . The total power radiated from the supernova in light scalars can be written as where V SN is the volume of the supernova core. The maximum allowed power radiated is P max ≈ 3 × 10 52 erg/s = 1.2 × 10 31 GeV 2 [131]. To produce the disfavored region shown in Fig. 10, we use the parameters of the fiducial model in [52]: V SN = 4/3 πR 3 SN with R SN = 10 km, T SN = 30 MeV, and ρ SN = 3.0 × 10 14 g/cm 3 .
However, several simplifying assumptions have been made here, such as the treatment of the stellar nuclei as degenerate, and a less approximate treatment (such as that done for dark photons in [52]) would produce refined results. As emphasized in [52], the robustly excluded region is smaller than the region excluded by the fiducial model, thanks to the uncertainty on the properties of the progenitor star, i.e., T SN , R SN , and ρ SN . Unlike in the case of the vector model, a naïve application of (5. 17-5.19) for the different progenitor star models in [52] leave no regions that are robustly excluded. For these reasons, we present the fiducial SN1987A bound with a dashed line in Fig. 10.
Cosmological constraints -Finally, for the lowest values of m s in the WIMP next door parameter space, the scalar can become cosmologically long-lived as all of its accessible decay modes are suppressed by tiny Yukawa couplings. Here, the deposition of macroscopic amounts of entropy into the SM during and after BBN can lead to unacceptably large decreases of the neutrino temperature relative to the photon temperature (as measured by N eff [53,132]), or unacceptably large disagreements between the BBN and CMB determinations of the baryon-to-photon ratio η B . A precise determination of these constraints requires a careful treatment of the temperature evolution of the scalars, which in this regime have thermally decoupled from the SM prior to their decay and thus undergo cannibal behaviour [133], and is beyond the scope of this paper. We have indicated the region where we estimate these cosmological constraints to be important with a dashed boundary in Fig. 10.
LHC searches for exotic Higgs decays -ATLAS and CMS have developed a program of searches for Higgs decays to two singlet scalars, s, or pseudoscalars, a. These decays give rise to the sig-natures h → ss (aa) → 4τ [118], h → ss (aa) → 2µ2τ [134,135], h → ss(aa) → 4µ [117], h → ss (aa) → 2µ2b [119], and h → ss (aa) → 4b [136]. Furthermore, searches for invisible Higgs decays constrain the invisible width of the Higgs boson to be smaller than 24% [137]. Finally, the ATLAS and CMS Run I combination of Higgs coupling measurements constrain the total exotic width of the Higgs boson to be smaller than 34% [138]. We show these constraints on the right panel of Fig. 10, having fixed m χ = 2m s . Presently, the constraints from global fits of Higgs properties (gray region) are generically stronger than the constraints from direct searches for exotic Higgs decays, with the exception of searches for h → ss(aa) → 4µ for masses below ∼ 5 GeV (see the blue region in Fig. 10 right). In the figure, we also show the possible prospects for probing the decays h → ss(aa) → 2µ2b (dashed green line), and h → ss(aa) → 4µ (dashed blue line) with 3000 fb −1 LHC data as studied in [139] and [73], respectively. 5 Also shown is the expected bound on the exotic Higgs width from Higgs fits (dashed gray line). At the HL-LHC, the 2µ2b signature is also expected to set more stringent bounds on light (pseudo-)scalars than the indirect bound on the exotic Higgs width. As the Higgs branching ratio into scalars is proportional to α Y sin 2 θ/m 2 χ (5.7), and we also have that the correct relic abundance gives roughly α Y ∝ m χ increasing the dark matter mass globally weakens the sin θ bounds in Fig. 10 right by approximately m χ /2m s .
A summary of all constraints is shown in Fig. 11 as a function of m χ and m s /m χ , where we present the contours for sin θ allowed by direct detection experiments. Also shown in the plot are the CMB constraints from bound state production (purple lines) and the regions where other constraints on the mediator, s, (collider, beam dump, cosmology, etc.) supersede the direct detection constraints (shaded green region). In blue, we indicate the region where the HS and the SM sectors were never in thermal equilibrium, and within the dashed blue-gray region the Higgs portal interaction equilibrated the two sectors but then fell out of equilibrium again prior to DM freezeout. As evident from the figure, over a sizable region of parameter space, current direct detection limits are sensitive enough to exclude values of the scalar portal coupling close to the thermalization floor.

Beyond the minimal model of dark scalar interactions
As with the vector portal model, the phenomenological features of our HSFO-HP model are representative of the behavior of more complicated Higgs portal WIMPs next door (see Sec. 4.4), under the assumption that dark matter is a fermionic state (and that CP is conserved).
One consequence of our decision to minimize the number of free parameters is that the mass of the dark matter is directly tied to the VEV of the scalar and therefore to the dark-visible Higgs mixing angle sin θ. This choice specifies a relationship between the portal coupling, , and the scalar-Higgs mixing angle, sin θ, for any values of of m χ , m s once the Yukawa coupling y = m χ /v s is determined through the relic abundance. Relaxing this assumption, e.g., by considering a model with a less stringent symmetry structure, allows the DM to have a bare mass m 0 independent of the scalar VEV. Unless there is some theoretical reason to expect m 0 , to be at the same scale as yv s , one would generally expect either one or the other to dominate. The case m χ ∼ yv s m 0 is well-described by the results shown above. However even when m χ ∼ m 0 yv s , the value of the Yukawa coupling y determined in our freezeout calculation will not change much: freezeout is dominated by the tand u-channel processes in (5.8), which are independent of v s . Thus the thermal relic result for y is largely determined by m χ , regardless of the origin of this mass. Therefore, at a fixed m χ , m s , and sin θ, introducing an independent bare DM mass results in a smaller v s , and thus a larger value for (5.6). As discussed in Appendix A, most of the scattering rates important for thermalization are dependent on sin θ, but the process ss → ff depends instead directly on . Decoupling the DM mass from v s thus makes the ss → ff process more important relative to the other processes. As can be seen from Fig. 1 right, the ss → ff processes peak near 30 GeV and die off sharply afterwards, so that adding a bare DM mass means that for T ∼ 30 GeV, thermalization becomes more efficient for the same value of sin θ.
Viewed in terms of , the case treated in detail here where m χ = yv s results in the largest parameter space above the thermalization floor, i.e., it allows thermalization for smaller . The case where m χ yv s has smaller sin θ for a fixed and thus typically requires larger to thermalize. Insofar as direct detection cross-sections depend on sin θ, and for most temperatures thermalization is controlled by sin θ, the ability of direct detection experiments to probe the thermalization floor is largely unaffected by the introduction of a bare Majorana DM mass. The exception is in regions T ∼ 30 GeV where thermalization is dominated by the ss → ff process so that the thermalization floor is located at smaller values of sin θ.
It is worth bearing in mind that our reference model, strictly speaking, predicts the scale of its own symmetry breaking phase transition. As discussed in Appendix A, this phase transition happens long before DM freezeout in our reference model; however, this is an important caveat to keep in mind when considering dark Higgs sectors with a different symmetry structure.
Unlike in the vector case, making the dark matter Dirac instead of Majorana only provides quantitative shifts and no qualitative changes, as this change does not alter the leading spin-independent matrix element for direct detection. A much more substantial change arises when one considers scalar dark matter: in this case the leading annihilation channel χχ → ss is now s-wave, and indirect detection signals become important and constraining [51,56]. (An unsuppressed s-wave annihilation cross-section could also arise for fermionic DM in the presence of CP-violation [140], or if the dark sector contains a light pseudo-scalar a in addition to the light scalar s, such that χχ → as can be an important annihilation channel [141].) Broadly speaking, the indirect detection signals and constraints in the presence of an s-wave annihilation cross-section are generally similar to the results found for the HSFO-VP model. In particular, constraints from the CMB are nearly identical, while limits from cosmic ray searches are qualitatively similar, achieving sensitivity to annihilation cross sections at the same order of magnitude.
Again, adding additional light states can allow for invisible mediator decays (and, in the presence of a leading s-wave annihilation cross section, mediator decays into stable dark states would allow indirect detection signals to be re-suppressed). Also as before, these additional dark sector states can face strong cosmological constraints; these can become especially acute for low-mass unstable dark sector states, which can become very long-lived when the Higgs portal is the leading interaction between sectors. Importantly, Higgs portal WIMPs next door are amenable to larger temperature drifts between the SM and hidden sector temperatures (see right panel of Fig. 2), which allows more scope for relativistic species at BBN. Stable dark radiation species would give visible signals in CMB-S4.
To summarize, one of the major differences between vector portal and Higgs portal WIMPs next door is that the Higgs portal models offer more opportunities to include dark radiation. There is slightly more model-dependence in the detailed location of the thermalization floor, as different choices in constructing the dark Higgs sector can alter the relationship between sin θ and . Straightforward extensions and variations of our minimal HSFO-HP model allow for the introduction of an s-wave annihilation cross-section and thus reintroduce indirect detection signals and constraints, while leaving direct detection signals qualitatively undisturbed.

Summary and conclusions
In this paper, we have comprehensively assessed the current constraints on and discovery prospects for a class of minimal hidden sector freezeout models, where the DM relic abundance is set by the freezeout of a single DM species, χ, into a dark mediator, φ. We consider the well-motivated scenario where the leading interaction between the dark sector and the SM is renormalizable, for two simple reference models, HSFO-VP, where the mediator is a dark photon kinetically mixed with SM hypercharge, and HSFO-HP, where the mediator is a dark scalar that mixes with the SM Higgs boson. In both cases, the interaction with the SM renders the mediator cosmologically unstable. Experiments across the cosmic, intensity, and energy frontiers provide complementary information about the nature and cosmic history of hidden sector DM in these reference models, and leading signals can can differ substantially from models of more traditional WIMP-like dark matter. We have carefully considered the cosmology of these HSFO models. When the interaction between the SM and dark sectors is sufficiently strong to ensure that the two sectors achieve thermal equilibrium prior to dark matter freezeout, the cosmology is highly predictive. When the leading interaction between the dark sector and the SM is renormalizable, this minimal cosmology is additionally UV-insensitive: the scattering that works to equilibrate the two sectors becomes increasingly important, relative to the Hubble rate, as the universe expands. Thus, as we emphasize, our simple models of HSFO define a minimal and robust dark cosmology. We define the WIMP next door as dark matter that freezes out from a dark radiation bath in thermal equilibrium with the SM. One major consequence for WIMPs next door is the existence of a new cosmological lower bound on the portal coupling, the thermalization floor min (T f ): the minimum value of that allows the dark sector to reach thermal equilibrium with the SM before DM freezes out at the temperature T f . We provide an initial computation of this thermalization floor for both Higgs and vector portal couplings. This bound is significantly more stringent than the bound from BBN over the vast majority of parameter space, and can be terrestrially interesting. While obtained in the context of our minimal models, these results should generally serve as a good guide for thermally populated DM in more general 'next door' hidden sectors.
WIMPs next door provide a sharply predictive scenario for hidden sector DM. Requiring that the dark radiation bath attains thermal equilibrium with the SM prior to DM freezeout enforces a relationship between the temperatures of the two sectors, so that the parameter space is bounded and clearly defined. Both the DM and mediator masses are bounded from below by BBN and from above by perturbativity, while the coupling of the mediator to the SM is bounded from below by min (T f ).
Dark matter direct detection experiments can access a significant portion of the parameter space for WIMPs next door, and in some regions can probe all the way down to the thermalization floor min (T f ). Low-mass (m 10 GeV) WIMPs next door with unsuppressed s-wave annihilation crosssections (like in the HSFO-VP model) can be robustly excluded from their impact on the CMB. On the other hand, WIMPs next door with a velocity-suppressed p-wave annihilation cross-section (like in the HSFO-HP model) predict interesting direct detection signals for DM candidates in the lowmass range that offer an attractive target for low-threshold direct detection experiments. Additionally, unlike in the case of standard WIMPs, a notable fraction of the viable WIMP next door parameter space dwells underneath the coherent neutrino scattering floor, providing a target for future direct detection experiments that would need some ability to distinguish these signals from the neutrino background.
The leading accelerator signature of WIMPs next door is the production of dark mediators. A variety of experiments currently constrain both kinetically-mixed vectors and Higgs-mixed scalars.
Both existing experiments, such as LHCb and NA62, and proposed experiments, like SHiP, MATH-USLA, CODEX-b, and FASER, project sensitivity to significant regions of unexplored territory, and will either lead to a revolutionary discovery or greatly improve the constraints on this parameter space. At higher masses, the multi-purpose LHC experiments have the best opportunities to discover the mediators in exotic Higgs decays (Higgs portal) or through direct production (vector portal). By contrast the traditional LHC mono-X searches give little hope of finding DM.
One avenue for future work is improving on our estimates of the thermalization floor, most critically through the SM's two phase transitions. In particular, a careful treatment of dark photon production through the chiral phase transition could be important for understanding windows for dark radiation and thus low mass (m 10 GeV) vector portal WIMPs next door.

A Thermal (De)coupling
In this appendix, we describe in detail our estimate of the minimum portal coupling necessary to thermalize the hidden sector with the Standard Model in the early universe, and point out some ways to improve on our treatment. We always assume that the particles of the hidden sector (i.e., the DM, χ, and the dark mediator, φ) rapidly thermalize among themselves and can be characterized by a single temperature. As we focus on the regime where the mediator forms a radiation bath at the time of equilibration, this assumption is well justified.
We begin with some general comments. First, we are mainly interested in the thermal interaction rates between particles at temperatures T m, where classical statistics do not apply. Once final state blocking and/or enhancement factors can no longer be neglected, the evaluation of collision terms becomes significantly more technically involved. Fortunately, the SM thermal bath is dominated by fermions: empirically, classical statistical ("Maxwell-Boltzmann") treatments of relativistic scattering processes involving fermions provide a reasonable approximation to the full quantum statistical expressions, agreeing within a factor of 2 (see, e.g., [11]). Thus we employ classical statistics to evaluate the rates for 2 ↔ 2 scattering processes like φf → (g/γ)f , φf → hf , and the crossed processes φ(g/γ) → ff , φh → ff , etc.
The other major simplifying approximation we make is to neglect 2 ↔ 2 scatterings with EW gauge bosons. This is a good approximation thanks in large part to the sheer numerical dominance of quarks in the SM plasma, combined with α s > α 2 and the larger color factors present in QCD scattering amplitudes. The processes φf → (W, Z)f and their crosses are thus numerically unimportant compared to φf → gf at high temperatures at our level of precision. At temperatures T m W , the W , Z masses render these scatterings irrelevant. Meanwhile all-bose processes such as φV → V V are only important for a small range of temperatures at and below the electroweak crossover T c ≈ 160 GeV [142] before Boltzmann suppression kicks in. A study of dark mediator production from electroweak boson scattering in this regime is interesting, but involves a careful treatment of (evolving, nonperturbative) thermal masses, and is beyond the scope of this paper. Above the electroweak crossover, the leading scattering processes that mediate thermalization have a different structure, as discussed further for each model below.
We incorporate three-loop running of α s above the chiral phase transition (everywhere in this work, we use Λ QCD ≡ 300 MeV). However, below the QCD phase transition, 2 ↔ 3 pion processes, e.g., π + π 0 → π + π 0 Z D , dominate. Due to the qualitative similarities, these are lumped into our 2 ↔ 2 processes in the discussion below. In order to estimate these processes, we expand the chiral Lagrangian to leading order in {p, m π }/4πf π to compute the relevant cross-sections. As the thermally averaged cross-section receives important contributions from values of s where this expansion is no longer reliable, we introduce a simple regulator that ensures σ(s) has physical high s behavior (σ(s m π ) ∝ s −1 ). As the extremely broad QCD σ (f 0 ) resonance would be expected to perform the bulk of the unitarization of pion scattering, we define our multiplicative regulator to be where the exponent ξ is chosen to enforce the desired UV behavior and m σ (Γ σ ) are set to 500 (600) MeV. While this assumption is grounded in physical expectations, it is a somewhat arbitrary choice and a different regulator could modify the results significantly. Further, while pion processes are included, kaons and other light QCD resonances, have been neglected. These states have a higher Boltzmann suppression, which should generally make their contribution subleading, but their contribution may still be considerable for scalar production thanks to the larger strange Yukawa coupling. Lastly, near the phase transition the strong self-interactions of the hadronic plasma likely make significant thermal corrections to the mediator production rates. For these reasons, our min (T eq ) results in this region are much more uncertain than the roughly factor of two precision that we have elsewhere. In Figs. 1 and 2, we shade temperatures near both the electroweak (T ∼ 160 GeV) and QCD (T ∼ 300 MeV) phase transitions to highlight the large uncertainties in these areas. At the level of precision we are using, we can check whether a process is in equilibrium simply by comparing it to the Hubble rate, requiring Here the effective number of relativistic degrees of freedom g * (T ) includes degrees of freedom from the dark mediator as well as those of the SM. We have checked that this simple equilibration criterion reproduces the results of a full numerical treatment of the energy transfer rate between sectors to within an O(1) factor; see also [143,144].
The ratio of scattering rates Γ 1↔2 /Γ 2↔2 in both HSFO-VP and HSFO-HP models is shown in Fig. 12. Regions below the dashed pink line in Fig. 12 have a radiation bath at freezeout. The figure demonstrates that 2 ↔ 2 scattering processes dominate thermalization in almost all of this region.

A.1 Vector Decoupling
In our vector model, the 2 ↔ 2 rate is dominated by the scattering processes Z D f → g/γf , Z D g/γ → ff , Z D f → hf , and Z D h → ff . For each of these cross-sections, σ 12→34 , the thermal average given by Maxwell-Boltzmann statistics is [145] σ 12→34 v T n eq 2 = ∞ √ s min and K n is the modified Bessel function of the second kind.
This expression for the cross section is related to the total reaction rate by , with the larger splitting for m χ ∼ 10 TeV and the smaller splitting for m χ ∼ 1 GeV. 6 The rate of inverse decay processes can potentially be larger than the 2 ↔ 2 rate. We continue to neglect final state blocking and stimulated emission factors in evaluating this rate. Here this approximation is reasonable as the 1 ↔ 2 rate becomes most important in comparison with the 2 ↔ 2 rate as T ∼ m Z D . We include all SM species in evaluating this rate, including EW gauge bosons. We compute the 1 ↔ 2 rate as using the Bose-Einstein distribution f eq (E) for the Z D . This yields where Γ Z D ∝ 2 is the zero-temperature width of the dark vector. Here in the top line we have used the Bose-Einstein result for the relativistic n Z D ,eq , while the bottom line gives the Maxwell-Boltzmann result. At low temperatures, the Bessel function ratio asymptotes to unity. As the interaction rate for the 2 ↔ 2 processes scales as Γ Z D int,2↔2 ∝ T in the UV and Γ Z D ∝ m Z D , this results in a rough parametric scaling of so decays and inverse decays are typically unimportant for thermalization unless T f m Z D . A comparison of the minimum allowed value for only the 2 ↔ 2 processes (which do not depend on the vector mass) and only the 1 ↔ 2 process for fixed values of m Z D /T f is shown in the left panel of Fig. 13. As the width of the dark photon rapidly increases for fixed at m Z D ∼ m Z , much smaller values of can thermalize the two sectors in this region. Below the QCD confinement scale, 2 ↔ 3 pion processes briefly dominate, but after they become subdominant near temperature of ∼ 50 MeV, 2 ↔ 2 process becomes proportional to α EM , and decays and inverse decays become relatively more important for thermalization. For m Z D ∼ m χ T f , the large width of Z D allows for this to dominate the thermalization.
The narrow green region in Fig. 2 corresponds to a scenario where at higher temperatures near 200 MeV, the hidden sector and standard model were in thermal equilibrium, but have since decoupled. The temperatures of the two sectors are then allowed to drift apart. For the vector model, this is a very narrow region of parameter space, so we defer our detailed discussion of this interesting region to the scalar decoupling section near (A.11).
Finally, depending on the origin of the dark vector mass, the dark vector model may implicitly contain a symmetry-breaking phase transition where the mass of the dark vector is generated. In this case, when both SM and dark sectors are in the unbroken phase, the leading scattering process responsible for bringing the two sectors into thermal equilibrium is ff → χχ [146]. As this process depends on the dark Yukawa coupling, and α D α S , for the purposes of determining the minimal portal coupling that can yield thermalization, the unbroken UV interaction rate is unimportant compared to the interaction rate after both electroweak symmetry breaking and dark symmetry breaking.

A.2 Scalar Decoupling
The dark scalar model predicts the critical temperature, T c , of its phase transition from the symmetric vacuum ( S = 0) to the broken vacuum where S develops a VEV. This phase transition occurs comfortably prior to DM freezeout, T c T f , as we now demonstrate. This model exhibits a secondorder phase transition, so T c occurs when the second derivative of the thermal potential at the origin changes sign. To estimate the critical temperature and understand its relation to other mass scales in the dark sector, it suffices to consider the one-loop approximation to the thermal effective potential for S, yielding Until the electroweak phase transition, the scalar's only tree-level interactions with the SM are with the Higgs multiplet. In the unbroken phase, above both the electroweak and dark sector phase transitions, there is thus a single process controlling the equilibration of the two sectors, ss → hh 7 . As this process involves only bosons with masses m i T , to accurately determine the interaction rate it is necessary to use Bose-Einstein statistics.
Using the techniques of [147], the thermally averaged scattering rate can be expressed as an integral over the total CM energy-squared, s, and p, the magnitude of the three-momentum of the CM frame in the rest frame of the plasma. Defining, as usual, the scattering rate Γ U V as the collision term 7 The crossed process sh → sh also assists thermalization by contributing to the energy transfer rate between sectors.
While this process does not change the number of dark particles, there exist rapid number-changing s self-interactions that serve this function. The sh → sh scattering rate is ∼ 10% − 20% larger than the ss → hh scattering rate, but provides a highly subleading contribution to the energy transfer rate once the temperatures of the two sectors are similar. Thus it is sufficient to estimate thermalization based on the scattering ss → hh. divided by the (equilibrium) number density of one of the initial state particles, we have Here p 0 ≡ s + p 2 , β = 1/T , and β i ≡ 1 − 4m 2 i /s. We evaluate this integral numerically. It diverges as the lower limit on s is taken to zero, s 0 → 0, reflecting the divergence in the Bose-Einstein distribution f (E) as E → 0. This divergence is regulated by the thermal masses of the scattering particles, s 0 = 4 max(m H (T ) 2 , m s (T ) 2 ). With m i (T ) ∝ T , we find Γ U V ∝ T , as we must. The size of the UV scattering rate thus depends indirectly on couplings internal to the two sectors through their role in determining the thermal masses of S and H. Smaller thermal masses cut off the divergence at a lower value of s 0 , and hence increase the rate. The relatively large couplings of the SM Higgs to the top quark, electroweak gauge bosons, and (to a lesser extent) itself ensure that m H (T ) determines s 0 , making Γ U V relatively insensitive to the detailed couplings of s within the hidden sector. We find that, for a fixed value of , the two sectors will thermalize in the unbroken phase only if they could also thermalize within the broken phase as well. In other words, the lower bound on that we find from requiring Γ U V (T ) = H(T ) at high temperatures is subdominant to the lower bound found by requiring Γ IR (T ) = H(T ) at temperatures below the phase transitions. Thus to understand the process of thermalization through the s-h interaction, it suffices to study the scattering rates in the broken phase.
At temperatures below the electroweak phase transition, but before the dark phase transition, the dominant processes are ss → ff processes via a SM-like Higgs mediator. By far the most important contribution here is ss → bb, which for temperatures near T ∼ 30 GeV benefit from the s-channel enhancement for transit through a nearly on-shell Higgs (regulated with a width set to the SM value Γ h = 4.15 MeV, with thermal effects neglected). 8 As this process depends on rather than sin θ, decreasing the dark matter mass enhances this rate, and for temperatures T = O(10 GeV) it becomes the dominant factor in determining whether the two sectors have ever been in equilibrium; see Fig. 1. After both the SM Higgs and the hidden sector scalar have VEVs, many processes can contribute to thermalization. Here the dominant ones are sh → ff , sf → hf , sg/γ → ff , sf → g/γf , as well as the ss → ff processes discussed above, which are unaffected at O (sin θ) by the dark phase transition and thus remain important in the broken phase. For temperatures T m t , mediators may also be produced through sg ↔ gg. This rate is logarithmically divergent, and we estimate it by cutting off the log divergence with a finite thermal mass for the gluon. Our estimate indicates sg ↔ gg is a subdominant contribution. Below the QCD phase transition, ππ ↔ ππs processes 8 The treatment of ss → ff thus also includes the important process where ss → h; as this process is not considered separately, there is no problem with double counting. dominate briefly. Our treatment of these thermal scattering rates follows that outlined in the case of the vector model above.
Once again, decays and inverse decays become important when m s T f , as can be seen in the right panel of Fig. 13. Our treatment of the 1 ↔ 2 scattering rates in the Higgs portal model follows the treatment described for the vector portal above. Again, we include all SM contributions to the scalar width.
As the hidden sector scalar couples to SM fermions with strength proportional to the fermion masses, the 2 ↔ 2 interaction rate Γ drops rapidly after crossing a fermion mass threshold as these massive particles drop out the thermal bath. Thus, it is possible that for some other T > T f , (A.2) is satisfied, but not at T f , so that the SM and the hidden sector were at one point in thermal equilibrium, but have since decoupled. When this happens, their temperatures drift apart as This region is shown in green in the right panel of Fig. 2. Tracking the detailed temperature evolution of the hidden sector in this region, which can involve cannibal behavior when the scalar is sufficiently massive and long-lived [148], is interesting, but beyond the scope of this paper.

B Sommerfeld Enhancement
When DM can interact via the long-range exchange of light mediators, the annihilation rate can exhibit a large enhancement over the tree-level rate, especially at low DM velocities [149][150][151][152]. This Sommerfeld enhancement is most pronounced when three basic scenarios are satisfied: the dark finestructure constant, α D , is large; the DM velocity, v, is small; and the mediator is much lighter than the dark matter, R = m φ m DM 1, all of which can be realized by heavy, thermal dark matter with a light mediator. It is common to define the Sommerfeld enhancement through the factorized formula [8,153,154], σv = S(v) (σv) tree , (B.1) where σv is the full cross-section, (σv) tree is the tree-level cross-section, and S(v) is the velocitydependent Sommerfeld enhancement.
To evaluate the Sommerfeld enhancement to DM annihilations, we make use of the analytic approximation obtained by replacing the Yukawa potential with the Hulthén potential [155], where δ = π 2 m φ 6 .

(B.4)
These analytic results from the Hulthén potential provide a good approximation to scattering from the true Yukawa potential [155] except in the resonant regime where disagreements can become numerically larger. We incorporate the Sommerfeld effect in two different ways. First, the Sommerfeld enhancement can become important during freezeout, especially at large DM masses. In this case, the increased annihilation from the Sommerfeld effect will reduce the size of the coupling constant necessary to achieve the correct relic abundance. We include Sommerfeld enhancement during freezeout by numerically solving the equation using an adaptive fifth-order Cash-Karp Runge-Kutta technique. Here σv s = σ 0 S 0 (v) and x . We further use the approximations S 0 (v) ≈ S 0 (v c (x)) and S 1 (v) ≈ S 1 (v c (x)).

(B.6)
In Figs. 4, 5, 7, 8, and 11 we will denote (in brown) the region where thermal freezeout with and without the inclusion of the Sommerfeld enhancement gives results for the dark fine structure constant that disagree by more than a factor of 2, Outside of this region, the approximation used in (B.6) proves very accurate. Deep within this region, the true coupling is typically smaller than that predicted by the freezeout calculation of (B.5), and in our approximation Sommerfeld resonances will be improperly positioned, by an even larger amount than from the use of the Hulthén potential. We further note that for very large couplings and/or very near resonances, the Sommerfeld enhancement as estimated in (B.3) can violate partial wave unitarity [156], and, again, the condition (B.7) reliably insulates us from this region. The Sommerfeld enhancement also may greatly affect the indirect detection of dark matter. Both today and at the era of recombination, dark matter moves very slowly and the Sommerfeld enhancement can substantially increase the annihilation rate. As the Sommerfeld enhancement decreases monotonically with increasing velocities, we will conservatively err on the side of assuming larger velocities. Dark matter in the Milky Way have relative velocities on the order of 10 −3 , and we use the conservative value of v GC = 1.7 × 10 −3 , (B.8) which corresponds to relative velocities of 500 km/s for determining the Sommerfeld enhancement for AMS-02. Dark matter in the smaller dwarf galaxies have characteristic velocities on the order of 10 −5 -10 −4 [157]. In our treatment of dwarf galaxies, we will conservatively use a uniform relative velocity across dwarfs of v dwarf = 10 −4 , (B.9) in computing the Sommerfeld enhancement. Especially in the case of faint dwarf spheroidal galaxies with a small half-light radius, such as Draco II and Segue I recently discovered by the Dark Energy Survey [158,159], this choice of characteristic velocity could significantly underestimate the Sommerfeld enhancement. Lastly, at the time of recombination, the characteristic velocity of dark matter was still dictated by its red-shifted temperature rather than by virialization within a structure. In particular, after the dark matter decouples from the thermal bath in the early universe at T KD , its velocity can be expressed as where T CM B ∼ 0.27 eV, and we have imposed the bound from Lyman-alpha forest data requiring T KD 100 eV [160]; the parameter space of interest in this work yields values for T KD well above this lower limit. 9 For simplicity, in this work we fix v CM B = 10 −7 , (B.11) which for mass ratios of interest falls well into the regime where the Sommerfeld enhancement does not grow any further with decreasing velocity.

C Bounds from dwarf galaxies
In this Appendix we discuss the procedure we use to set limits on our DM models from Fermi's search for DM annihilations in dwarf galaxies. We consider the 41 dwarf galaxies within the nominal sample of [38]. The Fermi collaboration provides a log-likeihood ratio (LLR) for a signal + background assumption to background only as a function of the injected signal for each of the 41 dwarfs in each of the 24 common energy bins. To use these LLRs in combination to constrain a different signal model, it is necessary to account for the correlation of the systematic uncertainty on the J-factors between dwarfs. As this information is not provided, we model it by considering a 0.5σ downward shift from the naïve central value (using σ unmeasured = 0.6), as this was determined to replicate constraints on the bb and τ + τ − annihilation models fairly reliably at lower masses, while slightly underestimating constraints at higher masses (see Fig. 14). In particular, the region of HSFO-VP parameter space currently excluded by Fermi dwarfs is reliably determined by this choice. After this shift, for each energy bin the 41 dwarf measurements are combined to form a net LLR as a function of the signal injected into that bin. For each point in the m χ vs m Z D parameter space, a dark matter annihilation process, χχ → Z D Z D → {all}, is generated within Pythia 8 (v8223) 9 See Ref. [161] for discussion of related models in the low TKD regime. [126] 10 and the resulting gamma ray spectrum is tabulated into the 24 energy bins. The population of these energy bins scales with σv . Combining these bins into a χ 2 with one d.o.f. (the overall DM annihilation rate) places limits on σv which are shown in Fig. 4.