Unitarity Bound on Dark Matter in Low-temperature Reheating Scenarios

Model-independent theoretical upper bound on the thermal dark matter (DM) mass can be derived from the maximum inelastic DM cross-section featuring the whole observed DM abundance. We deploy partial-wave unitarity of the scattering matrix to derive the maximal thermally-averaged cross section for general number-changing processes $r\to 2$ (with $r\ge 2$), which may involve standard model particles or occur solely within the dark sector. The usual upper limit on the DM mass for $s$-wave annihilation is around 130 TeV (1 GeV) for $r=2$ (3), only applies in the case of a freeze-out occurring in the standard cosmological scenario. We consider the effects of two nonstandard cosmological evolutions, characterized by low-temperature reheating: $i)$ a kination-like scenario and $ii)$ an early matter-dominated scenario. In the first case, early freeze-out strengthens the unitarity bound to a few TeVs for WIMPs; while in the second case, the WIMP DM can be as heavy as $\sim 10^{10}$ GeV due to a large entropy dilution.


Introduction
The omnipresence of nonbaryonic dark matter (DM) has been confirmed from celestial observations at different scales, from individual galaxies to clusters of galaxies, and even at the cosmological scale [1].Several properties of DM also emerge from such observations, albeit only by exploiting its gravitational interaction.However, observational evidence has not yet answered the basic questions of the nature, properties, and other interactions of DM.Also, in the absence of any hints of complexity and interactions with and within the dark sector, the general consensus deems DM as particles of a fundamental nature beyond our known form of matter profoundly established in the standard model (SM) of particle physics.Such a DM is required to be electrically neutral, stable at the cosmological time scale, and non-relativistic at the time of matter-radiation equality to permit structure formation.Furthermore, the observation of the cosmic microwave background (CMB) established that DM holds a 27% share of the total energy budget of the present universe, corresponding to a relic density Ωh 2 ≃ 0.12 [2,3].
With a decade-old theoretical development in modeling different particle DM candidates, the extensive possibility of a wide mass range of viable DM candidates has been opened [4].This, in turn, pushed the community to develop ingenious approaches for the next generation of DM experiments, covering ultralight to very heavy ones.Among this wide range of DM mass, the model-independent lower limit of around 10 −22 eV comes from the de Broglie wavelength of bosonic DM that can be confined within a dwarf galaxy [5][6][7].This paradigm of ultralight fuzzy DM candidates has received significant attention lately due to its ability to alleviate tensions between anomalies ranging from theoretical modeling and simulation structure formation and vast astrophysical data.Further refinement of the DM mass bounds is influenced by the specific properties and nature of DM.Notably, the lower limit for the DM mass is more constrained in the case of fermionic DM, because of the Pauli exclusion principle.Furthermore, observations of the Lyman-α forest play a crucial role in establishing constraints on DM.These observations set a lower limit of around 5.3 keV [8] for warm DM by constraining its free streaming length.In contrast, a model-independent upper limit of the DM mass can be derived in the range of 10 3 solar mass from the stability of stellar clusters in galaxies [9,10].Once again, observations of the Lyman-α forest play a role in further constraints from Poisson noise [11].
The consideration of a specific DM production paradigm in the early stage of the universe may further constrain the mass range for a viable DM candidate.For instance, the number-changing pair annihilation of DM to SM particles determines its present mass density, where it maintains the chemical and kinetic equilibrium with the thermal soup in the early universe.Interestingly, the requirement of the unitarity of the S-matrix sets a modelindependent upper bound on the DM mass for this scenario [12][13][14].The implication of unitarity offers the maximum inelastic cross section, which fixes the minimum number density of the frozen-out DM.Using this number density, one can establish the maximum allowed DM mass by fulfilling its observed relic density.In the theories of DM with long-range forces, DM bound states can form which lessen the unitarity bound than the picture without bound states, since the formation of bound states reduces the effective inelastic annihilation rate [15][16][17].In addition, the dark sector with particle-antiparticle asymmetry enforces a nonzero equilibrium chemical potential for DM, which further constrains the unitarity limits by demanding an increased effective DM number density at the time of freeze-out [15,18].Furthermore, different indirect searches for DM may put a lower limit on the DM mass for some specified scenarios.A strong model-independent lower bound for thermal DM that is annihilating to visible states through an s-wave process is about 20 GeV [19].In addition, a more restrictive lower limit has recently been found.It has been shown that the lower bound is 200 GeV, considering H.E.S.S. and other updated observational data [20].
In particular, all the DM scenarios mentioned so far pay attention to the 2 → 2 numberchanging process where a DM pair annihilates into a pair of SM particles, that is, the Weakly Interacting Massive Particle (WIMP) paradigm [21,22]. 1 Moreover, it is not necessary that the number-changing processes involve SM particles, so they may also occur within the dark sector.The minimalist realization of this scenario is the 3 → 2 process, where this kind of number-changing reaction involves a single DM species.In general, such processes arise in DM theories with new sizable self-interactions, and in several contexts as self-interacting DM [29][30][31], the Strongly Interacting Massive Particle (SIMP) paradigm [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49], or even the ELastically DEcoupling Relic (ELDER) scenario [50,51].Now, the question arises: What would be the consequences of the unitarity of the Smatrix for a general DM number-changing process of the kind r → n with r > n ≥ 2? Here, we may focus on a subset of the reaction with the type r → 2 instead of the generalized r → n interaction to avoid the complexity of handling the partial-wave decomposition for the r-body initial state.However, the partial-wave analysis for a two-body initial state is simple.The thermally averaged cross section ⟨σ r→2 v r−1 ⟩ can be easily obtained from ⟨σ 2→r v⟩ utilizing the fact that both cross sections are related in thermal equilibrium.The implication of unitarity helps to calculate the maximum inelastic cross section for the 2 → r process once the total cross section, calculated using the optical theorem and the elastic scattering cross section for the 2 → 2 process are known.
It is essential to mention that the early history of the universe plays a crucial role in the genesis of DM, since the decoupling of thermal DM occurred at that time.Generally, the studies of DM consider the standard cosmological picture in which the radiation energy density is assumed to dominate the energy budget before the Big Bang nucleosynthesis (BBN).However, there is no direct evidence for the energy content at very high temperatures.Therefore, it is vital to look at the effects of modified cosmology on the production of DM.In recent times, the evolution of the DM in the period of non-standard expansion usually triggered by the decay of a long-lived massive particle [48,[52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71] or by Hawking evaporation of primordial black holes  is receiving increasing attention.2 All such studies point towards the fact that nonstandard cosmology alters the value of the thermally averaged cross section needed to satisfy the observed relic of DM.Such a modification in the thermally averaged cross section may also change the unitarity mass bound of the DM. Rcently, the impact of early matter domination on unitarity limits has been studied [14].
The present work explores how such bound modifies in the context of two different non-standard cosmological pictures.One is the kination-like scenario in which a species, ϕ, dominates the energy density in the pre-BBN era.Here, the energy density of ϕ maintains the following redshift behavior, ρ ϕ ∝ a −(4+n) with n > 0 where a is the scale factor [112].Another is late-time reheating, where delayed reheating occurs after inflation due to the suppressed interaction of the inflaton field [52].In general, it can be realized as a scenario of fast expansion with entropy injection.We demonstrate that the kination-like scenario restricts the unitarity limits compared to the standard case.Since one needs a larger cross section than a standard picture due to the early freeze-out of DM.However, the picture is interestingly opposite in the case of late-time reheating.Although the DM decouples earlier because of the fast expansion, the unitarity bounds are relaxed here because of the entropy injection after the decoupling of DM from the thermal soup.
This article is decorated as follows.In Sections 2 and 3, we present the detailed derivation of the maximum thermally-averaged cross section allowed by the unitarity of the Smatrix.We discuss two different non-standard cosmological pictures: kination-like and latetime reheating in Section 4. Section 5 shows the analytical expressions for freeze-out and cross sections for the radiation-dominated universe and the mentioned modified cosmologies, and we also demonstrate our results.Finally, we summarize our findings in Section 6.
2 S-matrix: Unitarity and its Consequences S-matrix theory gives the information about the probability amplitude of various outcomes of a scattering process [113,114].Interestingly, the unitarity of the S-matrix dictates the conservation of the probability of interaction.Although the unitarity of the S-matrix appears to be apparently simple, it has some remarkable implications.We would recollect some consequence of it in the case of DM annihilation, since our primary interest is to put a bound for DM scattering.Here, we consider the multiparticle momentum eigenstate normalization convention [13] and closely follow the derivation in Ref. [14].One of the outstanding outcomes of the S-matrix is the optical theorem, which relates the imaginary part of the scattering amplitude (A ββ ) with the total cross section (σ total ) and can be expressed in the center-ofmass frame as with β describes the two-particle initial state.Here, | ⃗ k| and E cm represent the magnitude of the three momenta of each particle in the initial state and the total energy of the initial state in the CM frame, respectively.The state α is delineated by the three momenta, the z component of the spin (or helicity) of the individual particles along with all other internal quantum numbers, and σ tot ≡ γ σ β→γ where γ corresponds to the possible final states.It is important to mention that our focus will be on the interaction of nonidentical scalar particles for simplicity, and generalization to other spin can be done following Refs.[13,113,114].Here, we would work on a general basis state instead of the basis of momentum eigenstates, which can be done as a result of rotational invariance.
One needs to know the amplitude first to calculate the cross section, and to do so, let us focus on the operator of the S-matrix (S).It is obvious that the initial and final states are identical in the absence of any interaction.Therefore, the operator S is simply the identity 1, as nothing has happened.If an interaction occurs, one needs to subtract the identity from S, where all the information about the interaction is encoded in the nontrivial part S − 1.Now, we can express the operator of the S-matrix using the definition of amplitude as where k β and k γ represent total four momenta.Furthermore, the S-matrix can be expressed as a function of the quantum number of the total angular momentum l, the z component of the total angular momentum quantum number m and the total energy E tot .Using the relation of Eq. (2.2), the general expression of the matrix element for 2-to-2 scatterings where q and q ′ represent the channel indices for the two particle species.Here, Y lm ( k′ 1 ) and Y lm ( k1 ) correspond to the spherical harmonics in the direction of the momentum unit vector k′ 1 and k1 .Then, the total cross section, the sum of the elastic and inelastic cross sections, in a given channel can be obtained with the help of the optical theorem in Eq. (2.1) and the expression of the matrix element in Eq. (2.3) Then, we obtain the elastic scattering cross section σ el in the CM frame from the channel q after performing the integration over dΩ( k′ 1 ) while choosing k1 along the z direction Now, one can easily obtain the inelastic cross section by subtracting the elastic in Eq. (2.5) from the total cross sections in Eq. (2.4) Finally, we obtain the upper limits on the total inelastic cross section taking into account that The above expression can also be expressed in terms of velocity using the fact that, in the nonrelativistic limit, the momentum can be expressed as , where v refers to the relative velocity of the annihilating particles.Finally, we note that for a collision among identical particles, an additional multiplicative factor of 2 is required in the cross section to avoid double counting [113].

Dark Matter Annihilation and Unitarity Bound
We have obtained the upper bound on the inelastic cross-section, and now we are interested in calculating the unitarity bound on the inelastic reaction rate with the help of the Boltzmann equation for DM.The Boltzmann equation for the number-changing DM annihilation of the i th particle in the r → 2 process, where r ≥ 2, i.e., (1, 2, ) in the isotropic and homogeneous expanding early universe, can be manifested as with ∆n i being the net change in the i th particle number in the reaction, and where the thermally-averaged cross section is given by where σ r→2 represents the cross section, summed over final spins and averaged over initial spins, for the r → 2 process, v is the relative velocity between each particle pair, and f eq i refers to the equilibrium distribution function of the i th particle species.More specifically, n eq 1 n eq 2 • • • n eq r ⟨σ r→2 v r−1 ⟩ is the equilibrium rate density for the annihilation.Similarly, one can perform the thermal average to obtain ⟨σ 2→r v⟩.Using the relation of detailed balance, n eq 1 n eq 2 • • • n eq r ⟨σ r→2 v r−1 ⟩ = n eq a n eq b ⟨σ 2→r v⟩ (true for any individual process in equilibrium), Eq. (3.1) boils down to the following simplified form Using the maximum value of the inelastic cross section from Eq. (2.7), one can calculate the maximal value ⟨σ 2→r v⟩ max for the thermal average integral in Eq. (3.2) for the 2 → r process, considering equal masses for all r + 2 particles participating in the interaction, as4 [14] where x ≡ m/T and T is temperature of the SM bath.It is evident that ⟨σ r→2 v r−1 ⟩ max contains an exponential suppression factor, e −(r−2)x , which reveals the expense of phase space to produce each additional particle for r ≥ 3. Now, we can obtain the maximum thermally averaged rate for the r → 2 process using Eq.(3.4) and the detailed balance equation, as [14] where g stands for the internal degrees of freedom of the DM.Therefore, for r = 2 and l = 0 (s-wave), the general expression of the thermally-averaged cross section (for 2 → 2 process) brings us to the familiar result [12] Similarly, the maximum value of the thermally averaged s-wave annihilation cross section for 3 → 2 can be written as Note that the maximum inelastic cross section for identical initial-state particles is twice that of the non-identical scenario.However, a multiplicative factor 1/2 must be added as a symmetry factor when performing the thermal averaging integral in ⟨σ 2→2 v⟩ for identical particles in the initial state.As a result, both the expressions in Eqs.(3.4) and (3.5) are valid for identical and non-identical initial-state particles.Now, we pause to mention that the expressions in Eqs.(3.6) and (3.7) can be used to put a unitarity bound on the DM mass, as will be shown in Section 5.

Low-temperature Reheating
In the standard cosmological paradigm, between the end of inflationary reheating and the beginning of BBN at T = T BBN ≃ 4 MeV [116][117][118][119][120], the energy density of the universe is dominated by SM radiation with an energy density ρ R given by where T corresponds to the temperature of the SM bath.It follows that the Hubble expansion rate H is therefore with M P ≃ 2.4 × 10 18 GeV being the reduced Planck mass.Additionally, the conservation of the entropy S ≡ s a 3 of the thermal bath, with a being the cosmic scale factor, and the SM entropy density implies that the temperature of the SM bath scales as Here, g ⋆ and g ⋆s are the numbers of relativistic degrees of freedom that contribute to the energy density of the SM and the entropy of the SM, respectively [121].However, it is interesting to emphasize that the standard cosmological scenario is not guaranteed and that alternative cosmologies could also have occurred [122].In the following, we focus on cases characterized by low-temperature reheating. 5This reheating could correspond to the period just after the end of inflation, or to a secondary period in which an extra component beyond SM radiation dominated the energy density of the universe.In particular, two scenarios will be reviewed: one where the extra component ϕ that dominated the expansion of the universe has an energy density that gets diluted faster than radiation and does not decay (that is, a kination-like scenario), and the other where ϕ scales as non-relativistic matter and decays into SM particles (that is, an early matter-dominated scenario).These two scenarios will be described below.

Kination-like
In this scenario, the universe was dominated by a component ϕ whose energy density redshifts faster than free radiation [112], as with n > 0. A typical example of this scenario corresponds to kination [126,127], where n = 2.However, larger values for n are also possible, appearing, for example, in the context of ekpyrotic [128,129] or cyclic scenarios [130][131][132][133]; see also Ref. [134].Interestingly, as this component naturally tends to become subdominant, it is not mandatory to enforce its decay.Let us call T rh the SM bath temperature at which equality ρ R (T rh ) = ρ ϕ (T rh ) occurs and from which the universe is dominated by SM radiation (that is, the standard cosmological scenario is recovered).The Hubble expansion rate is therefore where we have taken into account that, as ϕ is not decaying, the SM entropy is conserved, and therefore the SM temperature follows the standard scaling shown in Eq. (4.4).

Early matter domination
Alternatively, the universe could have being dominated instead by a component ϕ with an energy density that scales like non-relativistic matter: ρ ϕ (a) ∝ a −3 .As this component tends to dominate the total energy density of the universe, it has to decay.Again, T rh is the bath temperature when the equality ρ R (T rh ) = ρ ϕ (T rh ) occurs, defining the beginning of the SM radiation dominance era. 6Taking into account that the decay of ϕ gives rise to a nonadiabatic epoch, the SM temperature scales as [52,64,135] for T ≥ T rh , g ⋆s (T rh ) g ⋆s (T ) It follows that the Hubble expansion rate is Interesting, the SM thermal bath can reach temperatures higher than T rh , up to T = T max [52].
Having settled the evolution of the background, in the next section the dynamics of the thermal DM in such alternative cosmological scenarios, and in particular the impact on the unitarity limit, will be carefully studied.

Freeze-out with a Low-temperature Reheating
In this section, two cases for the DM freeze-out are considered.The first corresponds to the visible freeze-out, where a couple of DM particles annihilate into a couple of SM states, with a total thermally-averaged annihilation cross section ⟨σv⟩.The evolution of the DM number density n can be described with the Boltzmann equation [21] where the DM number density at equilibrium n eq is for non-relativistic DM particles, with g and m being the number of internal degrees of freedom and the mass of the DM particle, respectively.Alternatively, DM could have been generated through a dark freeze-out (that is, a cannibalization process), where r DM particles annihilate into two DM particles with an interaction given by ⟨σ r→2 v r−1 ⟩.In that case, the evolution of n is given by [32,39,44] where, ∆n represents the net change of the DM number in the reaction.It is important to recall that in the present analysis, s-wave annihilations are always considered, which implies that ⟨σv⟩ and ⟨σ r→2 v r−1 ⟩ are temperature-independent quantities.The freeze-out temperature T fo corresponds to the temperature at which DM exits chemical equilibrium.It is important to note that here we assume that the kinetic equilibrium between the dark and visible sectors is maintained at least until the end of chemical freeze-out, so that at freeze-out the two sectors share the same temperature.The freeze-out temperature can be estimated by the equality between the Hubble and the interaction rates or in the case of a visible or dark freeze-out, respectively.

Kination-like
In this case, as the SM entropy is conserved, it is convenient to rewrite Eqs.(5.1) and (5.3) in terms of the comoving yield Y (T ) ≡ n(T )/s(T ) and x ≡ m/T as (5.6) where Y eq ≡ n eq /s.
The freeze-out temperature T fo = m/x fo can be estimated by comparing Eqs.(4.6) with (5.4) or (5.5), and can be expressed in the convenient way (5.8) with x rh ≡ m/T rh , and where W −1 correspond to the branch −1 of the Lambert W function.In the case of a visible freeze-out r = 2, one has to use that ⟨σ r→2 v r−1 ⟩ = ⟨σv⟩, while in the radiation-dominated era n = 0.
In the following, Eqs.(5.6) and (5.7) will be analytically solved in the context of a kination-like cosmology.For convenience, we start with the case corresponding to the dark freeze-out.

Dark freeze-out
If the freeze-out occurs during the radiation-dominated era, Eq. (5.7) can be analytically solved, from the DM freeze-out until today (i.e., small temperature and therefore large x) where Y 0 is the DM yield at late times, and we have used the fact that Y fo ≫ Y 0 .It follows that the cross section ⟨σ r→2 v r−1 ⟩ is (5.10) To match the whole observed DM relic density, it is required that with ρ c ≃ 1.05 × 10 −5 h 2 GeV/cm 3 being the critical energy density, s 0 ≃ 2.69 × 10 3 cm −3 the present entropy density [136], and Ωh 2 ≃ 0.12 the observed DM relic abundance [2].Alternatively, if the freeze-out happens during reheating (5.12) the integral has been split into two pieces, to emphasize the two regimes of H in Eq. (4.6).Therefore As expected, Eq. (5.10) can be recovered from Eq. (5.13) by taking n = 0.

Visible freeze-out
The case of the visible freeze-out in Eq. (5.6) can be computed following the same procedure presented in the previous subsection.However, one could also derive it by fixing r = 2 in Eqs.(5.10) and (5.13), which gives for the freeze-out in the radiation-dominated era, or Figure 1.Kination-like.Required freeze-out temperature T fo = m/x fo (left) and cross section ⟨σv⟩ (right) for 2 → 2 annihilations to fit the observed abundance of DM, in the case of radiation domination (n = 0) and low-temperature reheating with n = 2, 4 and 6.The blue curves correspond to T rh = 10 −2 GeV, while the red curves correspond to T rh = 10 2 GeV.In the red band the freeze-out is relativistic, whereas in the gray band unitarity is violated.during reheating.
The required freeze-out temperature and thermally-averaged annihilation cross section to fit the observed DM abundance, in the case of standard cosmology (n = 0) and a kinationlike scenario with n = 2, 4 and 6, are shown in Fig. 1.The blue curves correspond to T rh = 10 −2 GeV, while the red curves correspond to T rh = 10 2 GeV.In the case of radiation domination, the usual results are recovered: x fo ∼ O(25) and ⟨σv⟩ ∼ O(10 −9 ) GeV −2 , which correspond to a few picobarns, for WIMPs in the GeV ballpark [137].Alternatively, in kination-like scenarios, the universe expands faster, and therefore an early freeze-out is required.High values of n induce small values of x fo , and a potential relativistic freeze-out (i.e.x fo ≲ 3), depicted as a red band in the left panel.Additionally, larger cross sections ⟨σv⟩ are also required, potentially in tension with the unitarity constraint; cf.Eq. (3.6), shown with a gray band in the right panel.It is worth mentioning that, instead of a single constraint, there is a constraint per cosmological scenario, in this case for each choice of n and T rh .However, in Fig. 1 all the different constraints basically overlap.
The impact of the unitarity constraint in the case of a 2 → 2 scattering becomes more clear in Fig. 2, where the gray bands show the excluded regions by unitarity following Eq.(3.6), for n = 0, 2, 4 and 6, in the plane [m, T rh ].The standard upper bound of ∼ 1.3 × 10 5 GeV for the DM mass becomes tighter and depends on T rh , if the freeze-out happens in a kination-like era, as expected from Fig. 1.For example, if T rh = T BBN , the upper bound on the DM mass can be as strong as ∼ 3 × 10 3 GeV or ∼ 1.2 × 10 3 GeV for n = 2 and n = 4, respectively.Interestingly, for the case n = 6, DM may not reach chemical equilibrium in very low-temperature reheating scenarios and therefore freeze-out cannot occur, corresponding to the red dashed line.Finally, the red band corresponding to T rh < T BBN is in tension with BBN.
Additionally, the details for the dark freeze-out through 3 → 2 scatterings and the unitarity bound following Eq.(3.7) are shown in Fig. 3.In this case, smaller DM masses, in the MeV ballpark, are required.Therefore, the impact of low-temperature reheating is much milder, once the BBN bound on T rh is imposed.In a radiation-dominated scenario, the unitarity bound implies m ≲ 1 GeV; however, if freeze out occurs in a kination-like epoch, the bound becomes more stringent, implying that if T rh = T BBN , m ≲ 7 × 10 −1 GeV, m ≲ 5 × 10 −1 GeV, or even m ≲ 4 × 10 −1 GeV, for the cases with n = 2, n = 4 and n = 6, respectively.We note that the most stringent limits on DM mass occur for the kination-like scenario for T rh = T BBN and large values of n.

Early matter domination
In this case, in order to have a successful reheating, ϕ has to decay into SM particles and hence the SM entropy is not conserved. 7Therefore, instead of the yield Y , it is convenient to use the comoving DM number density N (a) ≡ n(a) × a 3 .Equations (5.1) and ( 5.3) can be rewritten, respectively, as (5.16) (5.17) with N eq ≡ n eq × a 3 .The freeze-out temperature can be estimated by comparing Eqs.(4.8) and (5.4) or (5.5), and is given by x 4 rh m 10−6r M 2 P ⟨σ r→2 v r−1 ⟩ 2 1 11−3r for x fo ≪ x rh , (5.18) where it is always assumed that T max ≫ T fo .Interestingly, Eq. (5.18) is also valid for the case of a visible freeze-out, taking r = 2 and replacing ⟨σ r→2 v r−1 ⟩ by ⟨σv⟩.Furthermore and as expected, if the freeze-out occurs in the standard radiation-dominated era, Eq. (5.8) (n = 0) and Eq.(5.18) (x fo ≫ x rh ), are equivalent.
Next, analytical solutions are presented for Eqs.(5.16) and (5.17) in the context of an early matter-dominated scenario.We will begin with the case corresponding to the dark freeze-out for convenience.

Dark freeze-out
If the freeze-out occurs during the radiation era, the solution of Eq. (5.17), or equivalently of Eq. (5.7), is the one presented in Eq. (5.10).Instead, if it occurs during the reheating period, one has that where N fo ≡ N (a fo ), N 0 is the asymptotic value of N (a) at large values of a, much after the freeze-out, and we have used the fact that N fo ≫ N 0 .Here again, the integral has been spit into the two regimes of H in Eq. (4.8).The thermally-averaged cross-section is, therefore, (5.20)

Visible freeze-out
If the freeze-out occurs during radiation domination, the solution of Eq. (5.16) is the same as the one of Eq. (5.14).Alternatively, if it occurs during reheating, one has instead ⟨σv⟩ ≃ 9 2π 5 2 simply corresponding to the limit r = 2 of Eq. (5.20).
The left panel of Fig. 4 displays the freeze-out temperature needed to fit the observed relic of DM as a function of its mass for radiation-dominated universe (thick black solid line) and reheating temperatures, T rh = T BBN (black solid lines), 10 0 GeV (blue dotted line), 10 3 GeV (blue dot-dashed line), 10 6 GeV (blue dashed line) where the DM number changing process is 2 → 2, that is the visible WIMP mechanism.The dilution by large entropy injection has to be overcome by the overproduction of DM at freeze-out, which results in an earlier freeze-out and therefore a small x fo .Interestingly, chemical equilibration of DM requires x fo ≳ 6.5, which corresponds to the red band.Smaller values for x fo could lead to the observed abundance of DM, but not through the WIMP mechanism, but instead through the FIMP paradigm [138][139][140][141][142][143][144].An additional red band corresponds to x fo ≤ 3 and represents the relativistic freeze-out.The right panel of Fig. 4 contains equivalent information, but is now presented in the plane [m, ⟨σv⟩].High freeze-out temperatures correspond to earlier decouplings and, therefore, to lower cross sections.Again, ⟨σv⟩ is bounded from below by the requirement of reaching chemical equilibrium and corresponds to the "No freeze-out" region.It is important to note that the cross section decreases with increasing DM mass for a fixed T rh .The reason is that the freeze-out of heavier DM occurs earlier (that is, at a larger T fo ), demanding a smaller ⟨σv⟩.Additionally, the gray-shaded region corresponding to high values of ⟨σv⟩ is disallowed by unitarity following Eq.(3.6). 8In low-temperature reheating scenarios with large injection of entropy, thermally averaged cross sections much In the bottom panel of Fig. 4, we have shown the allowed parameter space in the reheating temperature and DM mass plane where the white space is free of any constraints.The red band dictates that T rh lower than T BBN is in disagreement with the cosmological observation.The red and gray-shaded region displays the "No freeze-out" and the unitarity constraints discussed earlier.Note that the unitarity constraint is independent of the reheating temperature for m ∼ 1.3 × 10 5 GeV, reflected by the straight vertical shape in the unitarity bound.This is the DM mass where the cross section needed to fit the observed DM abundance (in a radiation-dominated background) and the maximum allowed cross section by unitarity match exactly.However, low-temperature reheating opens up the possibility of higher masses.Interestingly, masses higher than m ∼ 1.3 × 10 5 GeV are severely constrained from above due to unitarity and from below due to the requirement of chemical equilibrium, resulting in a narrow allowed region for T rh .
For completeness, Fig. 5 depicts results equivalents to the ones shown in Fig. 4 but now for 3 → 2 annihilations.In this case, the presence of late-time reheating opens up the DM mass to ∼ 10 8 GeV with ⟨σ 3→2 v 2 ⟩ ∼ 10 −36 GeV −5 .

Summary and Conclusion
The requirement of the de Broglie wavelength of dark matter (DM) to hold it inside galaxies and the stability of the stellar cluster in galaxies collectively put a broad allowed mass range for DM by providing the lower and upper bound of DM mass, respectively.Specifying some properties of DM can further tighten the mass range.Interestingly, one can put the model-independent upper bound by specifying the thermal production of DM in the early universe.The observed DM abundance and unitarity of partial waves from the scattering matrix jointly place an upper limit on DM mass.Primarily, the upper limits on the inelastic cross section for a general number-changing process 2 → r can be derived with the help of the optical theorem, the matrix elements and the elastic scattering cross section for the process.After that, one can obtain the thermally averaged cross section for the r → 2 process by invoking the principle of detailed balance.Finally, the unitarity bounds on the thermally averaged cross section translate to the upper limits on the DM mass satisfying the relic density constraints.It is known that the maximum allowed DM mass for the 2 → 2 and 3 → 2 DM annihilation processes is around 130 TeV and 1 GeV, respectively.However, these bounds depend not only on the particle physics model, but also have a strong dependence on the cosmological evolution of the universe, being valid only if the universe followed the so-called "standard cosmological scenario".
Instead, this article explores the DM mass bound in non-standard cosmological setups characterized by low-temperature reheating.In particular, we focus on i) kination-like scenarios, where the early universe was dominated by a fluid with an energy density that gets diluted faster than free radiation, and ii) early matter-dominated scenarios, where a component with an energy density that scales as nonrelativistic matter dominates the early universe and eventually decays into SM particles.
First, we study the kination-like universe, which demands a larger thermally-averaged annihilation cross section to saturate the observed abundance of DM compared to the standard radiation-dominated picture since, in this case, freeze-out occurs early.As a result, the upper bound on the DM mass becomes more stringent than in the standard case.For example, if the reheating temperature is as low as a few MeVs (corresponding to the start of the Big Bang nucleosynthesis epoch), the usual bound of the DM mass m ≲ 130 TeV can be reduced to a few TeVs for WIMPs.
Second, we also consider the picture of early matter domination, which dictates a fast expansion with entropy production.Although the DM freezes out early in this scenario, here one needs a smaller thermally averaged cross section to feature an observed relic of DM than a radiation-dominated picture to compensate for the dilution of the relic due to the huge entropy injection.Thus, the presence of an early matter domination relaxes the mass bound, and eventually the allowed mass range of the DM is enhanced.Here, the usual bound of the DM mass m ≲ 130 TeV can be relaxed, making the WIMP DM masses up to ∼ 10 10 GeV viable.
Before closing, we want to emphasize that the evolution of the early universe is largely unknown.The standard assumption of a universe dominated by standard-model radiation from the end of cosmological inflation until matter-radiation equality, together with a transition from an inflaton-dominated to a radiation-dominated universe occurring at a very early time, cannot be taken for granted.Having that in mind, here we have studied the impact of the unitarity bound on DM in the case of low-temperature reheating scenarios.

Figure 4 .
Figure 4. Early matter domination.The same as in Figs. 1 and 2, for 2 → 2 annihilations, but for early matter domination.Additionally, the "No freeze-out" constraint is shown by the red-shaded region.

Figure 5 .
Figure 5. Early matter domination.The same as in Fig. 4, but for dark freeze-out through 3-to-2 annihilations.