Roadmap to Thermal Dark Matter Beyond the WIMP Unitarity Bound

We study the general properties of the freezeout of a thermal relic. We give analytic estimates of the relic abundance for an arbitrary freezeout process, showing when instantaneous freezeout is appropriate and how it can be corrected when freezeout is slow. This is used to generalize the relationship between the dark mater mass and coupling that matches the observed abundance. The result encompasses well-studied particular examples, such as WIMPs, SIMPs, coannihilation, coscattering, inverse decays, and forbidden channels, and generalizes beyond them. In turn, this gives an approximate perturbative unitarity bound on the dark matter mass for an arbitrary thermal freezeout process. We show that going beyond the maximal masses allowed for freezeout via dark matter self-annihilations (WIMP-like, $m_{\rm DM}\gg\mathcal{O}(100~\rm TeV)$) predicts that there are nearly degenerate states with the dark matter and that the dark matter is generically metastable. We show how freezeout of a thermal relic may allow for dark matter masses up to the Planck scale.


I. INTRODUCTION
Overwhelming evidence, accumulating from the beginning of the last century, suggest the existence of dark matter (DM), yet we still don't know much about it.All we know is inferred from its gravitational interactions, while its particle nature remains a mystery: we still don't know its mass, forces it interacts with, and most importantly -how we may discover it.
We expect that the current DM abundance can be related to its microscopic properties such as mass and interaction cross-section, analogous for instance, to how the abundances of the light nuclear elements are determined [1,2].Furthermore, the possibility that the DM is a thermal relic is particularly exciting, namely that the DM interactions brought it to equilibrium with the Big Bang bath, and its interactions determine its abundance [3,4].Being a thermal relic, the abundance today will be insensitive to any initial conditions before the equilibrium state.Thermal relics are also very predictive, since the relic abundance is determined by size of the dominant interaction during freezeout.
The dominant paradigm of DM has been Weakly-Interacting-Massive-Particles (WIMPs), whose departure from equilibrium is determined by DM-DM annihilations into light Standard Model (SM) bath particles. 1This paradigm has led to the belief that a thermal elementary DM candidate, within standard cosmological history, has an upper bound on its mass of ∼ 100 TeV [17].However recent works [14,18] have shown number of exceptions (for heavy dark matter in other scenarios see the review in [19] and references within).In this work we put these examples in a general framework and derive bounds for a family of such mechanism.
In section II, we present a general framework to study thermal freezeout processes, focusing on mechanisms which provide DM masses higher then the WIMP unitarity bound (superheavy).We emphasize the characterizing properties like slow freezeout and solve for the mass -coupling relation showing excellent agreement with numeric results.
In section III, we use the above results and derive an approximate pertubative unitarity bound for a general process.We show that a consequence of a DM candidate χ beyond the maximal mass allowed by the WIMP unitarity bound is that there is another particle close to the mass of χ, which we will call χ ′ .If the mass of χ ′ is above the WIMP unitarity bound, then there must be a nearly degenerate particle to χ ′ as well, resulting in a chain of superheavy particles.The chain can end only either when we reach a particle below the WIMP unitarity bound or with an unstable particle.Such cases are exemplified in scenarios within a standard cosmological history [14,18] and within an early matter dominated scenario [20].Finally, we elaborate on this concept and analyze two examples of chains in section IV.

II. THERMAL FREEZEOUT
In this section we consider a DM candidate in thermal equilibrium with the bath and in chemical equilibrium via the n → m process of the form The general form of the Boltzmann equation for the χ number density is neglecting Pauli-blocking and stimulated emission.Here ⟨σv⟩ is the thermally averaged cross-section of the n → m process in Eq. ( 1).On dimensional grounds, Table I.Parameters and the maximum mass allowed by perturbative unitarity for different freezeout processes.Here χ is the DM, ψ is particle close to mass with χ, and ϕ is a particle much lighter than both χ and ψ.The definition of ∆ is implicitly defined in the Masses column for each process.Both ψ and ϕ are assumed to be in thermal and chemical equilibrium throughout the freezeout process.In the 2 → 2 freezeout cases, we assumed s-wave annihilations for calculating β.For coannihilations, the process ψψ → ϕϕ controls the abundance, where χχ ↔ ψψ is assumed to be in equilibrium throughout freezeout.
this thermally averaged cross section will have the form for some coupling α eff , and parameters b ′ and c ′ .Here we define the dimensionless time parameter, x = m χ /T .The power of α eff is chosen to be the form expected if the amplitude is determined solely by cubic interactions, such as gauge or Yukawa interactions, but this need not be the case.Usually, the temperature dependence can be well approximated at freezeout by a power law or exponential dependence on the temperature, which occurs for instance for forbidden channels.Defining the yield Y = n/s, the Boltzmann equation takes the familiar form (4) Typically, this is the form of the Boltzmann equation that is solved numerically for the asymptotic value of the DM abundance, Y χ (∞).An estimate for the relic abundance of χ can be obtained by the instantaneous freezeout approximation.In this approximation, it is assumed that the DM density departs from equilibrium and instantly stops annihilating, where its co-moving density freezes to the value when it departed equilibrium.Departure from equilibrium, which defines the temperature x d , occurs roughly when 2 2 Note that the definition of x d varies in the literature.In self annihilation scenarios (n = k = 2) it is common to define x d by H(x d ) = c⟨σv⟩n eq χ , where c is chosen such that the instantaneous freezeout approximation provides the best fit to numerical solution of the Botlzmann equations [21,22].In this Letter we focus on k = 1 process, where the definition in Eq. 5 gives The annihilation rate is given as: for some constants, a, b, and c, which are defined by this equation.The values b in different freezeout processes studied in the literature can be seen in Table I.For instance, for WIMP like χ-χ annihilations, using Eqs.( 3) and ( 6), we see that b = 1 and c = 3/2 corresponds to swave annihilation while c = 5/2 corresponds to p-wave annihilation.
Requiring that the co-moving density at the time of freezeout matches the one observed today by Planck [23] yields the constraint where T eq = 0.8 eV is the temperature at matterradiation equality.Assuming instantaneous freezeout, namely that the DM density abruptly changes from an equilibrium density to a fixed co-moving density at for a DM particles freezing out while non-relativistic; g ⋆ (g ⋆s ) is the relativistic (entropy) degrees of freedom and g χ is the number of χ internal degrees of freedom.Using Eqs. ( 5) -( 8), one finds that the mass-coupling better results (i.e.Ref. [14,16]).Such choice of definition is motivated analytically by the Yχ prefactor in Eq. 4, and gives close results both for the k = 1 and k > 1 cases, as shown in the right panel of Fig 1.
relationship that gives the observed abundance is where m pl = 2 × 10 18 GeV is the reduced Planck mass and x d is the solution to Overall numerical factors are left out for simplicity.These factors can be recovered here and later with the replacement where Ω B and Ω DM are the observed baryonic and DM relic fractional densities, respectively.The factor is typically an O(1 − 10) correction to the mass-coupling relationship.
Instantaneous freezeout is a good approximation in many scenarios, because the annihilation rate is dropping exponentially with time (or temperature).However, when the exponential suppression is weaker (b < 1), instantaneous freezeout may no longer be a good approximation for estimating the relic abundance, as shown in Ref. [14,18].For b ≪ 1, corresponding to an annihilation with a much lighter particle, the correction can be large.
Using the parameterization above, b < 1 is possible only when k = 1.Therefore, we consider here corrections to instantaneous freezeout when there is only one DM particle in the initial state, although the method here can be applied to any scenario.The Boltzmann equation ( 4) for k = 1, using the definitions above, can be written as There exists an integral solution to this Boltzmann equation.However, an easy way to get a good estimate for the correction to Eq. ( 9) is to evolve the equation from the moment χ departs equilibrium, which we approximate as occurring at x = x d .For x > x d the inverse reaction can be dropped.The solution is x −c+1 e −bx d (x−1) dx .
(13) We define a parameter β as a measure of the deviation from instantaneous freezeout as For β = 1, the instantaneous freezeout approximation is perfect, while Eq. ( 13) provides an additional possible exponential suppression compared to the instantaneous freezeout case given in Eq. ( 8), where In this form, the mass-coupling relationship in Eq. ( 9), can easily be found as before, but using Eq. ( 14) instead of Eq. ( 8).One finds where Again, overall numerical factors can be restored with the replacement Eq.(11).Notice that Eq. ( 15) and Eq. ( 17) still provide two equations needed to solve for x d and β.However, the integral β can usually be done semi-analytically.For b ≳ 1, where instantaneous freezeout it expected to be good, one finds that β ≃ 1 + 1/(bx d ).For bx d ≪ 1 (and further more for c → 0), β can be quite large.In the extreme b → 0 case, β ≃ 1 + 1/(c − 2).Note that for b = 0 and c ≤ 2, freezeout never occurs, because the annihilation rate is decreasing slower relative to the Hubble rate; this occurs for decay processes and can also occur for co-scattering with a massless mediator.The values of β for various processes are given in Table I. 3In left panel of Fig. 1 we show a schematic of freezeout, exemplifying the meaning of the different parameters defined in this section.In the right panel of Fig. 1 we show the mass vs coupling relationship that matches the observed abundance.The solid line indicates the results from the numerical simulation, while the dashed line indicates the analytical solution in Eq. ( 16).The numerical and analytical results are in excellent agreement across many different freezeout processes.10) and (17).After decoupling, the abundance continues to reduce by a factor e −x d (β−1) until it completely freezes out to a final yield, mχYχ(∞) = 0.55Teq.In gray we also show the abundance curves for for a mχ = 1 GeV DM candidate whose relic abundance is determined by a coscattering or 3 → 2 (SIMP) freezeout process.For coscattering, decoupling happens very early, freezeout is slow and β is quite large.On the other hand, SIMP DM freezes out very fast and instanaeous freezeout works as a better approximation of the final relic abundance.Right: Mass vs. coupling relationship that matches the observed abundance.The solid line shows the results from the numerical simulation, while the dashed line shows the analytical solution in Eq. ( 16).The definitions of ∆ for co-scattering, zombie, and forbidden annihilations are given in Table I.

III. UNITARY BOUND
Partial-wave unitarity can be used to derive upper limits on the size of cross-sections.Applying this to DM-DM annihilations, Ref [17] showed that partialwave unitarity gives an upper limit on the mass of the dark matter from thermal freezeout.The partial-wave unitarity bound on the cross-section for an arbitrary n → m process can be found, and then extended to other freezeout processes.However, partial-wave unitarity bound on cross sections is similar in magnitude to the perturbitivity bound on the couplings.Therefore, we simply use perturbitivity as a proxy for partial-wave unitarity, requiring the maximum coupling be α eff ∼ 4π.
From Eq. ( 9), the perturbative unitarity bound for WIMP-like DM-DM annihilations is immediately apparent by substituting b = 1, m = n = 2, and α eff ∼ O(10): For perturbative couplings, one sees that the mass of a WIMP DM candidate is never greater than a few hundred TeV.However, going beyond this mass is still possible for a thermal relic, provided that b < 1-in other words, that the annihilation rate is less exponentially suppressed as the DM density is annihilating away. 4 4 There are another two interesting possibilities to get an exponentially fast rate, but they require species to be out of chemi- The unitarity bound for an arbitrary interaction is given by Thus for sufficiently small b and c ≤ 1, one can expect that there might be a perturbative thermal candidate for DM all the way up to the Planck scale.As we will now show, going above the perturbative WIMP unitarity bound generically implies that there are additional nearly degenerate particles to the DM, and that the DM is metastable.To go beyond WIMP unitarity, first consider the process setting the abundance.We need only consider the case with one DM in the initial state, otherwise b ≥ 1.The rate for this process scales as cal equilibrium.One is to consider enhancement via stimulated emission if the DM annihilates to a boson with large occupation number.Another possibility is if DM annihilates with a particle that has already frozen out and therefore has a large chemical potential.See Refs.[11,12,14,26,27] for freezeout mechanisms that involve particles chemically decoupled from the bath.

Unstable DM Forbidden Channels
Here i mi is the sum of the initial particle masses excluding one χ particle, and f m f is the sum of the final particle masses.The green region shows the parameter space given by Eqs. ( 22) and ( 23) where b < 1 and dark matter can be above the WIMP unitarity bound.The gray region shows where Eqs. ( 24) is not satisfied and the χ is expected to decay too fast to be the observed DM.The blue and yellow regions show where b > 1, corresponding to 'Forbidden Channels' where the total mass in the final state is heavier than the initial state, and 'Rare Channels' where the DM is annihilating off of particles (cumulatively) rarer than itself.Also shown are contours (red-dashed) of constant b for b = 0, 0.5, 1, 2, and the corresponding upper bound on the DM mass.The upper bound is found is using α = 10 and c = 3/2.
where the second case is for forbidden DM.Requiring b < 1 imposes Finally, we require that the DM is stable from the decay induced by moving all initial particles to the final state These constraints are shown in Fig. 2. The green area shows the region where DM can be above the WIMP unitarity bound.A simple constraint that can be found by combining these conditions gives one of the main results of this Letter, Although Eq. ( 25) is not the strongest constraint that can be derived (see Fig. 2), a direct result of Eq. ( 25) is that if the DM mass is beyond the WIMP unitarity bound, then there is another particle close to its mass, which we will call χ ′5 .Since χ ′ mass exceeds the WIMP unitarity bound as well (m χ ′ ∼ m χ ), its abundance will be too large, unless it has some process that efficiently depletes its abundance.If it is unstable, then it may potentially cause the DM itself to be unstable.If it is stable, it must annihilate away via a non-WIMP like process, and Eq. ( 25) applies to χ ′ as well, resulting in a second particle, χ ′′ , with mass close to χ ′ (m χ ′′ ∼ m χ ′ ).Repeating this argument iteratively for χ ′′ and beyond, leads to a chain of DM particles and interactions.The chain will end either when the mass of last particles goes below the WIMP unitarity bound, or when it reaches a particle that decays in equilibrium with the bath.In the former case the DM may be absolutely stable, but in the the latter case, the DM will be metastable.
As we have shown, a DM chain is a generic consequence of stable or metastable thermal superheavy DM within standard cosmological history.In the next section, we present two examples of such DM chains: a zombie-type chain that gives raise to stable DM, orders of magnitude heavier than the WIMP unitarity bound, with only small number of particles in the chain.We include a second example of an inverse decay chain, which result in metastable DM with very high masses.
Eq. ( 25) gives a constraint on the masses of the particles involved in the freezeout process, when we require that the DM mass is above the WIMP unitarity bound (b = 1).A stronger constraint can be obtained if we require that the DM mass be even heavier.For instance, we can impose that the DM is heavier than the mass given by the unitarity limit in Eq. ( 19) for arbitrary b.In Fig. 2, the dashed orange curves show contours of constant b, and their corresponding unitarity bound on the dark matter mass.

IV. FREEZEOUT CHAINS
Having established a roadmap to going beyond the WIMP unitarity limit for thermal DM, we now turn to new mechanisms that predict super heavy DM.We will focus on two chain mechanisms allowing for orders of magnitude larger mass than the WIMP unitarity bound.First we discuss a zombie chain, in which absolutely stable DM is possible well beyond the WIMP unitarity bound with only a small chain.Next we study an inverse decay chain, where DM is metastable, but whose mass may be much larger.

A. The Zombie Chain
Next we study a chain based on zombie typeinteractions.It was shown in [14] that a zombie process can allow for a heavy DM candidate without a chain, but with an unstable DM particle.Here we introduce a chain of zombie interactions, and show that it supports exponentially larger mass for each additional particle in the chain, up to the bound set by Eq. ( 16), while allowing for completely stable DM.
The zombie chain consists of N DM particles, χ i (i = 1...N ) with zombie-type nearest neighbor interactions For simplicity, we will assume that each chain interaction has the same strength and that the mass hierarchy is constant going from the heaviest to lightest.Namely, we take ⟨σv⟩ i,i+1→i+1,i+1 ≡ α 2 /m 2 χi and m i+1 /m i independent of i.The last particle in the chain, χ N , is assumed to also have direct annihilations into the SM bath particles Since the abundance of χ N is determined by standard self-annihilations, and is assumed to be stable, its mass must be less that the WIMP unitarity bound, i.e., m χ N ≲ 100 TeV.
For the zombie interactions, the annihilation rate is and here b zombie = m i+1 /m i .Following Eq. ( 25) and Eq.(??), we see that in order to go beyond the WIMP unitarity bound, each link in the chain must satisfy The highest mass is obtained for the smallest allowed value b = 1/3; below this value, the DM is unstable.Furthermore, for the interactions above, each χ i is absolutely stable, although each will have an exponentially smaller relic abundance than the previous.Precise solutions to the relic abundance can be obtained by solving the coupled Boltzmann equations.However, the freezeout occurs sequentially: first χ 1 freezes out from χ 2 , then χ 2 freezes out from χ 3 , and so on.At each stage each particle freezes out with its neighbors still in chemical equilibrium.For each χ i , we can simply use Eqs.( 9) and (10) to determine the abundance.
In the left panel of Fig. 3, we show the numerical solution for N = 18, b = 0.337, α eff = 1, and m χ1 = 3.6 × 10 9 GeV, which exceeds the WIMP unitarity bound Eq. (18).The thermal evolution of the comoving density is shown for each χ i .The abundances of χ 1 through χ N −1 freezeout via zombie processes and have decreasing relic abundance m χi Y χi /m χi+1 Y χi+1 ≃ 70.The abundance of χ N freezes out via self annihilations Eq. ( 27), and therefore has a relatively large relic abundance, but much smaller than the χ 1 abundance.

B. The Inverse Decay Chain
Consider freezeout of χ via inverse decay (INDY) process [12,16] where γ ′ is a light particle in equilibrium with the bath.It is easy to see from detailed balance that the rate for depleting χ is where Γ ψ→χ+γ ′ is the partial width of ψ.For inverse decays, in terms of the definition of the rate in Eq. ( 6), and The suppression of the rate in temperature is only due to the mass-splitting between ψ and χ.Therefore, the rate only decreases when the mass splitting is of order the temperature, i.e. x∆ > 1.For sufficiently small ∆, decoupling and freezeout can occur arbitrarily late.For these reasons it is expected that arbitrarily high DM masses are possible.However, for the reasons discussed in the previous section, if m χ is well above the WIMP unitarity bound, then there needs to be a way for ψ to deplete in order to not be over abundant.If it is unstable while decaying in equilibrium, then certainly χ will be too short lived to be DM.Thus, we need to consider ψ depleting by inverse decays as well.To this end, we consider a chain of inverse decays in order to achieve arbitrarily high DM mass.The chain must be long enough so that the phase space in the χ decay is small enough to stabilize the DM on cosmological scales.Consider the chain of inverse decay processes for i = 1 . . .N , and the decay for χ N directly into the Standard Model bath The decay rate χ N is taken to be fast, so that it is an equilibrium process, i.e., Γ χ N →sm+sm > H(m χ N ).
For simplicity we consider the case where the decay width for each process in the chain is constant Γ ≡ Γ χi+1→χi+γ ′ and a constant mass splitting m χi+1 = m χi (1 + ∆).
The relic abundance can be found by solving the N coupled Boltzmann equations: . . .
where x = m χ1 /T .We show a numerical solution in the right panel of Fig. 3 for N = 30, ∆ = 0.01, α eff = 1, and m χ1 = 7.35 × 10 11 GeV, which exceeds the WIMP unitarity bound, Eq. ( 18).One can see that the analytic mass-coupling relationship derived in Eq. ( 16) does not well-describe the results here.Similar to what was observed for the coscattering chain in Ref. [18], the chain can weaken the strength of the DM depletion.For scattering, it was shown that the effective rate for annihilation is suppressed by N 2 .In other words, one should replace α eff → α eff /N 2 in order to obtain an estimate for correct coupling needed to match the observed abundance.Exact analytic solutions to the coupled set of equations are difficult to obtain, and we leave this to future work [28].
Finally, we comment on the stability of the DM for the example shown.The DM candidate, χ 1 , can decay to N SM particles via N − 1 off-shell χ i particles.We estimate the phase-space and choose N = 30 as a safe value that guarantees that the phase space in the χ decay is small enough to stabilize the DM on cosmological scales.

V. OUTLOOK
In this Letter we have derived general formula for the thermal relic abundance of dark matter for an arbitrary freezeout process, and have outlined the roadmap to obtain a thermal relic with mass above the standard WIMP perturbative unitarity bound of a few hundred TeV.
Thermal freezeout provides a well-motivated framework for studying the relic abundance of DM.The relic abundance can be determined in the instantaneous freezeout approximation for most cases with Eq. ( 9), or with Eq. ( 16) which corrects for when freezeout is slow.Depending on the process controlling freezeout, there is a different bound on the dark matter mass from perturbative unitarity, which can be substantially higher than that of the WIMP standard lore.To go beyond the standard WIMP unitarity bound, we find that there are particles degenerate with the DM, possibly constituting a long chain of DM interactions.Additionally, if the DM is a heavy thermal relic within a standard cosmological history, then generically it is expected be meta-stable.This would lead to a striking signal of ultra-high-energy cosmic rays (UHECR) from decaying DM, that can be searched for in dedicated UHECR detectors such as Ice Cube [29] and Pierre Auger Observatory [30] (in addition to gamma-ray satellites such as FERMI-LAT [31]).
Much work remains to be done to better understand heavy dark matter in context of a chain and to populate the model parameter space that achieves such heavy thermal relics.We leave this to future work [28].

Figure 1 .
Figure 1.Left: mχYχ vs.x for DM freezeout.The orange curves shows the solution of the Boltzmann equation for a mχ = 1 GeV WIMP.The red point shows the value x d , given by Eqs(10) and(17).After decoupling, the abundance continues to reduce by a factor e −x d (β−1) until it completely freezes out to a final yield, mχYχ(∞) = 0.55Teq.In gray we also show the abundance curves for for a mχ = 1 GeV DM candidate whose relic abundance is determined by a coscattering or 3 → 2 (SIMP) freezeout process.For coscattering, decoupling happens very early, freezeout is slow and β is quite large.On the other hand, SIMP DM freezes out very fast and instanaeous freezeout works as a better approximation of the final relic abundance.Right: Mass vs. coupling relationship that matches the observed abundance.The solid line shows the results from the numerical simulation, while the dashed line shows the analytical solution in Eq. (16).The definitions of ∆ for co-scattering, zombie, and forbidden annihilations are given in TableI.

Figure 3 .
Figure 3. Thermal evolution of the comoving energy density mχ i Yχ i vs. temperature x = mχ i /T for the zombie chain (left) and inverse decay chain (right).The solid lines shows the numerical abundance for each particle, while the dashed lines represent equilibrium abundance curves.In each case, we indicate the DM mass, mass-splitting, effective coupling and number of particles in the chain on the plot.