Super heavy thermal dark matter

We propose a mechanism of elementary thermal dark matter with mass up to $10^{14}$ GeV, within a standard cosmological history, whose relic abundance is determined solely by its interactions with the Standard Model, without violating the perturbative unitarity bound. The dark matter consists of many nearly degenerate particles which scatter with the Standard Model bath in a nearest-neighbor chain, and maintain chemical equilibrium with the Standard Model bath by in-equilibrium decays and inverse decays. The phenomenology includes super heavy elementary dark matter and heavy relics that decay at various epochs in the cosmological history, with implications for CMB, structure formation and cosmic ray experiments.


INTRODUCTION
One of the biggest questions in fundamental physics is the nature of dark matter (DM). The possibility that DM is a thermal Weakly Interacting Massive Particle (WIMP), whose abundance is determined by 2 → 2 annihilations into Standard Model (SM) bath particles, is exciting, but has alluded detection thus far. The WIMP is particularly intriguing because it is very predictiveits abundance is determined only by its interactions with the SM, which informs us how it may be detected.
The WIMP paradigm has been a guide towards the properties of DM, such as its mass and interactions. In particular, within the WIMP freezeout mechanism, there is a an upper bound on the DM mass from perturbative unitarity, of O(100) TeV [1]. The reason is that the WIMP annihilation rate is proportional to an exponentially decreasing DM density, and so the amount of dark matter that can be annihilated away before freezeout is limited by the theoretical size of the cross-section.
In this Letter we present a new freezeout mechanism within a standard cosmological history. The DM is an elementary particle that is thermalized with the SM at high temperatures, its relic abundance is determined via its freezeout from the SM bath, and the DM mass can be as high as 10 14 GeV for s-wave processes, without violating the perturbative unitarity limit. In future work, we show how Planck-scale DM can be reached for velocity-dependent processes [24].
The general idea is as follows. The dark matter consists of N approximately degenerate states, χ i (i = 1..N ). These states co-scatter [25] off of the SM bath, but only in a chain of nearest-neighbor interactions while the N th state co-decays [26][27][28][29][30][31] in equilibrium with the SM, Here χ 1 is the DM candidate. This setup is summarized in Fig. 1. The processes in Eq. (1) can maintain chemical equilibrium much longer than annihilations can, because the interaction rate for the scattering Γ sctr = n SM σv sctr never becomes exponentially suppressed. The in-equilibrium decay allows for the whole system to have vanishing chemical potential for a long time-if the N th particle was annihilating with the bath, the system would inherit the unitarity bound from coannihilations. Finally, the chain, which will typically require N > 5 − 20, depending on the DM mass, ensures the stability of the DM χ 1 .
GENERAL IDEA Consider a DM particle χ 1 , whose density changes via scattering with a light SM bath particle, where χ 2 has similar mass to χ 1 . This process can maintain chemical equilibrium much longer than annihilations with the same interaction strength, because the interaction rate for the scattering does not depend on the DM density, and therefore the rate does not become exponentially suppressed. The DM is able to maintain equilibrium to smaller temperatures, becoming more Boltzmann Freezeout mechanism: the dark matter consists of many nearly degenerate particles which scatter with the Standard Model bath in a nearest-neighbor chain, and maintain chemical equilibrium with the Standard Model bath by in-equilibrium decays and inverse decays.
suppressed than a WIMP for the same size cross-section.
However, χ 1 can only reduce exponentially for as long as χ 2 reduces exponentially (for instance, by maintaining chemical equilibrium with the SM bath). Thus in order to go beyond the unitarity bound on annihilations, some other process is still needed to reduce the χ 2 density. If this process is annihilations (such as χ 2 +χ 2 ↔ sm+sm), then it will freezeout too early, and the unitarity bound on annihilations will apply again. If instead, χ 2 decays equilibrium with the SM bath, chemical equilibrium can be maintained for much longer. Thus, our proposed mechanism is a combination of co-decay and co-scattering dark matter.
However, one can easily see that the combination of the scattering process Eq. (3), and the in-equilibrium decay process χ 2 → sm + sm, will necessarily lead to the fast decay of χ 1 via an off-shell χ 2 . In the presence of a coscattering chain, such that scatters take place only for nearest neighbors which codecay in equilibrium with χ N , χ 1 can instead be long lived; here χ 1 is still the DM, but its decay width is suppressed due to the large phase space needed to decay to the SM (χ 1 decays to 2N SM particles). This is summarized in Fig. 1.

TWO PARTICLE CASE
As a toy example, we first work through the simplest case of N = 2 dark matter particles with mass m. While this example is not cosmologically viable, it demonstrates the basic idea of the mechanism. Consider two degenerate states χ 1 and χ 2 , which co-scatter [25] off of SM bath particles via the process The dark matter candidate is χ 1 , while χ 2 decays in equilibrium with the SM bath.
The Boltzmann equations for the system arė n 1 + 3Hn 1 = n sm σv (n 2 − n 1 ), where σv is a thermally averaged cross section for the scattering process, n eq is the equilibrium number density of χ, and Γ 2 is the decay rate of χ 2 in the thermal bath. This system is similar to the co-scattering scenario [25], but here the DM bath is kept in equilibrium with the SM bath via decays and inverse-decays, rather than annihilations. Ultimately, since decays and inverse-decays can stay in equilibrium much longer, this will allow for much heavier DM. Unlike the case for freezeout via annihilations, the instantaneous freezeout approximation will not give a good estimate of the relic abundance. This is because the rate for χ 1 scattering, Γ ∼ n sm σv , is not dropping off exponentially fast with the expansion, and therefore freezeout takes a long time. However, an approximate analytic solution to the relic abundance can still be determined from the Boltzmann equations, as we now detail.
If the decay rate Γ 2 is larger than the Hubble expansion parameter when the temperature T of the universe is equal to the DM mass, the number density of χ 2 closely follows its equilibrium value. An approximate solution to the χ 1 density can then be found by considering the single equation where Y i = n i /s with s the entropy density, x = m/T and λ = (n sm σv /H)| x=1 . We have also assumed that the thermally averaged cross-section is velocity independent and took the number of relativistic degrees of freedom g = g s to be constant. The asymptotic value of the relic abundance is where g χ is the number of internal degrees of freedom of a χ particle. From the Boltzmann equation (8), χ 1 departs equilibrium when λ/x 2 fo ∼ 1. The fact that the relic abundance then scales as Y 1 (∞) ∼ e −2x fo and not as e −x fo (as one finds in the instantaneous freezeout approximation), is the result of the slow freezeout. At this point χ 1 creation stops, but χ 1 can continue to scatter away. Neglecting the inverse process, solving x fo → e −2x fo . One also sees from this solution that the abundance stops changing significantly when x fin ∼ λ = x 2 fo . From the above estimation, one can find σv that reproduces the observed relic abundance of dark matter, while satisfying the unitarity bound. For instance, parameterizing the cross-section as σv ≡ α 2 /m 2 χ , one finds that for α 1, m χ = 6 × 10 14 GeV reproduces the correct relic abundance (assuming the DM scatters off of 4 SM degrees of freedom and that g s = 106.75).

N > 2 DEGENERATE CASE
Although the simplest N = 2 example does not work, it indicates a clear path forward to a viable setup: suppressing the decay rate of χ 1 . We thus consider N > 2 particles with the same type of nearest-neighbor interactions as before. Similarly, we assume that only χ N is able to decay into SM particles, that the masses are degenerate, and that the cross section for each interaction is the same. The equations of the system are where the second equation is valid for j = 2, · · · , N − 1, λ = (n sm σv /H)| x=1 , and λ d = (Γ N /H) x=1 . Here prime denotes a derivative with respect to x. We also ignore the time variation of g and g s , and assume that the cross section is constant in the non-relativistic limit. Note also that we are interested in parameter space where λ, λ d 1.
This system can be solved in similar fashion to the N = 2 case, by assuming Y N = Y eq and diagonalizing the differential equations. However, we can more simply obtain a solution with the new N scaling by taking the large N limit, and taking flavor space to be continuous. For N 1, the system behaves as a random walk between χ j states, with a reflecting wall at j = 1 and an absorbing cliff at j = N (corresponding to the decay). Taking the indices of j as a continuous variable = πj/[2(N − 1)], the system is described by the diffusion equation where τ = −1/x and D = π 2 λ/[4(N − 1) 2 ] is a diffusion coefficient, and the boundary conditions are From the form of the diffusion equation, we see that the system depends on the cross-section only through the combination D π 2 λ/4N 2 . From solving the boundary conditions, the asymptotic profile Y (∞) ∝ cos . The solution of the diffusion equations gives the relic abundance where Y ∞ is defined in Eq. (9). In Fig. 2 we plot the yield of the χ i particles for N = 10, solving the full Boltzmann equations, Eqs. (11). We find that the relic abundances match the analytical estimate of Eq. (14). Now, one must check that χ 1 is sufficiently long lived to be a viable DM candidate. A model-independent bound on the DM lifetime is τ > 5 × 10 18 sec [32] (for earlier studies, see Refs. [33][34][35][36][37][38][39][40][41][42][43][44]). However, if the DM decays to SM particles, constraints on the lifetime may be much stronger, τ > 10 27 sec (see e.g. Refs. [45][46][47]).
The decay of χ 1 to SM particles takes place only through N −1 off-shell χ i particles, into 2N SM particles, with decay width Treating the SM particles as massless, the total 2N -body phase space is [48] where S is a symmetry factor accounting for identical particles in the final state. We may approximate the squared matrix element of the decay as [49] where the factor S 2 accounts for expected number of diagrams. The decay rate of χ 1 is then estimated as Assuming the χ N decay is in equilibrium Γ N > H| x=1 , requiring the stability of the DM, τ > 10 27 sec, places a lower bound on the number of dark matter particles N . For instance, assuming that the DM decays to 2N identical SM particles corresponding to S = (2N )!, we find for m = 10 7 GeV that N ≥ 5, and for m = 10 13 GeV that N ≥ 14. In Fig. 3, we plot the minimum value of N as a function of m χ , satisfying the lifetime requirement.

NON-DEGENERATE MASSES AND COUPLINGS
Up until now, we have considered the case of exactly degenerate masses, but it is natural to consider small mass splittings between the particles. This will have two important effects: creating forbidden channels and allowing the heavier states to decay directly into the lighter states.
In standard two-particle forbidden annihilations, the mass splitting is negligible if it is smaller than the temperature when the abundance is still evolving, i.e., δm/m < x −1 fo . For the chain interactions, a similar condition is found: where m max and m min is the maximum and minimum particle mass in the chain. If the mass splitting is larger than this, then the evolution will enter a forbidden regime. The abundance of the lightest particle (the DM candidate) decouples when the Boltzmann suppression factor due to the mass splitting becomes significant and the chain reaction departs from chemical equilibrium. After this point, the abundance of the heavier particles continues its Boltzmann suppression, Y heavy ∼ e −∆m/T . Consequently, the relic abundances of the heavier states are much smaller than that of the lightest state when in the forbidden regime.
We leave a detailed study of this forbidden regime to future work.
The second important effect of the mass splitting is that the heavier states can quickly cascade decay to the lightest state (which is still very long lived). Depending on the size of the mass-splitting, the decays may be in-or out-of-equilibrium. If the decays are in-equilibrium, the calculation of the relic-abundance will mostly remain unchanged. If the decays are out of equilibrium, the heavier states will simply transfer their abundances to the lightest state, increasing the DM abundance by a factor O(N ).
It is also reasonable to consider that the interaction strengths may vary across the chain. In this case, the Boltzmann equations are best solved by diagonalizing the system of differential equations. The relic abundance is approximately given by the degenerate case, but with the cross-section replaced by the smallest eigenvalue crosssection in the chain.
Depending on the parameters of the theory (especially the size of the mass splittings), some of the χ i may decay at various cosmological epochs. For instance, a component of the DM could dissociate light elements during Big Bang Nucleosynthesis (see Ref. [66] and references therein). Additionally, DM decay can lead to spectral distortions in the CMB [67][68][69], which would be probed by the proposed PiXiE experiment [70]. Partially decaying DM has also been shown to alleviate several cosmological tensions, such as the Hubble tension [71,72] and small scale structure puzzles (see e.g. Refs. [73,74]). For further studies of multicomponent dark matter with varying lifetimes, see Refs. [75,76]. A detailed exploration of the phenomenology of our framework will be presented in upcoming work [77].

A TOY MODEL
Having described the general mechanism, we now present a simple model realization. Consider the Lagrangian where χ i (i = 1, · · · , N ) are left-handed Weyl fermions, S is a SM singlet scalar field, L is the lepton doublet, and H is the SM Higgs. The model is technically natural, since it respects a Z N 2 symmetry that is broken by δm. Any correction to off-diagonal masses must then depend on some power of δm. For simplicity, we take δm to be flavor independent; relaxing the assumption does not qualitatively alter the results.
The mass matrix is given by The above structure is similar to a tight-binding model along a one-dimensional wire in a quantum mechanical system. It is well-known that the wavefunctions in such system are localized at each site [78]. Assuming ∆m ≡ (m max − m min ) δm, the localization length is [79,80] ξ −1 loc ln ∆m 2δm − 1.
Due to the localization, the mass eigenstate ψ i can be approximated as In the mass basis, the generated nearest neighbor interactions L ⊃ ye −1/ξ loc Sψ i ψ i+1 , , .
and non-nearest neighbor interactions are exponentially suppressed.
The mixing can generate the direct decay ψ 1 → H + L, with width Taking Γ N = H(m), we find ξ −1 loc 1 N 62 + log m 10 10 GeV (26) to ensure that τ 10 27 sec. The ratio between ∆m and ∆m controls both the relic abundance and the lifetime of DM in this toy model. For the longevity of DM, a small localization length is preferred, while, for the correct relic abundance, the localization should be not so small as it could suppress the nearest neighbor interaction. In Fig. 3, we plot the minimum number of heavy particles χ i needed for the stability and the observed abundance of dark matter.

SUMMARY
In this Letter, we presented a new freezeout mechanism for super heavy DM that freezes out with the SM within a standard cosmological history. The relic abundance is determined solely via its interactions with the SM. For a velocity-independent cross-section, we showed the DM mass could be as large as m ∼ 10 14 GeV within the perturbative unitary limit. In an upcoming paper we show how velocity-dependent interactions, such as if the scattering was mediated by a light mediator, allows for DM to be as heavy as the Planck scale [24].