Unitarity limits on thermal dark matter in (non-)standard cosmologies

Using the upper bound on the inelastic reaction cross-section implied by S-matrix unitarity, we derive the thermally averaged maximum dark matter (DM) annihilation rate for general $k \rightarrow 2$ number-changing reactions, with $k \geq 2$, taking place either entirely within the dark sector, or involving standard model fields. This translates to a maximum mass of the particle saturating the observed DM abundance, which, for dominantly $s$-wave annihilations, is obtained to be around $130$ TeV, $1$ GeV, $7$ MeV and $110$ keV, for $k=2,3,4$ and $5$, respectively, in a radiation dominated Universe, for a real or complex scalar DM stabilized by a minimal symmetry. For modified thermal histories in the pre-big bang nucleosynthesis era, with an intermediate period of matter domination, values of reheating temperature higher than $\mathcal{O}(200)$ GeV for $k \geq 4$, $\mathcal{O}(1)$ TeV for $k=3$ and $\mathcal{O}(50)$ TeV for $k=2$ are strongly disfavoured by the combined requirements of unitarity and DM relic abundance.


Introduction
Model independent bounds on the mass of dark matter (DM) candidates are extremely weak [1]. Since DM particles must be confined to galaxies, for bosons their de Broglie wavelength should be smaller than the typical size of the DM-rich dwarf galaxies, which is about 1 kilo parsec. This leads to a lower bound of around 10 −22 eV on the mass of bosonic DM [2,3]. For fermionic DM the lower bound is much stronger, about 1 keV, as the Pauli principle sets a maximum value of the phase-space density [4]. The stability of stellar clusters in galaxies also requires a very large upper bound of 10 3 solar masses on the mass of any DM candidate [5,6]. The allowed range of DM mass may be further restricted only if specific properties of DM particles are assumed. One such well-known example is the assumption that DM was in kinetic equilibrium in the early Universe with the standard model (SM) thermal bath. The free-streaming of such thermal dark matter is constrained by the Lyman-α flux-power spectra data, resulting in a lower limit of around 5.3 keV [7].
Specifying a particular production mechanism of the DM abundance in the early Universe may indicate a viable DM mass range for that scenario, though the allowed values tend to be model dependent. For DM particles that were in kinetic and chemical equilibrium in the early Universe with the SM bath, and whose present mass density was determined by number-changing pair annihilations to the SM sector, a model independent upper bound on the DM mass can be obtained from the requirements of unitarity of the Smatrix [8,9]. As first shown by Griest and Kamionkowski [8], such a bound follows from the maximum inelastic reaction rate implied by unitarity, which sets the minimum frozen out number density of DM, and hence the maximum mass that saturates the present density. The presence of long-range interactions which may lead to bound-state formation reduces the effective number-changing inelastic annihilation rate, and thus reduces the unitarity upper limit compared to the scenario with no bound-state effects [10][11][12]. Furthermore, the presence of a particle anti-particle asymmetry in the DM sector necessarily implies a non-zero equilibrium chemical potential for such DM, thereby increasing the effective number density of the surviving species at freeze-out, and reducing the unitarity upper limits further [11,13]. In particular, as shown in Ref. [13], in a purely asymmetric DM scenario, generated by a semi-annihilation process with large CP-violation, the unitarity limit on the DM mass can be as strong as 15 GeV.
The above studies on unitarity limits have focussed on scenarios in which the dominant DM number changing reaction is of 2 → 2 type, in which a pair of DM particles annihilate either to a pair of SM particles (pair-annihilation), or to a DM and a SM particle (semi-annihilation). However, as pointed out by Carlson, Machacek and Hall [14], and subsequently revived in recent studies [15], the dominant number changing interactions may take place entirely within the dark sector as well. The minimal such reaction involving a single DM species is of 3 → 2 type, dubbed as strongly interacting massive particles (SIMP), due to their appearance in theories of the dark sector involving new strong interactions [16]. Generically, a simple low energy effective theory of a complex scalar particle with a cubic and quartic self-interaction will have 3 → 2 number changing reactions.
A natural question therefore is what are the implications of S-matrix unitarity in such scenarios, where the dominant DM number changing interaction is of the type k → n, with k ≥ n? A completely general formulation of this problem is challenging since the partialwave decomposition of a k−body initial state, for k ≥ 3, is rather involved. One might however restrict to a smaller subset of reactions of the k → 2 type, with k ≥ 2, since as we shall see, in a thermal bath the relevant thermally averaged reaction rates σv k−1 rel k→2 get related to σv rel 2→k in equilibrium, and one can of course easily perform a partial wave analysis for the two-body initial state. We therefore only need to find out the maximum value of the inelastic 2 → k cross-sections as implied by unitarity, which can be obtained using the optical theorem and the study of the corresponding 2 → 2 elastic scattering process [17].
A general expression for the thermally averaged maximum rate of k → 2 reactions σv k−1 rel k→2 is the primary result obtained in this paper, as detailed in Secs. 2 and 3. The modification to the upper bound for two identical initial state particles in the 2 → k reaction is also given. Having obtained the upper limit on the annihilation rate, we translate this to a bound on the maximum mass of the DM particle in a radiation dominated Universe in Sec. 4, using the Boltzmann equation for k → 2 reactions set up in Sec. 3. Sec. 4 provides the generalization of the bound of Ref. [8] for k → 2 reactions taking place within the dark sector. We then go on to consider the possibility of a modification to the thermal history of the Universe, taking up the frequently occuring scenario of an intermediate period of matter domination in the pre-big bang nucleosynthesis (BBN) era [18][19][20][21][22][23][24][25][26]. It is observed in Sec. 5 that though the unitarity limits are weaker in such scenarios primarily due to the dilution of DM density from late-time entropy production, there are strong implications of unitarity to the possible values of the reheating temperature at which the radiation dominated Universe is restored. We summarize our results in Sec. 6.

Implications of S-matrix unitarity
To begin with, we recall some of the basic results on the implications of S-matrix unitarity, following the treatment of Weinberg in Ref. [17]. We focus the discussion to the context of dark matter annihilations, and adopt a multi-particle momentum eigenstate normalization convention different from Ref. [17] 1 . A similar approach is also followed by Hui [9]. For a two-particle initial state α, unitarity of the S-matrix implies the optical theorem, which reads in the centre of momentum (CM) frame as [17] Im M αα = 2| p|E CM σ tot , where, | p| is the magnitude of the three momentum of each initial particle and E CM is the total initial state energy, both in the CM frame. The state label α is described by the three momenta, spin-z components (or helicity) of each particle, and all other internal quantum numbers that label the particles, and σ tot = β σ α→β , where β denotes all possible final states that can be obtained from the initial state α. Let us consider for simplicity collision of two non-identical spin−0 particles; the generalization to non-zero spin is straightforward [9,17] 2 . We shall discuss the case for identical particles subsequently. To utilize the rotational invariance of the problem, we perform the basis transformation from the momentum eigenstates | p 1 , p 2 , n to the states | P , E, , m, n , where P = p 1 + p 2 , E = E 1 + E 2 , and ( , m) are a pair of integers such that 2) where Y m l (p 1 ) are the spherical harmonics corresponding to the direction unit vectorp 1 . Here corresponds to the total orbital angular momentum of the two initial state particles in the CM frame, and m to its z−component. The channel index n stands for the two particle species' labels n 1 and n 2 , and similarly for n . The normalization of the scalar product is chosen to ensure that in the CM frame, the state vector |0, E, , m, n has the following inner product with a general state vector | P , E , , m , n P , E , , m , n |0, E, , m, n = δ (3) P δ E − E δ , δ m,m δ n ,n . (2.3) With this, we can write the matrix elements for the operator S − I as
where, we have subtracted out the trivial no-scattering part of the S-matrix, and consider only the connected part of the S-matrix operator, S −I. The matrix element is independent of m and is only a function of and E, as the operator S commutes with the generators of rotation. With S βα = δ (β − α) + (2π) 4 δ (4) (p β − p α ) (iM βα ), we then obtain for a 2 → 2 scattering with | p 1 , p 2 , n → | p 1 , p 2 , n , in the CM frame, Choosingp 1 along the z−direction, and integrating over dΩ(p 1 ), one obtains for the elastic 2 → 2 cross-section in the CM frame in the channel n (2.6) Now using the optical theorem in Eq. 2.1 and the general expression for the matrix element in Eq. 2.5, we find the total cross-section from the channel n, which includes both elastic and inelastic processes, as σ total = π | p 1 | 2 (2 + 1) 2 Re (1 − S n,n ( , E)) . (2.7) Finally, subtracting the elastic cross-section in Eq. 2.6 from the total cross-section in Eq. 2.7 gives us the total inelastic scattering cross-section from the channel n (2.8) As |S n,n ( , E)| 2 ≥ 0, we obtain the upper bound on the total inelastic cross-section implied by unitarity σ inelastic ≤ π | p 1 | 2 (2 + 1) , (2.9) as also derived in [9]. For the annihilation of a pair of non-relativistic DM particles of mass m χ such that | p 1 | m χ v rel /2, with v rel being the relative velocity of the colliding particles, Eq. 2.9 implies the well-known result [8] assuming that the 2 → 2 process under consideration saturates the unitarity upper limit. So far, the discussion was focussed on non-identical initial state particles. For identical particles in the initial state, we modify the normalization in Eq. 2.2 with an additional factor of √ 2, (for identical initial particles), (2.11) such that the angular momentum eigenstate normalization remains the same as in Eq. 2.3. With this the expression for elastic scattering cross-section in the channel n becomes σ elastic = 2π | p 1 | 2 (2 + 1) | (S n,n ( , E) − 1) | 2 (for identical initial particles). (2.12) Similarly, the expression for total cross-section from the channel n also gets modified to σ total = 2π | p 1 | 2 (2 + 1) 2 Re (1 − S n,n ( , E)) (for identical initial particles), (2.13) thus modifying the total inelastic cross-section from the channel n and its corresponding upper bound 2π | p 1 | 2 (2 + 1) (for identical initial particles). (2.14) Our results for the maximum value of the inelastic scattering cross-section for identical initial state particles does not agree with the corresponding comment in Ref. [9], but agrees with the appearance of this extra symmetry factor (of 2, for 2 identical particles in the initial state) with Ref. [29]. However, the maximum values of the inelastic k → 2 reaction rates obtained in Refs. [29,30] do not agree with our results in the next section, as the authors in Refs. [29,30] did not maximize the inelastic rate while taking into account the fact that for a non-zero inelastic reaction rate, the elastic scattering cross-section is always non-zero as well [9,17,31]. As seen in this section, our derivation of the maximum rate of inelastic reaction cross-sections uses only the optical theorem and the matrix element and cross-section for the relevant 2 → 2 elastic scattering process.

General Boltzmann equations for k → 2 DM annihilation processes
We now consider the implications of the upper bound from unitarity on the total inelastic reaction rate, on the number-changing dark matter annihilations in the early Universe. The evolution of the number density of dark matter particles in the expanding universe is described by a set of coupled Boltzmann equations. With general k → 2 collision reactions, for k ≥ 2, the number density of the i−th particle n i satisfies the equation (3.1) where H is the Hubble expansion rate of the Universe. The product over α includes the momentum integral factors dΠ α for all the (n + 2) particles, and the sum over the reaction channels indicates all possible reactions involving the production and destruction of the i−th particle, with the net change in the i−th particle number being ∆n i in a reaction.
Appropriate symmetry factors should be included in the momentum integrals to take into account the presence of identical particles. Here, we define and the number density of each species n α (t) is given by an integral over the distribution The function ω ( p α ), ignoring Pauli blocking and Bose enhancement factors, is given as (3.4) where, the matrix elements have been summed over both initial and final spins. We can rewrite Eq. 3.1 in terms of the unpolarized cross-sections of the above reactions as where, the thermally averaged reaction rate is given by where, σ k→2 is the unpolarized (summed over final spins and averaged over initial spins) cross-section for the k → 2 process, v rel is the relative velocity of each particle pair, and f eq α is the equilibrium distribution function for the particle species α. The appropriate symmetry factors appearing in Eq. 3.1 should now be included in Eq. 3.6. The thermal average σ 2→k v rel can be similarly performed. Eq. 3.5 can be further simplified by noting that in equilibrium the following relation is satisfied for each individual reaction channel: Thus, we can express the Boltzmann equations only in terms of σ 2→k v rel by using Eq. 3.7 in Eq. 3.5 ..n k n eq 1 n eq 2 ...n eq k − n a n b n eq a n eq b (3.8) We can easily perform the thermal average integral in Eq. 3.6 for the 2 → k reactions with the maximum value of the inelastic cross-section as in Eq. 2.9 as input, and obtain, with all k + 2 particles having the same mass, where, x = m χ /T . For the DM scenarios to be considered subsequently, the k + 2 particles will be of the same species, including anti-particles, and will thus have the same mass, m χ . Here we have assumed that the unitarity limit on the total inelastic cross-section is saturated by the specific 2 → k reaction under consideration. With that, for k = 2 and = 0 this reduces to the well-known result [8] For k ≥ 3 we see that there is an exponential suppression factor in σ 2→k v rel max , namely e −(k−2)x , due to the phase-space cost for producing each extra particle. We can now use Eq. 3.7 to obtain the maximum value of the thermally averaged rate σ k→2 v k−1 rel as follows: Here g χ is the number of spin degrees of freedom of the DM particle. Thus, for example, the maximum value of the thermally averaged s-wave cross-section for a 3 → 2 reaction is given by Similarly, the maximum value of the thermally averaged s-wave cross-section for a 4 → 2 reaction is given by As mentioned at the end of the previous section, our results for the thermally averaged reaction rates in Eqs. 3.11, 3.12 and 3.13 do not agree with the results in Refs. [29,30]. For identical particles in the initial state, we found in Eq. 2.14 that the maximum inelastic cross-section for a 2 → k reaction is a factor of two larger than the non-identical case. However, the thermal averaging integral in σ 2→k v rel will in this case have a symmetry factor of 1/2 for the two identical particles in the initial state. Therefore, for identical initial state particles, σ 2→k v rel max as shown in Eq. 3.9 remains valid, and similarly Eq. 3.11 remains the same as well.

Unitarity limits on thermal DM mass: radiation dominated Universe
We can now apply the results on the unitarity upper bound for the thermally averaged annihilation rates to find out the limits on the DM mass. To begin with, we consider the standard scenario of a radiation dominated Universe in the early epochs before big-bang nucleosynthesis. In such a Universe, the Boltzmann equation for a dark matter particle thermalized with the SM bath, going through k → 2 annihilations within the same species (i.e., when all the k + 2 particles involved in the collision are the same), can be written using Eq. 3.8 as where, we have defined Y = n χ /s and x = m χ /T , with s being the entropy density, with the assumption that it is conserved during the evolution considered above, i.e., s(t)a(t) 3 = constant, where a(t) is the scale factor describing the expansion of the Universe in the Friedmann-Robertson-Walker cosmology. Taking the maximum value of the reaction rate using Eq. 3.9, we then arrive at the equation where, λ R is given by Here, σ 0 = 4 √ π/m 2 χ , g * ,s (T ) and g * (T ) are the effective number of relativistic degrees of freedom relevant for the entropy density and Hubble rate, respectively, and s(m χ ) and H(m χ ) are the entropy density and Hubble rate evaluated at the temperature T = m χ .
Eq. 4.2 can be solved analytically with the approximation that for large values of x, which correspond to the late time Universe, we may ignore the Y k−2 eq term compared to the Y k−2 term, and that the relativistic degrees of freedom do not change appreciably during the DM freeze-out. Then, with the boundary condition Y ( where, σ 2→k v rel is to be evaluated at x = x F . With this boundary condition, the approximate solution for the present DM yield is obtained to be (4.5) We can thus obtain the approximate value of the DM relic abundance as Ω χ = m χ Y (x → ∞)s 0 /ρ c , where s 0 is the present entropy density and ρ c is the critical density.
In order to obtain accurate values of the upper bound on the DM mass, we solve the Boltzmann equation 4.1 numerically in several simple scenarios in which the dominant DM number changing topology is of k → 2 type, with k ≥ 2. In each case a minimal scenario with (real or complex) scalar is considered, with a minimal Z N DM stabilization symmetry. The scenarios are as follows: 1. χ + χ → SM + SM: Here, χ is a non-self-conjugate DM candidate with a distinct anti-particle χ, stabilized by a Z 2 symmetry. The Boltzmann equation in this case is the familiar one 2. 3χ → 2χ: Here, χ is a complex scalar DM stabilized by a Z 3 symmetry. Assuming CP-conservation, the Boltzmann equation for the DM number density, taking into account all reactions allowed by the Z 3 symmetry, can be written as 3. 4χ → 2χ: Here, χ is a real scalar DM stabilized by a Z 2 symmetry. The Boltzmann equation in this case, taking into account all reactions allowed by the Z 2 symmetry, is 4. 5χ → 2χ : Here, χ is a complex scalar DM stabilized by a Z 5 symmetry. Assuming CP-conservation, the Boltzmann equation here, taking into account all reactions allowed by the Z 5 symmetry, is given as: We show the resulting unitarity upper limits on the DM mass in the above scenarios in a radiation dominated Universe, in Table 1, by requiring that the DM particle saturates the observed DM abundance of Ω DM h 2 0.12 [32]. The upper bounds are shown for three cases in each scenario, with the dominant annihilation mode being from (1) s-wave, (2) p-wave and (3) both s-and p-wave initial states. As the rate for each higher partial wave is expected to be further suppressed by powers of a non-relativistic relative velocity [8], we restrict our considerations to the s-and p-wave contributions only. As we can see from this Table, the well-known upper bound of around 130 TeV for non-identical s-wave DM pair annihilations to SM states is reproduced. The upper bound for the 3 → 2 scenario is around 1 GeV, while for the 4 → 2 and 5 → 2 scenarios it is around 7 MeV and 110 keV, respectively. For all the k → 2 scenarios with k ≥ 3, the bounds do not change significantly on inclusion of higher partial wave contributions to the annihilation rate, since the maximum value of the mass scales as (2 + 1) 1/k . For higher values of k, the annihilation rates are flux factor suppressed, thereby reducing the unitarity upper limit on the annihilation rate, thus increasing the resulting DM number density, and lowering the upper bound on the DM mass that saturates the observed DM abundance. For non-zero DM spin, angular momentum conservation will impose further selection rules, allowing only certain annihilation topologies, and within each topology possibly requiring specific values of .   (3) both s-and p-wave initial states.

Unitarity limits on thermal DM: intermediate matter dominated Universe
of matter-radiation equality at a temperature of around 1 eV. The successful predictions of big-bang nucleosynthesis (BBN) require the Universe to be RD at temperatures of the order of 1 MeV [1]. However, in the pre-BBN era, the energy density could be dominated by a non-relativistic matter field, which eventually decays to radiation sufficiently before the BBN, restoring back an RD Universe. This requires the heavy matter field (Φ) to be very feebly interacting with the SM bath, such that it is long-lived and does not thermalize with the SM sector, during the epoch of interest. For this section, we further assume that the Φ field does not possess sufficiently strong number-changing interactions within its own sector either, and explore the consequences of relaxing this assumption elsewhere. Such fields may, for example, interact with the SM sector only gravitationally, as often encountered in extensions of the SM of particle physics.
If a sufficient density of Φ is generated at the end of the inflationary reheating along with radiation, then eventually due to the faster dilution of the radiation energy density from the cosmic expansion, compared to that of non-relativistic matter (namely, ρ R ∝ 1/a 4 and ρ Φ ∝ 1/a 3 in the absence of any net creation or destruction), the Φ field energy becomes equal to that of radiation at a time t = t MR , say. This results in the onset of a stable matter dominated era which continues during t MR t t MD , when at around t t MD the decays of the Φ particles start to produce significant amount of radiation, and thus entropy. This leads to a Universe with energy density dominated by a decaying matter field, during t MD t t RH , when at around t t RH most of the Φ density decays back to radiation, thus restoring the RD phase before the BBN. We schematically show this timeline of cosmological events in Fig. 1, in which this alternative thermal history with an intermediate matter dominated (IMD) phase is contrasted with the standard thermal history, assuming an inflationary scenario at the earliest epoch of the Universe.
The scenario described above is governed by a set of three coupled differential equa- tions, namely, the Boltzmann equations for the energy densities of radiation and Φ, where, Γ Φ is the total decay width of Φ, which is assumed to decay entirely to radiation. This system of equations have been studied extensively in the literature within different contexts, including that of inflationary reheating, see for example, Refs. [1,18,20,21,33]. We shall adopt a simple analytical approximation, which captures the essence of the modified cosmic history involving the interplay of radiation and Φ densities, and then treat the evolution of the dark matter density accurately with this cosmology as an input. Our approach largely follows Ref. [1], with certain improvements in the analytical treatment as will be clear from the discussion below. Eq. 5.2 is easily solved to obtain the time-evolution of ρ Φ (t) as where t I is some initial time after inflationary reheating, when the Φ energy density is ρ Φ (t I ). It is convenient to choose t I in the range t MR t I t MD , such that during this time the Hubble parameter H(t) is approximately determined by ρ Φ , with Φ being an essentially stable non-relativistic matter field. Thus, using Eq. 5.3, with H(t) 2/(3t), we find ρ Φ (t I ) 3M 2 P H 2 (t I ) 4M 2 P /(3t 2 I ), where we have defined the reduced Planck mass by 8πG = 1/M 2 P . With a(t) ∝ t 2/3 during this epoch, we can write using Eq. 5.4 With Eq. 5.5 as an input to the Boltzmann equation 5.1 for the radiation energy density, with the boundary condition ρ R (t MR ) = ρ Φ (t MR ), we obtain the approximate solution for ρ R (t) where, the first term in the RHS stems from the decay of the Φ field, while the second term is the contribution of the radiation density present before Φ decay. For t I ∼ O(t MR ), Γ Φ t I << 1, and hence e Γ Φ t I ∼ 1.
The decay of the Φ field to radiation generates entropy, which we shall now compute. We define t MD to be the time when i.e., at t = t MD , the contribution to the radiation energy density from Φ decay equals the pre-existing radiation energy density in Eq. 5.6. The time scale t MD is approximately given in terms of Γ Φ and t MR as Most of the entropy is essentially generated by the time t = t RH , which is defined by Γ Φ t RH = 1, which is also the time scale for which Γ Φ ∼ H(t). At this time, we find Since the entropy density is dominated by relativistic species, we can express the ratio of the total entropy S(t) in a co-moving volume a(t) 3 between an initial time t i and a final time t f as (5.10) Using Eqs. 5.7 and 5.9 in Eq. 5.10, with t i = t MD and t f = t RH ∼ t MD + δ, with δ → 0, we obtain the ratio of total entropies in a co-moving volume as In the instantaneous decay approximation used above, the relativistic degrees of freedom and the scale factors do not change significantly and therefore cancel out in the ratio in Eq. 5.10.
We can express the product Γ Φ t MR in Eq. 5.11 in terms of the temperatures for matterradiation equality, T MR , and for reheating back to RD, T RH . For this we use the relations (5.12)

DM freeze-out in constant entropy phase
If the entropy production happens mostly after the DM freeze-out, i.e., if the DM particle freezes out during the stable matter dominated epoch, or the preceding radiation dominated epoch, then we can express the ratio of the present dark matter number density (n(T 0 )) and entropy density (s(T 0 )) in terms of the corresponding ratio at freeze-out, and the entropy dilution factor as follows Thus for these two scenarios of DM freeze-out, we can solve the Boltzmann equation determining the DM number density in an iso-entropic Universe Eq. 4.1, and then obtain the present density using Eq. 5.13, with the input Hubble parameter for this epoch (5.14) Since T MR > T in this case, the Hubble expansion is faster than in the corresponding radiation dominated Universe at the same temperature, resulting in an earlier DM freezeout. The Boltzmann Eq. 4.2 is then modified for the adiabatic matter dominated Universe as This equation can be approximately solved for large values of x, which then including the entropy dilution factor gives

DM freeze-out in varying entropy phase
If the DM freezes out during the phase in which the decay of Φ is generating significant entropy, within the interval t MD t t RH , the temperature dependence of the Hubble expansion rate is modified to H Decay (T ) = π 2 g * (T RH ) 90 where, we have matched the Hubble rate to the one in a radiation-dominated Universe at the boundary T = T RH . Further matching the Hubble rate H Φ (T ) and H Decay (T ) at T = T MD , we find a relation between T MD , T MR and T RH . Assuming g * (T MD ) ∼ g * (T RH ), this relation is given by Furthermore, in this scenario, the DM freeze-out is no longer taking place in an iso-entropic Universe. In this case, we can estimate the present DM abundance in terms of the freezeout abundance as follows. For T < T FO , conservation of DM number in a co-moving volume implies that For a temperature range in which the variation of the relativistic degrees of freedom in the thermal bath is not significant, we can show that in the decaying matter dominated phase the scale factor a(T ) ∝ T −8/3 [20,21]. With this, and using s(T ) ∝ T 3 , Eq. 5.20 implies The Boltzmann equation 3.8 can now be solved with the modified Hubble rate in Eq. 5.18 to obtain the number density at freeze-out, which can then be used in Eq. 5.21 for computing the present DM abundance.

Unitarity constraints in the IMD scenario
With the above results, we now explore the implications of unitarity on the scenarios of DM freeze-out in which the IMD epoch plays a significant role. There are essentially three free parameters determining DM properties in a thermal history with an IMD phase, m χ , T MR and T RH , where the DM annihilation rate is given entirely in terms of m χ by unitarity. For each value of T MR and T RH , unitarity thus implies a value of m χ that satisfies the observed relic abundance, which is also the highest value of m χ allowed in that scenario. The highest value of m χ also depends upon whether the DM freeze-out takes place in a nearly constant entropy phase, or in the phase with significant entropy generation. For a given set of T MR and T RH values, this is determined by m χ , since the temperature at the approximate boundary between these two phases, T MD , is obtained through Eq. 5. 19. In order to understand the range of DM mass and annihilation rates allowed by both the relic abundance requirement and unitarity, we parametrize the 2 → 2 reaction rate as where, σ 2→2 v rel 0 is a free parameter, which takes a maximum value allowed by s−wave unitarity as given in Eq. 3.10. Similarly, we parametrize the 3 → 2 reaction rate as where, σ 3→2 v 2 rel 0 is a free parameter, which takes a maximum value allowed by s−wave unitarity as given in Eq. 3.12. Finally, the 4 → 2 reaction rate is parametrized by where again σ 4→2 v 3 rel 0 is a free parameter, whose maximum value allowed by s−wave unitarity is as given in Eq. 3.13. T RH =100 GeV We show in Fig. 2 the region in the σ 2→2 v rel 0 and m χ plane in which the observed value of the DM density Ω DM h 2 = 0.12 is reproduced. The scenario with T RH = 10 MeV is shown in the left panel, while the one with T RH = 100 GeV is shown in the right panel. In each figure we show the results for different values of the matter-radiation equality temperature T MR = 10 3 , 10 9 and 10 15 GeV using the red, blue and green lines, respectively, for the T RH = 10 MeV scenario. To generate a minimum amount of entropy dilution that can impact the DM properties, the corresponding T MR values chosen for T RH = 100 GeV starts from T MR = 10 5 GeV (red line in the right panel), the other two T MR values shown being the same. For comparison, the black solid line (left panel) shows the results in a purely radiation dominated scenario, with no subsequent IMD or entropy dilution. The grey shaded region is disallowed by the unitarity limit on the annihilation rate in Eq. 3.10. As we can see from this figure, for a given T MR and T RH , the required value of σ 2→2 v rel 0 changes with different slopes for increasing values of m χ , depending upon whether the DM freezeout takes place during radiation domination (shown by solid lines), nearly stable matter domination (dashed lines), or decaying matter domination epochs (dot-dashed lines). This is due to the differences in the expansion rate of the Universe in these regions, as shown by the modification to the Hubble rate compared to a RD scenario in Eqs. 5.14 and 5.18. For freeze-out during stable or decaying matter domination, the net DM yield depends on m χ explicitly, in addition to the dependence through x F , and higher values of m χ require lower cross-sections to saturate the observed density. The maximum cross-section allowed by unitarity also decreases with increasing m χ . The intersection point of the relic density allowed line and the unitarity constraint thus gives the maximum allowed mass for a given scenario. In the purely RD scenario (black solid line, left panel), this intersection point gives the unitarity limit quoted in Table 1. For a fixed T RH , increasing T MR leads to a higher dilution of the frozen out DM density. To compensate for that, we need higher DM densities to survive at freeze-out. This implies an earlier decoupling, thereby requiring lower cross-sections for the same mass to satisfy the DM abundance. This is why the unitarity limit on the DM mass that satisfies the relic abundance also increases for higher values of T MR .
Our results in Fig. 2 can be compared with the results in Ref. [34]. The required values of the DM mass and annihilation rates that satisfy the observed DM abundance are in agreement in the two studies. However, in Fig. 2, the unitarity constraints are shown accurately, while Ref. [34] presented only an estimate using O(1) coupling values, and hence the unitarity disallowed regions differ accordingly.  The results for the scenarios in which 3 → 2 and 4 → 2 annihilations are the dominant DM number changing process, are qualitatively similar to the 2 → 2 scenario discussed above, as shown in Figs. 3 and 4, respectively. Since the mass dependence of the maximum annihilation rates in these scenarios are stronger than in the 2 → 2 case, the implied upper bounds on the DM mass are also correspondingly smaller. In order for the IMD epoch to have any impact on the properties of DM, the freeze-out temperature T FO should be larger than T RH . This in turn implies a typical minimum value of the DM mass of around x F T RH , with the value of x F usually lying in the range of 20−30. However, the cross-sections required by the relic abundance condition might already be higher than the unitarity upper limit for this mass. Thus we find that with increasing T RH , the region allowed by the combination of relic density and unitarity conditions gets reduced. This is more pronounced for the 3 → 2 and 4 → 2 annihilation scenarios, in which, as we can see in Figs. 3 and 4, for T RH = 100 GeV, only a small region remains allowed even with T MR as high as 10 15 GeV. How high can T RH be then for this largest value of T MR = 10 15 GeV considered in this study? We show in Fig

Summary
In this paper, we have studied the implications of S−matrix unitarity on general k → 2 thermal DM annihilations, for k ≥ 2. We first derived the upper limit on the inelastic 2 → k annihilation cross-sections using the optical theorem, and the matrix element and cross-section for the elastic 2 → 2 annihilation process. These limits are obtained both for non-identical and identical initial state particles, where for the latter case they are larger by a factor of two, with two particles in the initial state. The thermally averaged k → 2 annihilation rates in the Boltzmann equation are expressed in terms of the 2 → k annihilation rates using detailed balance. On thermal averaging, σ 2→k v rel turns out to be the same for both the identical and non-identical initial state particles. The unitarity limits on the annihilation rates translate to upper limits on the DM mass, which generalize the well-known result for 2 → 2 DM pair annihilation to a pair of SM particles, with an upper bound of around 130 TeV for non-identical s-wave DM pair annihilations. The s-wave upper bound for the 3 → 2 scenario is found to be around 1 GeV, while for the However, on the other hand, increasing T RH leads to a reduction of the region allowed by the combined constraints of the observed DM abundance and unitarity, eventually strongly disfavouring values of reheating temperature higher than O(200) GeV for k ≥ 4, O(1) TeV for k = 3 and O(50) TeV for k = 2, with T MR 10 15 GeV.