Phenomenology of ELDER Dark Matter

We explore the phenomenology of Elastically Decoupling Relic (ELDER) dark matter. ELDER is a thermal relic whose present density is determined primarily by the cross-section of its elastic scattering off Standard Model (SM) particles. Assuming that this scattering is mediated by a kinetically mixed dark photon, we argue that the ELDER scenario makes robust predictions for electron-recoil direct-detection experiments, as well as for dark photon searches. These predictions are independent of the details of interactions within the dark sector. Together with the closely related Strongly-Interacting Massive Particle (SIMP) scenario, the ELDER predictions provide a physically motivated, well-defined target region, which will be almost entirely accessible to the next generation of searches for sub-GeV dark matter and dark photons. We provide useful analytic approximations for various quantities of interest in the ELDER scenario, and discuss two simple renormalizable toy models which incorporate the required strong number-changing interactions among the ELDERs, as well as explicitly implement the coupling to electrons via the dark photon portal.


Introduction
Cosmological observations at a variety of length scales, from individual galaxies to the Hubble scale, indicate that most of the matter in the universe is in the form of dark matter (DM). DM cannot consist of any of the known elementary particles, and its existence provides solid experimental evidence for physics beyond the Standard Model (SM). The microscopic nature of dark matter is one of the major mysteries in fundamental physics. For many years, both theoretical work and experimental searches for dark matter focused on a short list of possible candidates independently motivated by particle physics-primarily QCD axions and weaklyinteracting massive particles (WIMPs) realized within supersymmetry or other extensions of the SM at the weak scale. Despite decades of experimental effort, no evidence for these candidates has been found. While neither WIMP nor axion dark matter is ruled out and the experimental searches are ongoing, there has been renewed interest in exploring alternative particle dark matter candidates.
A promising new direction is to consider models in which dark matter particles have strong number-changing self-interactions [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. If the DM is a thermal relic, its current density in such models can be determined either by the cross section of the number-changing self-interaction processes ("Strongly-Interacting Massive Particle", or SIMP, scenario [18]) or by the cross section of elastic scattering between the DM and SM ("Elastically Decoupling Relic", or ELDER, scenario [7]). In both cases, the observed DM density is naturally obtained if the mass of the DM particles is parametrically close to the QCD confinement scale, m DM ∼ 10−100 MeV. This leads to an attractive particle physics framework: a "dark sector" of fields not charged under the SM gauge groups, containing a non-Abelian "dark QCD" gauge group that confines at a scale similar to Λ QCD . The proximity of the SM and "dark" confinement scales may be due to a discrete symmetry relating the dark QCD gauge coupling to the SM g 3 at a high energy scale [19][20][21][22][23]. The dark matter may then consist of mesons that emerge from dark QCD upon confinement [2]. If the dark sector also contains an Abelian gauge field, kinetic mixing between this field and the SM electromagnetic field naturally provides the requisite interaction between the dark matter particle and the SM, via the dark photon portal [4,8].
The goal of this paper is to study the above possibilities in more detail, in particular, the ELDER scenario proposed in Ref. [7]. In Ref. [7], we demonstrated the viability of this scenario in a general framework, without reference to a specific model of either the dark sector or the portal connecting it to the SM. Instead, we used a simple parametrization of the DM number-changing self-scattering and DM-SM elastic scattering cross sections. Moreover, the analysis of Ref. [7] was primarily based on numerical solution of Boltzmann equations. Here, we expand that analysis in several directions: • We provide an approximate analytic solution to the Boltzmann equations that describe the evolution of the ELDER dark matter density during the epoch of its kinetic decoupling from the SM. This in turn leads to precise analytic estimates of ELDER relic density, and hence the model parameters required to obtain the observed dark matter abundance. We also combine these estimates with unitarity considerations to obtain a model-independent upper bound on the ELDER dark matter mass. This is the subject of Section 2.
• We consider the phenomenology of the ELDER scenario with a dark photon portal mediating the interactions between the ELDERs and the SM. We find that the model makes a remarkably robust prediction for rates expected in direct-detection experiments. This prediction has no free parameters beyond the ELDER particle mass, and is completely independent of the details of the dark-sector self-interactions. The reason is that the ELDER relic density with this portal is determined by the cross section of elastic scattering of dark matter particles on electrons, which is precisely the same process used for direct detection in the MeV-GeV DM mass range. This feature is unique to ELDERs. Likewise, this scenario provides firm predictions for dark photon searches. Together with the well-known "thermal targets", the ELDER and SIMP predictions define a well-defined target region in the parameter space for direct-detection and dark photon experiments, bounded from all sides. These findings are reported in Section 3.
In particular, Figs. 3 and 4 encapsulate the main results of this paper.
• We discuss two simple perturbative models for the dark sector, which realize the ELDER scenario with the dark photon portal; see Section 4. These can be thought of as toy models that describe interactions among low-lying mesons created by confining gauge dynamics in the gauge sector.
Details of the Boltzmann equations, an approximate analytic solution for kinetic decoupling, and some useful formulas for thermally-averaged rates, are collected in the Appendices.

ELDER Dark Matter
Consider a particle χ with mass m χ 1 . The χ particles can undergo the following processes: 1. Elastic scattering: χ + SM ↔ χ + SM, where "SM" stands for any of the Standard Model particles. (In practice, the important SM states are those with mass below m χ ; for ELDERs, this will typically include electrons, photons, and neutrinos.) 2. Annihilations to SM: χ + χ ↔ SM + SM.
We assume that in the early universe at temperatures above m χ , all four reactions are "active", i.e. occur in the plasma at rates Γ > H. This means that the ELDERs have a thermal energy distribution (thanks to reaction 4), zero chemical potential (reaction 3), and temperature equal to that of the SM plasma (reactions 1 and 2), which we denote by T . The ELDER number density follows the equilibrium trajectory, n eq (T ). As the temperature drops below m χ , the ELDERs become non-relativistic, and the equilibrium density drops exponentially, n eq (T ) ∝ e −mχ/T . The rates of the reactions 2, 3, and 4, drop off exponentially, while the reaction 1 slows more gradually.
All the reactions eventually decouple, Γ < ∼ H, but the order of decoupling is crucially important in determining the relic abundance. It is natural for 3 → 2 self-annihilation to decouple before 2 → 2 self-scattering: the interaction strengths entering the two rates are generically of the same order (both involve interactions internal to the dark sector), but Γ 3→2 ∝ n 2 χ while Γ 2→2 ∝ n χ . On the other hand, the rate of annihilations to SM, Γ an , is controlled by the coupling between the SM and the dark sector, which can naturally be small. (For example, in the dark photon portal model considered below, this will be controlled by kinetic mixing between the SM and dark-sector U (1) gauge groups.) In this paper, we will consider the regime where annihilations to SM decouple first, while the 3 → 2 process is still active. This is the case in both the SIMP and ELDER scenarios.
The rate of elastic scattering Γ el is proportional to the SM density, which is not exponentially suppressed at T < m χ . However, the scattering cross section is suppressed by the small coupling between the SM and χ. Generically, this cross section is of the same order as that of annihilations to SM, and therefore decoupling of elastic scattering occurs after annihilations to SM are decoupled. Depending on the relative strength of the SM-χ coupling and χ self-couplings, the decoupling of the elastic scattering may occur either after or before the decoupling of the 3 → 2 self-annihilation. The former case corresponds to the SIMP scenario [18], while the latter is the ELDER scenario [7].

The Thermal History of ELDERs
After annihilations and elastic scattering with the SM decouple, but while the 3 → 2 and 2 → 2 self-interactions are still active, the ELDERs are still in thermal equilibrium at zero chemical potential, but their temperature T no longer has to be the same as the SM plasma temperature T . As shown in Appendix B, the two temperatures are related by Here ψ is the SM particle that couples to χ, with corresponding number of degrees-of-freedom g ψ and g χ , respectively; N ψ 3+n is a numerical constant given in Eq. (B.4). We assume that the effective number of relativistic degrees of freedom g * ,d remains constant throughout the decoupling process. (The case of varying g * can be handled numerically.) The "elastic scattering strength" c n is defined as the dimensionless coefficient of the leading term in the lowenergy expansion of the matrix element-squared of the elastic scattering process χψ ↔ χψ: where |M| 2 is averaged over initial and final-state degrees of freedom, including spin, color, and electric charge. (See Appendix A for details.) If χ couples to more than one SM particle, a summation over the relevant SM species is implied in the definition of a. The formalism presented here is applicable to SM particles that are relativistic at the time of χ decoupling, m SM T d ∼ m χ /10. SM particles with m SM T d are irrelevant to the decoupling process, while the case m SM ∼ T d can be studied numerically.
An approximate analytic solution to the temperature evolution equation can be found (see Appendix B): x = e t a n + 4 At small x, x ≈ x, corresponding to SM and ELDER sectors in thermal equilibrium. At large x, the asymptotic form of the solution is x ≈ 3 log(x) + a n + 4 (2.5) Identifying the "decoupling temperature" at which the ELDER and the SM thermally decouple, Eq. (2.5) can be rewritten as x This is precisely the behavior expected in the "cannibalization" regime [1], where ELDER temperature decreases only slowly (logarithmically with the scale factor) as the universe expands. The physical reason is that the kinetic energy released by 3 → 2 self-annihilations partially compensates for the energy lost when particle momenta are redshifted due to the expansion. This regime persists until the 3 → 2 process decouples, after which the ELDER density is frozen out. Note that the dark matter particles remain non-relativistic throughout the cannibalization period, so that from the point of view of Cosmic Microwave Background (CMB) and structure formation, ELDER is a Cold Dark Matter (CDM) candidate, consistent with observations. The evolution of ELDER temperature throughout the kinetic decoupling and freeze-out process is illustrated in Fig. 1. The ELDER-to-SM temperature ratio starts growing after kinetic decoupling due to cannibalization, reaching the maximum value of T /T ∼ 10 at the time of freeze-out. It drops rapidly after freeze-out since ELDERs are non-relativistic and T ∝ R −2 , while T ∝ R −1 , where R is the size of the universe. The analytic function (2.4) provides an excellent approximation to the numerical solution of the Boltzmann equations up until 3 → 2 freezeout.
We note that Eq. (2.7) can also be derived by assuming instantaneous kinetic decoupling between the dark sector and the SM at temperature T d , and using the conservation of comoving entropy in the dark sector after decoupling. This approach was taken, for example, in Ref. [7]. The alternative derivation presented here does not make the assumption of instantaneous decoupling, relying instead on the approximate solution for the evolution of T accurate throughout the decoupling process. Apart from being better justified physically, the distinct advantage of the new derivation is that it automatically provides the expression for T d in terms of the underlying model parameters, Eq. (2.6).
In the instantaneous freeze-out approximation, the asymptotic value of the yield Y χ = n χ /s 0 , where s 0 is the entropy density today, is given by where x f and x f are the temperatures of the SM and the ELDERs, respectively, at the time of freeze-out. The effective multiplicity at freeze-out, g * s,f , is strongly dominated by the SM degrees of freedom that are relativistic at that temperature, and the ELDER contribution to entropy is negligible; for typical ELDER parameters, g * s,f = 10.75. The ELDER relic density is given by where x d is the decoupling temperature defined in Eq. (2.6). The 3 → 2 self-annihilations freeze-out when n 2 χ σ 3→2 v 2 H. Let us parametrize (2.10) The freeze-out and decoupling temperatures can then be estimated by solving the equations Numerically, x d 17 and x f 25 for a typical ELDER model. The decoupling temperature is directly related to the strength of elastic scattering between ELDERs and SM particles, see Eqs. (2.2), (2.6). Once x d is found by solving Eq. (2.12), it is straightforward to compute the corresponding elastic scattering strength: where g * ,d is the effective number of relativistic degrees of freedom at T d , and ξ n = (n + 4)[Γ( n+3 n+4 )] −n−4 /N ψ 3+n is a numerical constant. (For future reference, ξ 0 0.08 and ξ 2 0.004.) Once a mechanism that mediates ELDER-SM scattering is specified, this formula can be used to make detailed, robust phenomenological predictions, as discussed in the next Section. Remarkably, such predictions are almost completely independent of the details of self-interactions of ELDERs, or their interactions with other dark sector states.

ELDER Mass Estimates
A model-independent upper bound on the ELDER dark matter particle mass can be obtained as follows. Self-consistency of the ELDER scenario requires (2.14) Here we see that ELDER dark matter is pushed to the strongly interacting regime (α 1). The thermally averaged 3 → 2 rate can be bounded above by unitarity, in similar spirit to the bound derived on the thermally averaged WIMP annihilation rate [24]. The optical theorem states that where M forward is the matrix element for forward scattering χχ → χχ, and dΠ X is the Lorentz invariant phase space. Picking the term with X = χχχ from the sum yields the inequality Close to this bound, the kinetic decoupling and freeze-out occur close in time, and the formulas derived in this Section, which assumed a clear separation between the two events, are not strictly applicable. The bound on α for "pure ELDER" regime, in which the separation is clear, is stronger by a about a factor of 2. For smaller α, a "mixed SIMP-ELDER" regime occurs, which does not lend itself to simple analytic estimates. Numerical analysis of this regime indicates a smooth connection between "pure SIMPs" and "pure ELDERs", see for example Fig. 2. Using this in the definition of the thermally averaged rate in Eq. (A.7), in the non-relativistic limit, the rate is bounded above by In the absence of light degrees of freedom, non-relativistic elastic scattering of scalar χ particles is typically dominated by the s wave. Partial-wave unitarity requires 3 |M forward | ≤ 8π √ s/p 48π/ √ 5, which in turn implies (taking into account the typical freeze-out temperature x f 20) an upper bound α < ∼ 73. (2.18) where α is defined in Eq. (2.10). Combining this bound with Eq. (2.14) yields This partial-wave unitarity bound is independent of the details of the dark sector. In specific models of dark sector self-interactions, other considerations, such as perturbativity of couplings, may impose stronger bounds. For example, in simple scalar models discussed in Section 4, the upper bound on the ELDER mass from perturbativity is about 200 MeV. There is also a lower bound on m χ . As the ELDER becomes non-relativistic, energy and entropy are transferred from the dark sector to the SM, reheating the SM degrees of freedom. This process continues until the decoupling of elastic scattering between ELDERs and the SM at temperature T d . If the energy and entropy transfer is active during or after Big-Bang Nucleosynthesis (BBN), it will generally result in modification of BBN predictions for lightelement abundances, and/or the effective number of neutrinos N eff inferred from the Cosmic Microwave Background (CMB) measurements; see e.g. Ref. [25]. This is certainly the case if the interactions between the ELDER and the SM are mediated via the dark photon portal, which, as argued in Section 3, is the most plausible renormalizable portal compatible with this scenario. The dark photon portal couples the ELDERs very weakly to neutrinos. If entropy transfer continues below the temperature of neutrino decoupling from the electron/photon plasma, non-standard N eff is produced. It is in principle possible that this bound could be avoided in a model in which electrons, photons and neutrinos are reheated equally. However in this paper we will adopt [25] m as a rough lower bound on the ELDER mass.

ELDERs, SIMPs and WIMPs, Oh My!
If in a given model c n <c n , defined in (2.3) and (2.13), the particle χ cannot account for the observed dark matter. On the other hand, if c n >c n , the correct relic density can still be 3 At √ s = 3m, the χ particles are moderately relativistic, β 2 ∼ 0.5, and corrections to s-wave scattering amplitude may be non-negligible. This will affect the unitarity bound at the level of order-one factors. Thus, this bound as well as the mass bound in Eq. (2.19) should be viewed as order-of-magnitude estimates.

ELDER ELDER
10 -8 10 -7 10 -6 10 -5 10 -2 10 -1 1 10 1 Figure 2: Regions of parameters corresponding to the observed relic density. For each mass, the vertical section of the line of the left/top corresponds to the elastically decoupling relic (ELDER) scenario proposed in this paper; the horizontal line to the SIMP scenario; and the vertical section on the right/bottom to the WIMP scenario. This figure, reproduced from Ref. [7], is based on a numerical solution of the Boltzmann equations for a model with g χ = 2, ψ = photon, n = 0, c 0 = 8π 2 . The same behavior is observed in other models, see for example Fig. 6 below.
achieved through the SIMP mechanism. In this case, dark matter and SM remain in kinetic equilibrium until the 3 → 2 interactions decouple and the χ density freezes out: x d > x f . The relic density is given by where the freeze-out temperature x f is found as a solution to After freeze-out, elastic scattering with SM no longer affects n χ ; thus in the SIMP regime, the relic density is determined by the self-interaction strength α, and is independent of c n . The SIMP value of α, is close to the lower bound on α required for the ELDER scenario, Eq. (2.14), and scales the same way with m χ . This gives a clear intuitive picture of the relation between the two regimes: for a given dark matter particle mass, the ELDER value of c n gives the lower bound on c n for SIMPs, while α SIMP is the lower bound of α for ELDERs. If c n is increased even further, eventually a point is reached where annihilations to SM decouple after the 3 → 2 interactions. At this point, the relic density is determined by the cross section of annihilations to SM, and is once again independent of α. Since this is the mechanism that sets the relic abundance of the conventional WIMPs, we refer to it as the "WIMP regime", even though the dark matter particle mass is still well below the weak scale, and a small coupling to SM is required to obtain the correct relic density. (For theoretically motivated realizations of such a scenario, see [26].) Figure 2 illustrates the three regimes. This figure, reproduced from Ref. [7], is based on a numerical solution of the Boltzmann equations for a model with g χ = 2, ψ = photon, n = 0, c 0 = 8π 2 , which was performed in that paper. The same behavior is observed in other models, see for example Fig. 6 below.

Dark Photon Portal and Phenomenology
It is well known that there are only three renormalizable interactions that can couple SM to dark sector states: "dark photon", "Higgs", and "right-handed neutrino" portals [27]. Of these, only the dark photon portal is compatible with the ELDER scenario in its simplest form. In the case of the Higgs portal, the interaction has the form S 2 H 2 , where S is a darksector field and H is the SM Higgs. In the case of ELDER, the decoupling temperature is at the MeV scale, and the relevant SM degrees of freedom are electrons, photons, and neutrinos. The couplings to these particles at MeV temperatures mediated by the Higgs are too weak to produce the elastic scattering of the strength required in the ELDER scenario. In the case of the neutrino portal, the interaction is of the form HLN , where N is a dark-sector fermion. The ELDER dark matter particle must possess 3 → 2 interactions, and thus must be a boson. Hence, we will focus on the dark photon portal as the most plausible mechanism for ELDER-SM coupling.

Dark Photon Portal
Specifically, we consider a complex scalar field χ, neutral under SM gauge symmetries but charged under an abelian U (1) D gauge group in the dark sector: where g D is the U (1) D coupling constant, and A is the corresponding gauge field. The A kinetically mixes with the SM photon: 4 where B and F D are the field strength tensors of the U (1) Y and U (1) D , and θ W is the Weinberg angle. Diagonalizing the kinetic terms yields the SM photon A, under which χ is uncharged, and the "dark photon" V , which couples to the SM electromagnetic current with strength γ e, and to the "dark" U (1) D current with strength g D . We further assume that U (1) D is broken, giving the dark photon mass m V . (For specific models that realize this setup, including ELDER self-interactions, see Section 4.) If dark photons have a significant abundance in the early universe at the time of ELDER decoupling and freeze-out, the physics of these processes becomes considerably more complicated: for example, co-annihilation processes may play an important role in transferring energy between the SM and the dark sector. To avoid these complications, we focus our attention on the "pure ELDER" case, when the dark photon is significantly heavier than the dark matter particle. For concreteness, we assume m V > 2m χ . Elastic scattering of ELDER on electrons is mediated by the t-channel dark photon exchange. In the language of Section 2, the dark photon portal model corresponds to ψ = e ± , g ψ = 4, n = 2, and the elastic scattering strength is given by 5 Here we defined the dimensionless combination where α D = g 2 D /(4π). This is the same combination of parameters that controls dark matter annihilations to the SM, as has been previously noticed in studies of the conventional scenario where such annihilations determine the relic density [28]. In the ELDER scenario, the value of y that corresponds to the observed relic density can be inferred from Eq. (2.13): where x d is the solution to Eq. (2.12). This is a robust prediction of the ELDER scenario with the dark photon portal, independent of the details of ELDER self-interaction dynamics. As discussed above, if the dark matter coupling to the SM is increased above the ELDER value, correct relic density can still be achieved by SIMP or WIMP mechanisms. In the dark photon portal model, the WIMP regime corresponds to the well-known "thermal target" value for y [28]: where ξ = 1 − 4m 2 χ /m 2 V , and x f,a is the temperature at which annihilations to SM freeze out. Any value of y between y ELDER and the thermal target is compatible with the SIMP mechanism, which can yield the correct relic density for appropriately chosen 3 → 2 selfscattering cross sections.
Before proceeding, let us briefly comment on the astrophysical constraints on this model. Dark matter pair annihilation into electrons is constrained by the CMB measurements [29][30][31], as well as indirect-detection searches. However, in the case of scalar dark matter in  [36][37][38][39], superconductors (10 meV threshold) [40,41], superfluids [42,43], scintillators [36,44] and graphene [45]. the relevant mass range, the s-wave annihilation cross section is suppressed by a factor of (m e /m χ ) 2 < ∼ 10 −2 , while the p-wave contribution is velocity-suppressed. As a result, ELDER dark matter is easily consistent with these constraints. Also, the reaction e + e − → χχ (with or without an on-shell dark photon) can provide an additional mechanism of cooling in supernovae, which is constrained by the observation of neutrinos from SN1987A (see e.g. [32,33]). We checked that in the ELDER region, the elastic scattering of χ on electrons is always sufficiently strong to prevent the dark matter particles from leaving the supernova core. The produced χ's become trapped in the core, and do not contribute to the cooling rate.

Direct Detection
Direct detection of sub-GeV dark matter has been an area of active recent investigations. Heavy nuclear recoils do not carry sufficient energy to be detected in this mass range, and direct detection is easier for dark matter scattering on electrons. Remarkably, in the ELDER scenario with a dark photon portal, it is precisely the same process that determines the DM relic density. The observed dark matter density completely determines the direct detection cross section, with essentially no free parameters other than the ELDER mass m χ . The direct detection cross section is given by Setting y = y ELDER in this formula defines a very sharp "ELDER target" for the direct detection experiments. This complements the "thermal target" [27,28], which in our language corresponds to y = y WIMP , while the region y ELDER < y < y WIMP corresponds to SIMP dark matter. Moreover, as discussed above, observational constraints and unitarity considerations restrict m χ to a range between roughly 5 MeV and 1 GeV. These considerations define the direct detection target region, shown in Fig. 3.
The predicted cross sections are well below the current XENON bounds [34,35]. However, novel experimental approaches that are currently being investigated have the potential to dramatically increase the sensitivity to DM-electron scattering in this mass range. Target materials under study include semiconductors [36][37][38][39], noble liquids [34,36], superconductors [40,41], superfluids [42,43], scintillators [36,44] and graphene [45]. Projected sensitivities of these experiments will allow them to test a significant part of the SIMP and ELDER target region, see Fig. 3.

Dark Photon Searches
Searches for a dark photon in the MeV-GeV range have also been an area of much activity recently. Existing experimental data has been used to place bounds on the dark photon, and several dedicated experiments are now running or in preparation. The ELDER, SIMP and WIMP scenarios with dark photon portal provide a well-defined dark photon target region for such experiments, shown in Fig. 4.
In the ELDER scenario, the dark photon mass m V must be large enough so that the process χχ * ↔ V V is not relevant throughout the χ kinetic decoupling and freeze-out process. For the discussion of this section, we assume m V > 2m χ . In this case, the decay V → χχ * is likely to be the dominant dark photon decay channel, since its amplitude is proportional to the dark sector gauge coupling g D , which is naturally of order one, while the amplitudes of competing decays such as V → e + e − are controlled by the small kinetic mixing parameter γ . As a result, the experiments relevant for constraining our scenario are those searching for invisible dark photon decays. There are two basic experimental approaches. First, one can search for missing mass or energy in collider events due to an invisible particle V . The strongest current constraints from this approach come from re-analysis of BaBar data [46], as well as, at low masses, the dedicated NA-64 experiment at CERN [47]. These searches do not yet constrain the ELDER scenario. In the future, the missing-energy LDMX experiment proposed at SLAC [28,48] will have sufficient sensitivity to test a significant part of the ELDER parameter space. Second, one can search for a dark matter particle that is produced  Figure 4: The dark photon target region predicted in the ELDER, WIMP and SIMP scenarios. For comparison, the current bounds and projected sensitivities of searches for dark photon decaying to dark matter particles [27] are also shown.
in dark photon decay and propagates through shielding material to a downstream detector. (This would in effect amount to "direct detection" of a dark matter particle produced in an accelerator.) This approach was recently pioneered by the MiniBooNE experiment [49], and dedicated experiments such as BDX [50] and SHiP [51] have been proposed. Such future experiments may be sensitive to ELDER and SIMP dark matter. A snapshot of the current and expected sensitivities of a variety of dark photon searches, collected in Ref. [27], and overlaid with the ELDER and SIMP regimes, is shown in Fig. 4.
We remind the reader that while the theoretical predictions of the dark photon target region are naturally defined in terms of the y variable, and are largely insensitive to variations of model parameters that leave y unchanged, the same is not true of experimental sensitivities, which depend on model parameters in different ways. For example, sensitivity of a missingmass experiment such as BaBar is completely independent of g D , as long as it's large enough so that the invisible branching ratio of V is close to 100%. Thus, additional assumptions have to be made in displaying experimental sensitivities in terms of y, as in Fig. 4; see Ref. [27] for further discussion.

Models of ELDERs
We argued above that strong self-interactions in the dark sector are required in the ELDER scenario, with 3 → 2 cross-section of order one in its natural units. At some level, this is welcome: Such strong self-interactions are indeed expected if the ELDER is a bound state of confining dynamics in the dark sector, a paradigm that can potentially provide a natural explanation of proximity of the ELDER mass to the QCD confinement scale. On the other hand, it does create obvious challenges for model-building. Moreover, strong numberchanging self-interactions tend to be accompanied by a large ELDER elastic scattering cross section, which can run afoul of observational constraints on dark matter self-scattering in galactic clusters such as the Bullet cluster. Fortunately, many phenomenological predictions of the ELDER scenario are independent of the details of dark sector self-interactions. This allowed us to completely sidestep these questions in the discussion of Section 3. We will now discuss two simple, renormalizable dark-sector models that explicitly realize the ELDER scenario. While not deeply rooted in strong gauge dynamics, they can be thought of as toy models representing interactions among the lightest mesons produced by such dynamics. They provide a useful illustration of the issues involved in dark-sector model building, and an "existence proof" demonstration that consistent models can be found.

χ 3 Model
Here we consider a simple model in which the dark matter is a complex scalar charged under an unbroken Z 3 symmetry [3,18,52]. Consider a dark sector consisting of a U (1) D gauge field with gauge coupling g D , and two scalar fields charged under it, Φ and χ, with Q(Φ) = +3 and Q(χ) = +1. The χ particle will play the role of dark matter. The scalar potential is where V (ψ) = m 2 ψ |ψ| 2 + λ ψ |ψ| 4 . We will assume m 2 Φ < 0, so that this field gets a vacuum expectation value (vev) Φ = w/ √ 2. We further assume that m 2 χ is positive. For simplicity, we consider the situation m χ < |m Φ |, with sufficient separation to ensure that the radial degree of freedom of Φ is sufficiently heavy to not play a role in the calculation of χ relic abundance. The effective Lagrangian for such calculation is then given by where we defined a dimensionless 3-point coupling The vev of Φ leaves a global Z 3 subgroup of the U (1) D unbroken, and the charge of χ under this discrete symmetry guarantees its stability, as required for a dark matter candidate. The U (1) D gauge boson gets a mass m V = √ 3g D w. The symmetry of the theory allows for kinetic mixing between the U (1) D gauge boson and the SM hypercharge gauge boson, as in Eq. (3.2). As long as there are states, at any mass scale, that are charged under both gauge groups, such kinetic mixing will be generated, with values of γ ∼ 10 −4 − 10 −2 being generic if no cancellations occur at the one-loop level [53]. Thus, this construction provides a stable scalar dark matter candidate with natural coupling to the electron via a dark photon portal. The matrix elements for non-relativistic 3χ → 2χ annihilations are given by Here we set λ χ = 0 for simplicity. This point is unexceptional (there is no enhanced symmetry associated with vanishing of λ χ ) and is sufficient to illustrate the important physical features of the model. This yields the thermally averaged cross section In the SIMP scenario, the coupling R can be inferred from the relic density as follows: The required coupling is quite large, consistent with the idea that SIMP/ELDER dark matter particle can be a bound state of dark-sector confining gauge group: in this scenario, the potential (4.1) can be thought of as a toy model representing the interactions among the two lightest mesons. The range of validity of the perturbative χ 3 model can be estimated as R < ∼ 4π. In the SIMP scenario, this gives an upper bound on the dark matter particle mass: As discussed in Section 2, the "pure ELDER" scenario requires larger 3 → 2 cross section than SIMP for the same m χ , and therefore the upper bound on m χ is somewhat lower for ELDERs. The dark matter elastic self-scattering cross section is constrained by observations of galactic clusters, such as the Bullet cluster [54][55][56], and halo shapes [57][58][59]: while in the ELDER scenario the cross-section is even larger (bounded from below by Eq. (4.9)). Thus, the simplest single-field χ 3 model cannot provide sufficiently strong self-interactions required in these scenarios, while being consistent with observational constraints. We will now show that adding another dark-sector scalar field can resolve this problem.

Choi-Lee Model
This model was originally introduced by Choi and Lee (CL) [9] in the context of the SIMP scenario. The dark sector contains a U (1) D gauge symmetry, with gauge coupling g D , and three complex scalar fields charged under this symmetry: φ, S, and χ, with charges q φ = +5, q S = +3, and q χ = +1. The most general renormalizable scalar potential consistent with these charge assignments is We assume that m 2 φ < 0, while the other two scalar fields have positive mass-squared. The vev Φ = w/ √ 2 breaks the gauge symmetry, giving the U (1) D gauge boson a mass m V = √ 5g D w. The Φ vev preserves a discrete Z 5 subgroup of the U (1) D , under which S and χ are both charged. The lighter of these particles, which we will assume to be the χ, is therefore stable, and can play the role of dark matter. The scalar interactions after spontaneous symmetry breaking are described by where we have omitted interactions with the Higgs component of Φ which play no role in the phenomenology considered here, and defined dimensionless couplings As in the χ 3 model, the dark gauge boson V kinetically mixes with the SM photon, providing a dark photon coupling between the dark sector and the SM. The 2χ ↔ 3χ scattering process is induced by the couplings in the first line of Eq. (4.11). For simplicity, we set λ 3 = 0; this point is unexceptional (there is no enhanced symmetry associated with vanishing of λ 3 ) and is sufficient to illustrate the features of interest to us. The key observation is that for m S ≈ 3m χ , the 2χ ↔ 3χ scattering is resonantly enhanced, while the 2χ ↔ 2χ process is not. This effect is illustrated by the left panel of Fig. 5, where we plot the thermally-averaged σ 3→2 v 2 at temperature close to ELDER kinetic decoupling. A dimensionless ratio of the number-changing and number-preserving cross sections, m 2 2→2 , can reach O(10 3 ). For comparison, in the χ 3 model studied in the previous section, this ratio is close to 1. Note that the values of couplings R i required in the SIMP/ELDER scenarios are fairly large, so that the S resonance is rather broad and no significant fine-tuning of m S /m χ is required to achieve significant enhancement of the 3 → 2 rate. This enhancement makes it possible to successfully implement SIMP and ELDER dark matter in the CL model without conflict with observational constraints from galaxy clusters and halo shapes.
Because of the resonance at √ s ≈ 3m χ , the quantity σ 3→2 v 2 has a non-trivial temperature dependence in the non-relativistic regime, making the parametrization of Eq. (2.10) inapplicable. To compute the relic density, we integrate the Boltzmann equations numerically. The relic density is controlled by the seven model parameters that enter the Boltzmann equations: particle masses m χ , m S , and m V ; and dimensionless coupling constants R 1 , R 2 , g D , and γ . To perform numerical analysis in this large parameter space, we made the following choices: • The ratio of S and χ masses was fixed close to the 3 → 2 resonance, m S /m χ = 3.1.
• As discussed in Section 3, the relic density depends on the three parameters of the dark photon portal only through a single dimensionless combination y, defined in Eq. (3.4). Therefore it is sufficient to fix two of these parameters, and vary the third one. We fix m V /m χ = 10 and g D = 1, and vary γ .
• The 3 → 2 matrix element is proportional to a product R 1 R 2 2 , so that the DM relic density primarily depends on these couplings through the "effective" 3 → 2 coupling, The relic density also depends on the width Γ S , which is proportional to R 2 2 . In practice, in the numerical analysis we fix R 2 (specifically, R 2 = 2 for m χ = 10 MeV and R 2 = 4 for m χ = 35, 100 MeV) and vary R 1 . However, we checked that in the parameter range of interest, the relic density is insensitive to variations of R 2 within broad ranges around these values, allowing us to present the results solely in terms of the effective coupling R eff .  Figure 6: Left: Regions of CL model parameter space with χ relic density consistent with the current best-fit ΛCDM value. Right: Constraints from galaxy cluster observations (the regions below the curves are allowed). In both plots, m S /m χ = 3.1, m V /m χ = 10, and g D = 1. In the right panel, we fix R 1 = 10 for illustration.
The three remaining parameters (m χ , R eff , γ ), are scanned over. The results are illustrated in the left panel of Fig. 6, which shows regions of parameter space consistent with the current best-fit ΛCDM dark matter density, Ω χ h 2 = 0.1188±0.0010 [60], in the γ −R eff plane for three values of DM mass, 10, 35 and 100 MeV. The (roughly) horizontal bands of viable parameter space correspond to the SIMP scenario, while the (roughly) vertical bands realize the ELDER scenario. The values of γ for the ELDER regime are in excellent agreement with the results of the analytic approach, Eq. (3.5). In the intermediate regime, the DM-SM elastic scattering and the DM number-changing self-scattering decouple at roughly the same time, and both processes play a role in determining the relic density.
As expected, realizing ELDER (or SIMP) dark matter in the CL model requires O(1) couplings among the scalars of the dark sector. The range of validity of perturbative CL model can be estimated as R 1 < ∼ 4π, R 2 < ∼ 4π. Combined with the relic density calculation, these constraints place an upper bound on the ELDER dark matter mass, m χ < ∼ 200 MeV. Furthermore, the χ self-scattering cross section is constrained by observations of galactic clusters and halo shapes, Eq. (4.8). The self-scattering cross section receives contributions from a quartic coupling λ χ as well as the S-exchange diagram controlled by R 2 , and partial cancellation of the two diagrams is possible. Combined with the perturbativity bound on R 1 , cluster observations place an upper bound on R eff , shown in Fig. 6. For m χ > 5 MeV, the values of R eff required in the ELDER scenario are compatible with observations.

Conclusions
In this paper, we studied the Elastically Decoupling Relic (ELDER) scenario for thermal dark matter. We presented an approximate analytic solution for the evolution of ELDER temperature throughout the kinetic decoupling epoch. This solution was used to provide explicit formulas relating various relevant quantities, such as, for example, the relic density of ELDERs and the cross section of their elastic scattering off SM particles. We also applied partial-wave unitarity constraint to obtain a bound on the allowed mass range for the ELDER dark matter candidate, 5 MeV < ∼ m χ < ∼ 1 GeV. These results are valid in a broadly modelindependent framework.
Further, we showed that a dark photon portal can naturally provide the coupling between the dark matter particles and SM of the strength required in the ELDER scenario. Within the dark photon model, the ELDER scenario provides unambiguous predictions for dark matter direct detection experiments and dark photon searches, shown in Figs. 3 and 4. These predictions have no free parameters other than the dark matter mass. They are also independent of the details of dark sector, as long as it provides sufficiently strong number-changing self-interactions to realize the ELDER scenario. Together with the well-known "thermal target" and predictions of the Strongly-Interacting Massive Particle (SIMP) scenario, the ELDER predictions delineate a well-defined target region in the parameter spaces relevant for direct-detection and dark photon searches, which will be explored by the next generation of experiments.
Both the ELDER and SIMP scenarios require O(1) strength (in natural units) of selfinteractions among dark matter particles. Here, we studied two simple scalar-field models that incorporate such interactions while remaining within perturbative regime, for m χ < ∼ 200 MeV. The models also naturally contain coupling to the SM via the dark photon portal. The simplest model, with just two scalar fields, exhibits tension with bounds on dark matter selfscattering cross section from observations of galaxy clusters. However, a slightly more complex model, with three scalar fields and a resonance structure, easily evades such bounds. These results indicate that there is no fundamental obstruction to finding dark sectors compatible with ELDER and/or SIMP scenarios.
An important motivation for SIMP and ELDER scenarios is the proximity of the predicted dark matter particle mass to Λ QCD , a well-established important scale in the SM. In the toy models studied in this paper, m χ ∼ Λ QCD is put in by hand. The natural next step in the model-building direction would be to construct models in which this relation, as well as the strong self interactions, emerge naturally from UV physics. tion (MP), and the Bethe Postdoctoral Fellowship (EK). YT is partially supported by the Visiting Graduate Fellow program at Perimeter Institute.

A Boltzmann Equations
The Boltzmann equation for the DM phase space distribution, f χ (p; t), in an expanding Universe is is the Hubble expansion rate, and C[f χ ] is the collision term. For ELDER dark matter, the relevant collision terms are the 3 → 2 self-annihilations and χ-ψ elastic scattering (ψ can be any light SM particle). The collision terms also includes annihilations to SM, but their effect in the ELDER scenario is negligible, and are omitted. Strong elastic self-scattering of ELDERs ensures that, throughout the kinetic decoupling and freeze-out process, the phase space distribution follows a thermal distribution: where µ χ (t) is the chemical potential, and T (t) is the temperature of the dark sector. Eq. (A.1) can be most easily solved by taking the first two moments, the DM number density n and energy density ρ: where g χ is the number of degrees of freedom in χ. These obey where n eq χ is the density of χ in chemical equilibrium (i.e. at zero chemical potential). The thermally averaged 3 → 2 annihilation and energy transfer rates are where is the Lorentz invariant phase-space integration volume. The squared matrix elements, |M| 2 , are averaged over initial and final degrees of freedom, including spin, color, and charge. 6 During the kinetic decoupling and freeze-out process, the χ particles, to a good approximation, follow a Maxwell-Boltzmann distribution. Then, Here 'eq' denotes the values of the variables in chemical equilibrium, µ χ = 0: and n χ /n eq χ = e −µχ/T . The Boltzmann equations (A.4), (A.5) then reduce to a system of coupled partial differential equations for T and µ χ (or equivalently T and n χ ).
In the epoch of interest, the entropy of the universe is dominated by relativistic SM degrees of freedom: s 0 = 2π 2 45 g * s T 3 , where T is the SM plasma temperature, and g * S = 10.75 at the relevant temperatures (0.5 MeV < ∼ T < ∼ 100 MeV). The contribution of ELDERs to the entropy is suppressed both because they are non-relativistic, and because the number of degrees of freedom is small compared to SM. Neglecting this contribution, the time variable in the Boltzmann equations can be conveniently traded for the SM temperature: (A.13)

B Kinetic Decoupling and Approximate Analytic Solution
Since the relic density of ELDER dark matter is primarily determined at the time of its kinetic decoupling from the SM, we would like to obtain analytic insight into this process. Kinetic 6 For the case of complex χ considered in this paper, we treat χ and χ * as two states of the same particle, and averaging over these two states for each initial and final dark matter particle. For instance, for self-scattering and for χ ( * ) e ± → χ ( * ) e ± in the dark photon portal, decoupling occurs before freeze-out of the 3 → 2 interactions, so that µ χ = 0 throughout the decoupling process. The ELDERs can then be completely characterized by their temperature T , whose evolution is dictated by Eq. (A.5). In this section, we describe an approximate analytic solution to this equation, which in turn yields an analytic estimate of the ELDER relic density.
In the limit that the non-relativistic χ particles are scattering off thermalized relativistic ψ particles, an approximate analytic form of the energy transfer rate integral, Eq. (A.7), can be found. In Ref. [61], this was achieved by expanding the integrand in small momentum transfer. Here we present the necessary equations, but refer the reader to detailed calculation in the Appendix of [61]. First the thermally averaged energy transfer rate is written in terms of the collision operator in the non-relativistic limit: where c n is the leading coefficient of the matrix element expanded in E ψ /m χ at zero momentum transfer and Taking f χ to be the Maxwell-Boltzmann distribution at temperature T and integrating over the collision operator yields n χ n ψ σ el vδE n χ c n g 2 ψ g χ m χ N ψ Note that when the two sectors have the same temperature, the energy transfer vanishes, which is expected for particles in thermal equilibrium. The energy transfer rate can be related to the more commonly used quantity σ el v , as follows: When the 3 → 2 process is active and the dark matter particles are non-relativistic, they follow equilibrium Maxwell-Boltzmann distributions, and the energy density Boltzmann equation (A.5) gives a differential equation for the temperature where a ≡ c n g 2 ψ g χ m χ N ψ 3+n 32π 3 H T =mχ . (B.8) On the right-hand side, there are two competing terms. The first term contributes to the cannibalization of the dark matter, which tends to increase the dark temperature relative to the SM. The second term, which comes from the elastic scattering term, pushes T → T . The scattering term falls faster with temperature, and at some point will no longer be able to compete. At that time, the dark matter will thermally decouple from the SM bath, and cannibalization will take over the evolution of the dark sector. This decoupling occurs roughly when the second term is of order one: 7 T d m χ a −1/(4+n) . (B.9) After decoupling, the second term can be ignored, and dark temperature grows only logarithmically relative to the SM temperature, until the dark matter freezes out. We can attempt to find the analytic asymptotic behavior of the Boltzmann equations. Recasting Eq. (B.7) in terms of x and x yields The kinetic decoupling temperature can also be estimated by observing that in equilibrium, the rate of energy transfer to the SM must keep up with the rate of kinetic energy release by 3 → 2 annihilations: ne σ el vδE ∼ −m 2 χ HT −1 . According to Eq. (B.6), σ el vδE ∼ σ el v T 2 /mχ. This approach, which was used in Ref. [7], gives a result consistent with Eq. (2.6).
There appears to to be no closed form solution to the above differential equation, but the following differential equation does have a closed form solution: (B.12) In the limit x x d , x = x so the two differential equations are approximately the same. Likewise, when x x d , the 2nd term is negligible in both equations, so the change is not relevant. The closed form solution to Eq. (B.12) is x = e t a n + 4 x (x → 0) = x , (B.14) x (x → ∞) = 3 log(x) + a n + 4 This agrees with the rough estimate of Eq. (B.9), and provides the precise value of the numerical coefficient. Assuming instantaneous freeze-out of the 3 → 2 annihilations at the dark sector temperature x f (corresponding to SM plasma temperature x f ), the dark matter yield is given by Here, the exponential dependence of relic density on the elastic scattering rate is manifest. We also note a logarithmic dependence on the temperature at freeze-out, which shows only a very minor dependence on the 3 → 2 rate, provided it is still active at decoupling.
If the matrix element has significant dependence on kinematics even in the non-relativistic regime, the remaining integrals have to be evaluated numerically. This is the case in the Choi-Lee model of Section 4.2, where a resonance at √ s ≈ 3m χ can lead to rapid change of the matrix element with s near threshold. Our analysis of that model is therefore based on numerical evaluation of Eq. (C.6). In most cases, however, the matrix element in the non-relativistic regime can be approximated as a constant, independent of kinematics. In this case, all integrals in Eq. (C.6) can be evaluated analytically. This yields