Real Scalar Dark Matter: Relativistic Treatment

TL;DR

This paper presents a fully relativistic analysis of a real scalar dark matter candidate, exploring its thermalization, freeze-out regimes, and observational constraints, highlighting its potential to be warm or cold dark matter.

Contribution

It introduces a relativistic treatment of scalar dark matter evolution, including a new Bose-Einstein analog of the Gelmini-Gondolo formula and comprehensive parameter space analysis.

Findings

Dark matter can be warm or cold depending on parameters.

Relativistic effects significantly influence freeze-out and thermalization.

New formula for thermal averaged cross section derived.

Abstract

A stable real scalar provides one of the simplest possibilities to account for dark matter. We consider the regime where its coupling to the Standard Model fields is negligibly small. Due to self-coupling, the scalar field can reach thermal or at least kinetic equilibrium, in which case the system is characterized by its temperature and effective chemical potential. We perform a fully relativistic analysis of dark matter evolution, thermalization conditions and different freeze-out regimes, including the chemical potential effects. To this end, we derive a relativistic Bose-Einstein analog of the Gelmini-Gondolo formula for a thermal averaged cross section. Finally, we perform a comprehensive parameter space analysis to determine regions consistent with observational constraints. Dark matter can be both warm and cold in this model.

Equations

$$V=\frac{m^2}{2}{S}^2+\frac{\lambda }{4!}{S}^4,$$

$$f(p)=\frac{1}{e^{\frac{E-\mu }{T}}-1},$$

$$n=\frac{1}{2{\pi }^2}{\int }_m^{\infty }\frac{(E^2-m^2{)}^{1/2}}{e^{\frac{E-\mu }{T^{\prime }}}-1}EdE,$$

$$V_{\varphi \mathrm{h}}=\sigma \varphi h^2,$$

$$V_{\varphi \mathrm{S}}=\frac{1}{2}{\lambda }_{\varphi S}{\varphi }^2{S}^2,$$

$$V_{\inf }=\frac{1}{2}{m}_{\varphi }^2{\varphi }^2,$$

$$n(a=1)=c{m}_{\varphi }^3,$$

$$E(a)=E(1)/a,n(a)=n(1)/a^3.$$

$$\Gamma =\frac{{\sigma }^2}{8\pi m_{\varphi }}\sim H\simeq \frac{m_{\varphi }{\varphi }_{\mathrm{end}}}{\sqrt{6}{M}_{\mathrm{Pl}}{a}^{3/2}},$$

$$T_R\sim 10^{-1}\sqrt{\Gamma M_{\mathrm{Pl}}},$$

$$\frac{E(a)}{T_R}\sim {\left(\frac{m_{\varphi }{\sigma }^2{M}_{\mathrm{Pl}}}{{\varphi }_{\mathrm{end}}^4}\right)}^{1/6}\ll 1$$

$$\frac{E(a)}{T_R}\sim \frac{T}{T_{\mathrm{SM}}}\ll 1.$$

$$\frac{n(a)}{n(T(a))}\sim \frac{c{\pi }^2}{\zeta (3)}.$$

$$\frac{dn}{dt}+3Hn=2({\Gamma }_{2\to 4}-{\Gamma }_{4\to 2}),$$

$${\Gamma }_{a\to b}=\int \left(\prod _{i\in a}\frac{d^3{\mathbf{p}}_i}{(2\pi {)}^{3}2E_i}f(p_i)\right)\left(\prod _{j\in b}\frac{d^3{\mathbf{p}}_j}{(2\pi {)}^{3}2E_j}(1+f(p_j))\right)|{\mathcal{M}}_{a\to b}{|}^2(2\pi {)}^4{\delta }^4(p_a-p_b).$$

$$f(p)=\frac{1}{{\exp }^{\frac{E-\mu }{T}}-1},$$

$$\displaystyle f(k_{1})~{}f(k_{2})~{}f(k_{3})~{}f(k_{4})~{}\bigl{(}1+f(p_{1})\bigr{)}~{}\bigl{(}1+f(p_{2})\bigr{)}=$$

$$\displaystyle f(p_{1})~{}f(p_{2})~{}\bigl{(}1+f(k_{1})\bigr{)}~{}\bigl{(}1+f(k_{2})\bigr{)}~{}\bigl{(}1+f(k_{3})\bigr{)}~{}\bigl{(}1+f(k_{4})\bigr{)}$$

$$\frac{d}{dt}n{a}^3=0.$$

$$3Hn\gtrsim 2{\Gamma }_{2\to 4},2{\Gamma }_{4\to 2},$$

$$\sigma (p_1,p_2)=\frac{1}{4F(p_1,p_2)}\int |{\mathcal{M}}_{2\to 4}{|}^2(2\pi {)}^4{\delta }^4(p_1+p_2-\sum _i{k}_i)\prod _i\frac{d^3{\mathbf{k}}_i}{(2\pi {)}^{3}2E_{k_i}}(1+f(k_i)),$$

$$f(p)=\frac{1}{e^{\frac{u\cdot p}{T}}-1},$$

$${\Gamma }_{2\to 4}=(2\pi {)}^{-6}\int d^3{\mathbf{p}}_1{d}^3{\mathbf{p}}_{2}f(p_1)f(p_2)\sigma (p_1,p_2)v_{\mathrm{M}\text{o}\mathrm{l}},$$

$$v_{\mathrm{M}\text{o}\mathrm{l}}=\frac{F(p_1,p_2)}{E_1{E}_2}.$$

$$p=\frac{p_1+p_2}{2},k=\frac{p_1-p_2}{2},$$

$${d^{3}{\bf p_{1}}\over 2E_{1}}{d^{3}{\bf p_{2}}\over 2E_{2}}=2^{4}d^{4}p~{}d^{4}k~{}\delta\Big{(}(p+k)^{2}-m^{2}\Big{)}~{}\delta\Big{(}(p-k)^{2}-m^{2}\Big{)}\;.$$

$$p^0=E\cosh \eta ,$$

$$p^1=E\sinh \eta \sin \theta \sin \varphi ,$$

$$p^2=E\sinh \eta \sin \theta \cos \varphi ,$$

$$p^3=E\sinh \eta \cos \theta ,$$

$$d^{4}p={\sinh }^2\eta E^{3}dEd\eta d{\Omega }_p,$$

$$\int \frac{d^3{\mathbf{p}}_1}{2E_1}\frac{d^3{\mathbf{p}}_2}{2E_2}G(p_1,p_2)=2{\int }_m^{\infty }dE\sqrt{E^2-m^2}{E}^2{\int }_0^{\infty }d\eta {\sinh }^2\eta \int d{\Omega }_{p}d{\Omega }_{k}G(p_1,p_2),$$

$$1+f(k_i)=1+\frac{1}{e^{(k_i^0\cosh \eta +k_i^3\sinh \eta -\mu )/T}-1}$$

Full text

1. 1 Introduction

2. 2 Inflationary motivation

3. 2.1 Initial conditions for dark matter evolution

4. 2.2 Example

5. 3 Boltzmann equation and reaction rates

6. 3.1 Relativistic reaction rates with the Bose–Einstein distribution function

7. 3.2 Conditions for thermal and kinetic equilibria

8. 4 Freeze–out

9. 4.1 Non–relativistic freeze–out

10. 4.2 Relativistic freeze–out

11. 4.3 Ultra–relativistic freeze–out

12. 5 Parameter space analysis

13. 5.1 Thermalized dark matter

14. 5.2 Dark matter in kinetic equilibrium and $\mu\not=0$

15. 6 Conclusion

1 Introduction

The nature of dark matter (DM) remains one of the outstanding questions of modern physics. The popular WIMP (weakly interacting massive particle) paradigm currently finds itself under pressure due to ever improving bounds on dark matter interaction with nucleons [1]. This motivates one to explore alternative scenarios in which the dark matter relic density is not related to the annihilation into Standard Model (SM) particles. In this work, we consider in detail one such example. In this case, the DM relic density is set by either the dark sector thermodynamics or postinflationary initial conditions.

We study the simplest dark sector consisting of a single real scalar field $S$ with the potential [2]

[TABLE]

where the self–coupling $\lambda$ is large enough for the $S$ quanta to reach kinetic equilibrium at least. We assume that the dark matter coupling to the Standard Model, $\lambda_{HS}H^{\dagger}HS^{2}$, is negligibly small, $\lambda_{HS}<10^{-13}$. In this case, the observable and dark sectors do not equilibrate, DM freeze–in production [3] is suppressed and the DM evolution is determined entirely by its self–coupling [4].

Although we focus on just one model of dark matter, the novelty of our approach is that we account for all the stages in the (thermal) dynamics, which requires a relativistic treatment with the Bose–Einstein distribution function. In particular, we derive relativistic expressions for the number conserving and number changing reaction rates generalizing the well known Gelmini–Gondolo result [5], and solve the corresponding Boltzmann equation coupled with the entropy conservation condition. This allows us to study the different freeze–out regimes, including the relativistic one, as well as trace the transition from the pre-freeze–out epoch to post-freeze–out by deriving the chemical potential evolution. We consider both fully thermalized dark matter and DM in kinetic equilibrium. Finally, we delineate parameter space consistent with all the observational constraints, for which a relativistic treatment is essential.

Regarding observational prospects, dark matter detection would be very challenging due its negligible interaction with the visible sector. One may potentially observe effects due to its self–interaction, but even that is not guaranteed because $\lambda$ is allowed to be very small and still consistent with the relic density constraints.

Previous analyses of this or closely related models include Refs. [6]-[11] where $S$ is treated as a WIMP–like particle, and non–relativistic studies of feebly interacting $S$ [12]-[16]. We go beyond these results in that we offer a fully relativistic treatment and also include the chemical potential. Relativistic effects in the Boltzmann equation within other contexts have been considered, for example, in [17]-[19]. These analyses are not however applicable to the model at hand.

The paper is organized as follows: after providing an inflationary motivation for our study, we derive the expressions for the relativistic reactions rates. We then apply them to derive the thermalization and freeze–out conditions, and delineate parameter space consistent with the observed DM relic abundance.

2 Inflationary motivation

In this section, we motivate the setting for our study. Our main assumptions are that the interaction between the SM fields and dark matter (DM) is negligible, and that DM reaches some degree of thermal equilibrium, be it kinetic or chemical one. We argue that the temperatures of the observable and dark sectors can differ dramatically, and that the latter may also be characterized by effective chemical potential.

2.1 Initial conditions for dark matter evolution

A bosonic system in thermal or kinetic equilibrium is described by its temperature $T$ and effective chemical potential $\mu$ via the momentum distribution

[TABLE]

where $E=\sqrt{m^{2}+\vec{p}^{2}}$. Here the chemical potential is not necessarily associated with the conserved Noether current, instead it reflects approximate conservation of a particle number in a given regime. The above distribution applies also to bosonic systems in an expanding FRW Universe [20, 21, 22].

In general, the dark matter temperature $T$ can be very different from that of the observable sector $T_{\rm SM}$, as long as the interaction between the two is very weak [4]. Furthermore, if the number changing processes in the dark sector are inefficient, dark matter can be assigned effective chemical potential $\mu$ which determines its number density at a given temperature,

[TABLE]

where $\mu<m$. Even in the simplest case of real scalar dark matter, both $T$ and $\mu$ are in general necessary to describe its state. These determine the initial conditions for the evolution of the system from the relativistic regime to freeze–out and its eventual state we observe today.

2.2 Example

Dark matter and observable matter may be produced by very different mechanisms resulting in different thermodynamic properties of the two sectors. In particular, their temperatures may differ by orders of magnitude and the hidden sector can be endowed with chemical potential if the number changing processes are inefficient.

Consider a simple possibility that the observed matter and dark matter are produced after inflation via direct interaction with the inflaton $\phi$. Suppose that the leading interaction terms are

[TABLE]

where the unitary gauge for the Higgs field $h$ has been assumed. For the purpose of illustration, let us also choose the simplest inflaton potential,

[TABLE]

with $m_{\phi}\sim 10^{13}\;{\rm GeV}\gg m$. The matter production mechanisms depend on the balance between $\sigma$ and $\lambda_{\phi S}\phi$. For relatively small $\sigma$, SM matter is produced at late times via perturbative decay of the inflaton, while dark matter can be produced efficiently right after inflation through parametric resonance.

Let us consider these mechanisms in more detail. During the preheating stage, the inflaton oscillates with a decreasing amplitude as $\phi(t)\simeq\phi_{0}/m_{\phi}t\;\sin m_{\phi}t$, where $\phi_{0}$ is close to the Planck scale. This creates an oscillating mass term for $S$ and leads to efficient dark matter production as long as $q=\lambda_{\phi S}\phi^{2}/m_{\phi}^{2}\gg 1$. When $q$ reaches 1, the resonance stops and the subsequent DM evolution is determined primarily by its self–coupling $\lambda$. If $\lambda$ is small enough, the total number of the DM quanta remains approximately constant.

For small $\lambda_{\phi S}\lesssim 10^{-8}$, the resonance is efficient for a short time of order 10$m_{\phi}^{-1}$ and produces relativistic DM quanta with typical momentum $k\sim m_{\phi}$ (see [23] for a recent study). It converts only a small fraction of the inflaton energy into radiation. The corresponding DM number density is conveniently parametrized by

[TABLE]

where $c$ depends on the strength of the resonance $q$ and $a$ is the FRW scale factor which we set to one at this stage. The energy of the DM quanta and their density evolve in time as

[TABLE]

As long as the energy density of the Universe is dominated by the inflaton, $a\propto t^{2/3}$.

Radiation domination sets in when perturbative inflaton decay $\phi\rightarrow hh$ rate becomes comparable to the Hubble rate,

[TABLE]

where $\phi_{\rm end}$ is the inflaton amplitude at the end of the resonance. This equation can be solved for $a$ which plays the role of the time variable. The resulting reheating temperature is found via the usual relation

[TABLE]

which assumes that nearly all of the inflaton energy converts into Standard Model radiation.111Here we neglect the $h$ production via tachyonic resonance. This is justified for $\sigma\phi/m_{\phi}^{2}<1$. At this point, the ratio between the average energy of the DM quanta $E(a)\simeq m_{\phi}/a$ and $T_{R}$ is given by

[TABLE]

for $\sigma\ll M_{\rm Pl}$. This ratio stays constant in time since both quantities scale as radiation. Thus, it determines the ratio between the hidden sector temperature $T$ and that of the observable sector $T_{\rm SM}=T_{R}$, once dark matter reaches kinetic equilibrium through self–interaction,

[TABLE]

While the temperature is fixed by the average energy of the quanta, the number density is an independent quantity determined by the production mechanism. It is instructive to compare the DM number density $n(a)$ to the equilibrium density at temperature $T(a)$ with zero chemical potential, $n(T)\simeq\zeta(3)/\pi^{2}\;T^{3}$,

[TABLE]

This quantity also stays constant in time if the number changing interactions are suppressed. Its magnitude depends strongly on the efficiency of parametric resonance. If the resonance is weak, $\lambda_{\phi S}<10^{-8}$, one finds $c\ll 0.1$ and ${n(a)/n(T(a))}\ll 1$.222We thank Stanislav Rusak for verifying this point with lattice simulations. This implies under–density of the DM quanta at a given temperature and thus the presence of negative chemical potential.

The thermalization of dark matter is controlled by the self–coupling $\lambda$. For small $\lambda\ll 1$, e.g. $10^{-6}$ for typical cases, dark matter reaches kinetic equilibrium at some stage, however the number changing processes are suppressed by a further factor of $\lambda^{2}$ and chemical equilibrium never sets in.

Note that integrating out the inflaton leads to a tiny Higgs–DM coupling of order $(\lambda_{\phi S}/8\pi^{2})\sigma^{2}/m_{\phi}^{2}$. At sufficiently small $\sigma$, it is irrelevant to both thermalization and freeze–in production of dark matter [3], making the effect emphasised in Ref. [18] negligible.

This example illustrates that dark matter and observed matter can be produced by very different mechanisms, in which case one expects different temperatures in the two sectors. Furthermore, the dark sector can be endowed with effective chemical potential, as long as the number changing interactions are inefficient.

3 Boltzmann equation and reaction rates

The particle density $n(t)$ evolution is described by the Boltzmann equation. In addition to the Universe expansion, $n(t)$ is affected by the particle number changing processes such as $SS\leftrightarrow SSSS$ (Fig. 1) and those of higher order. Keeping the lowest order terms, the Boltzmann equation in the FRW background reads (see e.g. [22, 12])

[TABLE]

where $H=\dot{a}/a$, the factor of 2 comes from the particle number change in the scattering process and the reaction rates per unit volume are

[TABLE]

Here ${\cal M}_{a\rightarrow b}$ is the QFT transition amplitude, in which we also absorb the initial and final state symmetry factors; $f(p)$ is the momentum distribution function. It can deviate from the corresponding thermal distribution, yet as long as the system enjoys kinetic equilibrium through efficient $2\rightarrow 2$ scattering, $f(p)$ takes the form [20],[22]

[TABLE]

where $\mu$ is the effective chemical potential. This can be understood from the fact that such a Bose–Einstein distribution maximizes the entropy with the approximately constant particle number, while the number changing interactions are relatively slow.

In thermal (chemical) equilibrium, the number changing reaction rates are much greater than the Hubble rate $H$ and $\Gamma_{2\rightarrow 4}$ equals exactly $\Gamma_{4\rightarrow 2}$. This is because at $\mu=0$

[TABLE]

due to energy conservation in the reaction $p_{1}p_{2}\leftrightarrow k_{1}k_{2}k_{3}k_{4}$, and $|{\cal M}_{a\rightarrow b}|=|{\cal M}_{b\rightarrow a}|$. As a result, the right hand side of the Boltzmann equation vanishes and the total number of particles is conserved,

[TABLE]

When the system departs from thermal equilibrium and a non–zero $\mu$ develops, the cancellation is no longer exact. As a result, the total particle number changes. Eventually, due to the Universe expansion, both $\Gamma_{2\rightarrow 4}$ and $\Gamma_{4\rightarrow 2}$ decrease to the extent that the right hand side of the Boltzmann equation becomes negligible again and the total particle number remains approximately constant. This happens roughly when

[TABLE]

which is called “freeze–out”. The freeze–out can take place in both relativistic and non–relativistic regimes, which we will consider separately in what follows.

In order to understand both regimes, we need to derive compact expressions for the rates which can be analyzed either analytically or numerically.

3.1 Relativistic reaction rates with the Bose–Einstein distribution function

Consider the reaction rate $\Gamma_{2\rightarrow 4}$. Following Gelmini and Gondolo [5], we find it convenient to express this rate in terms of the cross section $\sigma(p_{1},p_{2})$,

[TABLE]

where $F(p_{1},p_{2})=\sqrt{(p_{1}\cdot p_{2})^{2}-m^{4}}$ and $1+f(k_{i})$ are the final state Bose–Einstein enhancement factors. In this expression, the momentum distribution function can be put in a Lorentz–covariant form as (see e.g. [24])

[TABLE]

where $u_{\mu}$ is the 4–velocity of our reference frame relative to the gas rest frame in which $u=(1,0,0,0)^{T}$. Apart from the Bose–Einstein factors, the above cross section is manifestly Lorentz–invariant.

The reaction rate is then expressed as

[TABLE]

where the Møller velocity is defined by

[TABLE]

It is clear that the rate is proportional to the thermal average $\langle\sigma v_{\rm M\o l}\rangle$.

The cross section is easiest calculated in the center–of–mass frame.333The public code CalcHEP [25] computes the cross sections either in the center–of–mass or in the lab frame. Thus, for numerical analysis, it is necessary to convert our expressions into one of these frames. Hence, for each pair $p_{1},p_{2}$ in the gas rest frame, we find the center–of–mass frame and compute the corresponding cross section. The 4-velocity factor is thus a function of the momenta, $u_{\mu}(p_{1},p_{2})$.

Introduce

[TABLE]

such that

[TABLE]

The center–of–mass frame is defined by the requirement that $p$ has zero spacial components. Let us parametrize the timelike $p$ as

[TABLE]

where $E$ is the particle energy in the center–of–mass frame, $\eta$ is the rapidity and $\theta,\phi$ are the angular variables (see Appendix A). Then

[TABLE]

where $\Omega_{p}$ is the solid angle in $p$-space. Due to the $\delta$–functions, the $k$–integral reduces to angular integration over the solid angle $\Omega_{k}$ in $k$–space. We thus have, for any $G(p_{1},p_{2})$,

[TABLE]

where in the integrand one must set $k_{0}=0,|{\bf k}|=\sqrt{E^{2}-m^{2}}$ in $k$-dependent quantities. Note that $E$ is half the center–of–mass energy.

The cross section is calculated in the center–of–mass frame. Thus, we need to transform the final state momenta $k_{i}$ to that frame, $k_{i}\rightarrow\Lambda k_{i}$, where $\Lambda$ is the corresponding Lorentz transformation given in Appendix A. Due to Lorentz–covariance, $\sigma(p_{1},p_{2})$ retains the same form except for the Bose–Einstein enhancement factors which now become

[TABLE]

since $(\Lambda^{-1}u)\cdot k_{i}=k_{i}^{0}\cosh\eta+k_{i}^{3}\sinh\eta$ (see Appendix A). As a result, the center–of–mass cross section depends on $\eta$ as well, $\sigma_{\rm CM}(E,\eta)$. The angular integrations can be performed explicitly with the result

[TABLE]

Here, $\sigma_{\rm CM}(E,\eta)$ includes the Bose–Einstein factors for the final state and is non–zero for $E\geq 2m$. The cross section can be computed numerically with CalcHEP by absorbing $1+f(k_{i})$ into a momentum–dependent vertex. Note that in our convention, the symmetry factors due to the identical particles in the initial and final states, namely, $1/(2!4!)$, have been absorbed into $\sigma$. In the small $T$ limit, the final state enhancement factors can be neglected and one recovers the Gelmini–Gondolo result

[TABLE]

where $K_{1}(x)$ is the modified Bessel function and $\mu$ is set to zero.

The rate $\Gamma_{4\rightarrow 2}$ can be obtained from $\Gamma_{2\rightarrow 4}$ by noting that

[TABLE]

such that

[TABLE]

Therefore

[TABLE]

While obtaining $\sigma_{\rm CM}(E,\eta)$ for the $2\rightarrow 4$ and $4\rightarrow 2$ reactions in a closed form does not seem possible, the $2\rightarrow 2$ reaction is simple enough such that $\sigma_{\rm CM}^{2\rightarrow 2}(E,\eta)$ can be computed explicitly. The corresponding reaction rate is needed to describe kinetic equilibrium. We have

[TABLE]

where $|{\cal M}|^{2}=\lambda^{2}/(2!2!)$. As explained above, we include the symmetry factors for the final and initial states in the amplitude, so the large $E$ limit of this cross section differs from the standard result (see e.g. [26]) by $1/2!$. This is admissible since we are only interested in the thermal averages, which are convention–independent. The angular dependence comes entirely from $f(k_{i})$. Computing the integral, we obtain

[TABLE]

where we have set $\mu=0$. A non–zero $\mu$ is trivially included via the replacement $E\cosh\eta\rightarrow E\cosh\eta-\mu$. Plugging this result into an analog of (22), we obtain $\Gamma_{2\rightarrow 2}$.

Finally, let us note that at high temperatures it is important to include the thermal mass term, i.e.

[TABLE]

This regularizes the behaviour of the rates at $T\gg m$ and cures the infrared divergence as $m\rightarrow 0$.

3.2 Conditions for thermal and kinetic equilibria

Thermal or kinetic equilibrium is maintained if the relevant reaction rate is larger than the Hubble rate. It is convenient to express this condition in terms of the quantities appearing in the Boltzmann equation. For full thermal equilibrium, we require

[TABLE]

which implies that the number changing interactions are efficient and the Bose–Einstein distribution is realized. On the other hand, kinetic equilibrium is maintained as long as

[TABLE]

such that the scattering rate is high enough to define a temperature. This is a weaker condition since normally $\Gamma_{2\rightarrow 2}\gg\Gamma_{2\rightarrow 4}$. The system is still described by the Bose–Einstein distribution which maximizes entropy, however a non–zero chemical potential is necessary to account for approximate particle number conservation.

Throughout this paper we assume that before and around freeze–out, the energy density of the Universe is dominated by the SM fields. This can be either due to higher temperature $T_{\rm SM}$ in the observed sector or due to a larger number of SM degrees of freedom $g_{*}$. In this case,

[TABLE]

In the relativistic regime, $T\gg m$, the temperature ratio

[TABLE]

remains constant. We find that the ratio $\Gamma_{2\rightarrow 4}/(nH)$ is maximized around $T\sim 5m/\sqrt{\lambda}$, i.e. when the thermal mass becomes comparable to the bare mass (see Sec. 4.3). Imposing condition (37) at this temperature, we obtain a lower bound on $\lambda$ at each $m$. The allowed parameter space is shown in Fig. 2. The values of $\lambda$ above the green lines are $consistent$ with thermal equilibrium for a given $\xi$.444The constraint provides the necessary condition for thermal equilibrium, while the thermalization process depends on further details such as the initial distribution function. (Naively, since $n\propto T^{3}$ and $\Gamma_{2\rightarrow 4}\propto T^{4}$, one expects the constraint (37) at high temperature to be of the type $m\ll T<{\rm const}\times\lambda^{4}\xi^{-2}M_{\rm Pl}$; while this shows the right trend, i.e. the minimal $\lambda$ increases with $m$, in reality the $\lambda$–dependence is more complicated due to the thermal mass contribution.)

Regarding kinetic equilibrium, the ratio $\Gamma_{2\rightarrow 2}/(nH)$ is maximized at $T\sim{\cal O}(m)$. Using our result for $\Gamma_{2\rightarrow 2}$, we find the allowed parameter space in Fig. 2.555The reaction rate is affected by the chemical potential. However, if $|\mu|$ is small compared to the temperature, the effect is not very significant. The couplings above the purple lines are necessary for maintaining kinetic equilibrium in the relativistic regime. Since the rate involves a lower power of $\lambda$, the $m$ and $\xi$ dependences are stronger than those in the case of full thermal equilibrium.

Here we require thermalization in the relativistic regime. In the non–relativistic case, the bound on the coupling is stronger (see Appendix B), so Fig. 2 gives the necessary condition.

We see that thermal equilibrium requires self–coupling of order $10^{-4}$–$10^{-3}$ at $m\sim 1$ GeV, while for kinetic equilibrium $\lambda$ can be as small as $10^{-8}$–$10^{-7}$. Thus, there is a large window in which only kinetic equilibrium is maintained and the DM number is approximately conserved.

4 Freeze–out

In what follows, we consider separately the non–relativistic and relativistic freeze–out regimes. Here we assume that the system enjoys full thermal equilibrium such that $\mu=0$ initially. In this case, the relic DM density is affected by the number changing interactions, which we study in detail.

4.1 Non–relativistic freeze–out

In the non–relativistic regime, the expressions for the reaction rates simplify. The momentum distribution is given by the Maxwell–Boltzmann function $f(p)=e^{-(E-\mu)/T}$ and the final state enhancement factors can be neglected. In this case, the chemical potential dependence factorizes out and according to (15) we have,

[TABLE]

It is conventional to $define$

[TABLE]

which is momentum independent in the non-relativistic limit $k_{i}\simeq(m,\vec{0})^{T}$. As usual, we absorb the symmetry factor $1/(2!4!)$ into $|{\cal M}_{4\rightarrow 2}|^{2}$. We thus have

[TABLE]

where the equilibrium particle density is $n_{\rm eq}=(2\pi)^{-3}\int d^{3}{\bf p}f(p)_{\mu=0}$ and $\langle...\rangle$ denotes a thermal average at $\mu=0$ over the momenta of the incoming particles. Expressing the chemical potential in terms of the particle densities as $e^{\mu/T}=n/n_{\rm eq}$, we obtain the Boltzmann equation in the form

[TABLE]

An important feature of this equation is that $\langle\sigma_{4\rightarrow 2}v^{3}\rangle$ is temperature independent.

It is convenient to replace the time variable with the SM sector temperature and the number density with the total particle number, which stays approximately constant when the number changing interactions become inefficient. Let us define

[TABLE]

where $s_{\rm SM}$ is the entropy density dominated by the SM contribution,

[TABLE]

with $g_{*s}$ being the number of degrees of freedom contributing to the entropy. This number is a function of the temperature such that

[TABLE]

In terms of the new variables, the Boltzmann equation reads

[TABLE]

where the modified Hubble rate is defined by

[TABLE]

and $H=\sqrt{\pi^{2}g_{*}/90}\;m^{2}/(x^{2}M_{\rm Pl})$. Observe that the $x$–dependent prefactor on the right hand side of the Boltzmann equation falls off sharply with $x$, namely as $x^{-8}$. This implies particle number conservation at late times.

Our next task is to derive $Y_{\rm eq}(x)$. Indeed, there are two unknowns in our system: $T(t)$ and $\mu(t)$ which should be determined by 2 equations. The second constraint comes from entropy conservation in the dark sector, $sa^{3}=const$, or

[TABLE]

which is constant in time. The entropy density in the non–relativistic limit is given by

[TABLE]

where in $s$ we also include the subleading term proportional to $T$. This allows us to express $T$ as

[TABLE]

Since $Y_{\rm eq}=Ye^{-\mu/T}$, we obtain

[TABLE]

This is to be inserted in the Boltzmann equation, while $c$ is determined by the boundary condition at $\mu=0$:

[TABLE]

and $Y_{0}$ is fixed by the initial dark and observed sector temperatures, $T_{0}$ and $T_{{\rm SM}0}$.

Now the Boltzmann equation can be solved numerically. We assume that at the initial point defined by $T_{0}$ and $T_{{\rm SM}0}$, the system enjoys thermal equilibrium ($\mu=0$). Then, $Y(x)$ is found by solving the Boltzmann equation with this boundary condition. In the non–relativistic limit, we find (see Appendix C)

[TABLE]

where we have factored out the $1/(2!4!)$ symmetry coefficient associated with the initial and final state phase space. The resulting solution for a representative set of input parameters is shown in Fig. 3.

The upper left panel of Fig. 3 shows that $Y(x)$ follows $Y_{\rm eq}$ closely up until $x\sim 25$, at which point it freezes–out and the number density remains approximately constant. As we see from the right upper panel, the freeze–out is well described by

[TABLE]

After that, $\Gamma_{4\rightarrow 2}$ becomes suppressed compared to $3nH$, whereas $\Gamma_{2\rightarrow 4}$ turns negligible even faster.

The effective chemical potential becomes appreciable, of order $T$, around the freeze–out point after which it can be approximated by a linear function of $T$ asymptotically approaching $m$, $m-\mu\propto T$. This follows from the entropy and particle number conservation in the dark sector (see Eq. 51). In this regime, $T\propto T_{\rm SM}^{2}$ as required by Eq. 52. Before the freeze–out, the $T$ dependence on $T_{\rm SM}$ is only logarithmic. Non–relativistic behaviour of DM leads to the heating of the dark sector [4], which can be viewed as a result of appreciable $\Gamma_{4\rightarrow 2}$ at that stage.

Note that $Y_{\rm eq}(x)$ is not a solution to the Boltzmann equation since $Y_{\rm eq}^{\prime}(x)$ does not vanish. The true solution $Y(x)$ is close to the “equilibrium” value until freeze–out. This implies that the right hand side of the Boltzmann equation is necessarily non–zero which entails reduction of the total DM number $na^{3}$ during this period. We find that the reaction rate difference $2(\Gamma_{4\rightarrow 2}-\Gamma_{2\rightarrow 4})$ is indeed not far from $3nH$ such that the number reduction is tangible.

After the freeze–out, the particle number is approximately constant, typically within 10%. Hence one can approximate

[TABLE]

where $x_{f}$ is the freeze–out point. This is a slightly different condition compared to what is often used in the literature. Here, we do not neglect the $Y_{\rm eq}(x)$ term which makes the number reduction less efficient, especially in the vicinity of the freeze–out point.

4.2 Relativistic freeze–out

For $m/T\lesssim 1$, the relativistic effects are important and one should use the Bose–Einstein distribution function. As Fig. 4 shows, the resulting reaction rates differ from their Maxwell–Boltzmann analogs by $10\%$ to $100\%$ at $m/T\sim 1$, while at $m/T\sim 0.1$ this difference can reach two orders of magnitude. The Bose–Einstein rates are greater due to the enhancement factors for the low energy states. Therefore, the effect is sensitive to the value of the thermal mass and $decreases$ for larger couplings.

As in the non–relativistic case, the evolution of the system is determined by the chemical potential and the temperature. The Boltzmann equation and entropy conservation fix $\mu$ and $T$ as functions of $T_{\rm SM}$. Defining

[TABLE]

and $x\equiv m/T_{\rm SM}$, we rewrite the Boltzmann equation as an equation for $\mu(y)$ and $x(y)$:666This can be derived by expressing the left hand side of the Boltzmann equation in terms of $dY/dx$ with $Y=n/s_{\rm SM}$, and calculating this derivative explicitly in terms of $f(p)$.

[TABLE]

where

[TABLE]

with $E(p)=\sqrt{m^{2}_{\rm eff}+{\bf p}^{2}}\equiv\sqrt{m^{2}+\lambda T^{2}/24+{\bf p}^{2}}$, $n=\int d^{3}{\bf p}/(2\pi)^{3}f(p)$ and $I_{1,2},\Gamma_{2\rightarrow 4}$, $\Gamma_{4\rightarrow 2}$, $\tilde{H},n$ are to be expressed in terms of $\mu,x$ and $y$. The modified Hubble rate $\tilde{H}$ is defined by Eq. 49. Note that we have included the thermal mass correction which leads to an extra contribution in Eq. 60.

The second equation for $\mu(y)$ and $x(y)$ is provided by the entropy conservation condition

[TABLE]

with

[TABLE]

and $s_{\rm SM}$ from Eq. 46. Here the energy density $\rho$ and the pressure $p$ are given by the standard Bose–Einstein formulas, and are to be expressed in terms of $\mu$ and $y$.

The two coupled equations can be solved numerically.777We compute $\Gamma_{2\rightarrow 4}$ by integrating numerically the CalcHEP output. We present our results in Fig.5. The corresponding freeze–out temperature is $T_{f}=1.2$ GeV with $m=1$ GeV, which makes the freeze–out regime relativistic.

Compared to the non–relativistic case, we observe a few differences. First, the evolution of $Y$ and $\mu$ is slower. Second, there is no “warming” period in which the dark temperature decreases much slower than $T_{\rm SM}$. After freeze–out, $T$ decreases faster than $T_{\rm SM}$ does. This is due to approximate conservation of the particle number $n/T_{\rm SM}^{3}$: the increase in $\mu$ gets compensated by a decrease in $T$. Eventually, $T\propto T_{\rm SM}^{2}$ in the non–relativistic regime. For comparison, we present the evolution of $\xi$ for relativistic and non–relativistic freeze–out in Fig. 6.

Again, we find that $Y(\infty)$ can be well approximated by $Y$ at freeze–out.

4.3 Ultra–relativistic freeze–out

Freeze–out at $T\gg m,m/\sqrt{\lambda}$ is not possible in our model. In this regime, $T/T_{\rm SM}$ stays constant by virtue of entropy conservation and

[TABLE]

Thus, if the thermal equilibrium condition $3nH<2\Gamma_{2\rightarrow 4}$ is satisfied at some point, it will continue to hold as long as dark matter remains ultra–relativistic.888This is in contrast to the SM neutrino case, where the reaction rate involves a higher power of $T$. As a result, no freeze out is possible.

We note that the thermal mass effect is important in this regime. Due to the infrared singularity, the high temperature rates contain the $T/m_{\rm eff}$ factors. For instance, $\Gamma_{2\rightarrow 2}\propto T^{4}\ln{T\over m_{\rm eff}}\;,$ where $m_{\rm eff}$ includes the thermal correction (36). At $T>5m/\sqrt{\lambda}$, the effective mass is dominated by the thermal term and the expected behaviour $\Gamma_{2\rightarrow 2}\propto T^{4}$ is reproduced. (Here we neglect the usual log–running of the coupling constant which is insignificant for $\lambda$ in the range of interest.) Similar considerations apply to $\Gamma_{2\rightarrow 4}$. The effect of the thermal mass is clearly seen in the left panel of Fig. 7: while in the non–relativistic regime the rate scales as $\lambda^{4}$, at higher temperatures this is no longer true. The $T^{4}$–behaviour is recovered when the thermal mass dominates.

Fig. 7 (right) collects the different regimes in the $\Gamma_{2\rightarrow 4}$ behavior and presents an overall picture. At high temperatures, $\Gamma_{2\rightarrow 4}/n\propto T$ evolves slower than $H\propto T_{\rm SM}^{2}$ does. This changes when the thermal mass becomes subdominant, which is marked as “rel./semi-rel.” in the plot. Finally, in the non–relativistic regime $\Gamma_{2\rightarrow 4}/n\propto n^{3}$ is exponentially suppressed. The magnitude of the rate relative to $H$ is determined by the coupling constant. The plot makes it clear that freeze–out in the ultra–relativistic regime is impossible. Also, if the dashed line is above the solid line, thermalization is never achieved. This is determined by $\Gamma_{2\rightarrow 4}$ in the relativistic/semi–relativistic regime where $T$ is not too far from $m$. The resulting lower bounds on $\lambda$ are presented in Fig. 2.

5 Parameter space analysis

In this section, we delineate our parameter space and determine regions consistent with the dark matter relic abundance constraint as well as other relevant bounds. We consider separately the full thermal and kinetic equilibrium cases.

5.1 Thermalized dark matter

In this subsection, we assume that dark matter has been in thermal equilibrium and analyze the $(m,\lambda,\xi)$–parameter space consistent with the observed DM relic density. This constraint can be put in the form

[TABLE]

Let us discuss the main qualitative features of the model. Consider first the non–relativistic freeze–out regime. As discussed earlier, we define the freeze–out point by

[TABLE]

Solving this equation for $T_{f}$ and equating $Y(x_{f})=Y(\infty)$, we find that the correct relic density is reproduced along the curves with the approximate scaling

[TABLE]

where $\xi_{f}$ is $T_{\rm SM}/T$ at the freeze–out point and we have neglected the logarithmic terms. For a fixed $\xi_{f}$, the mass–coupling relation is approximately linear and the freeze out temperature decreases with the coupling,

[TABLE]

where the constant is positive. Thus, at sufficiently small $\lambda$, the non-relativistic approximation breaks down and a fully relativistic analysis is necessary.

Our numerical results with the relativistic reaction rates are presented in Fig. 8. We observe that in most allowed parameter space the freeze–out is relativistic, i.e. $m/T_{f}<3$ as marked by the green dotted line. In this region, the constant $\xi_{f}$ curves tend to approach a vertical at $T\gtrsim m$ since $Y$ is determined mostly by $\xi$ and rather insensitive to $\lambda$. These curves exhibit characteristic kinks due to a change in the SM degrees of freedom at QCD phase transition ($T_{\rm SM}\sim 10^{-1}$ GeV) and electron decoupling. It is interesting that $T_{f}$ does not vary monotonically along the constant $\xi_{f}$ curve: first, it decreases down to a minimum value and increases from there on.

For a given $\xi_{f}$, the correct relic density curve cannot be continued indefinitely to smaller $\lambda$: at some point, $3nH>2\Gamma_{2\rightarrow 4}$ for any $T$ (see Sec. 4.3). Thus, the relic density and freeze–out equations have no simultaneous solution. The excluded region is marked “no thermalization” in Fig. 8. Although different points on the border of this region correspond to different $\xi$, its shape is consistent with Fig. 2. We note that here the thermal mass effect is significant which makes analytical calculations more challenging.

Another constraint is imposed by the bound on DM self–interaction from the Bullet Cluster, $\sigma/m<1\;{\rm cm^{2}/g}$, with $\sigma=\lambda^{2}/(128\pi m^{2})$. It excludes light DM with significant self–coupling. Finally, perturbative unitarity is violated in the process $SS\rightarrow SS$ if $\lambda\gtrsim 8\pi$ [27]. We note that the nucleosynthesis constraint on the effective number of neutrinos is insignificant here since the dark matter contribution to the energy density in the relativistic regime is suppressed by $T^{4}/T^{4}_{\rm SM}\ll 1$.

In the non–relativistic regime, we find qualitative agreement with the results of Ref. [13] although there are numerical differences.

Fig. 8 shows that the correct relic density can be obtained for a wide range of DM masses: from 10 keV to 100 TeV, as long as the dark temperature is significantly below $T_{\rm SM}$. Thus, both the “warm” and the “cold” options are open.

5.2 Dark matter in kinetic equilibrium and $\mu\not=0$

For small enough self–coupling (see Fig. 2), the system reaches only kinetic but not thermal equilibrium. The notion of temperature is still well defined in this case, but the Bose–Einstein distribution involves an effective chemical potential from the start. The latter is determined by the initial number density. Assuming $T\gg m$, we have

[TABLE]

where ${\rm Li}_{3}(x)$ is the third degree polylogarithm, ${\rm Li}_{3}(x)=\sum_{n=1}^{\infty}x^{n}/n^{3}$. For a given initial $n$, $\mu$ is read off from (69).

We are mostly interested in $\mu<0$ which suppresses the dark matter density. In theories with antiparticles, the antiparticle distribution involves $-\mu$ which restricts $|\mu|<m$ [28, 29]. This does not apply to the model at hand and the $S$–gas can be arbitrarily dilute as long as kinetic equilibrium is maintained. At large $-\mu$, the density is exponentially suppressed, $n\propto e^{\mu/T}$.

It is straightforward to derive the temperature and chemical potential evolution. Since both the total number and entropy are conserved, away from the thresholds where the number of SM degrees of freedom changes, we have

[TABLE]

such that at $T\gg m$ one finds

[TABLE]

Our numerical results are shown in Fig. 9. It displays the constant $\xi$ curves producing the correct DM relic density in the $(m,-\mu/T)$ plane.999We focus on $\mu<0$ since positive values of $\mu$ are bounded by $m$ and omit the possibility of the Bose–Einstein condensate formation. Here the parameters $\xi$ and $\mu/T$ are defined in the relativistic regime ($T=10m$). The fixed relic density lines correspond approximately to $e^{-\mu/T}\propto m/\xi^{3}$. For a fixed $\xi$ and small $m$, the dilution factor $e^{-\mu/T}$ becomes insignificant and the relativistic dark sector starts making a substantial contribution to the energy density of the Universe. Throughout this paper we assume that the SM thermal bath dominates the energy density balance, hence we exclude the region $\rho_{\rm DM}/\rho_{\rm SM}>0.05$, which is also disfavored by cosmological bounds on the effective number of neutrinos [30].

The right panel of Fig. 9 shows that significantly larger couplings are required for kinetic equilibrium as the gas becomes more dilute. The rate $\Gamma_{2\rightarrow 4}$ is still smaller than $\Gamma_{2\rightarrow 2}$, so full thermal equilibrium is not reached for a significant range of the couplings. We find that the Bullet Cluster bound on the self–coupling is insignificant here and superseded by the constraint $\rho_{\rm DM}/\rho_{\rm SM}<0.05$.

Note that the dark sector in kinetic equilibrium is allowed to be significantly hotter than the observable one, $T\gg T_{\rm SM}$. This is consistent with the constraints due to the exponential suppression of $\rho_{\rm DM}$ for substantial $-\mu/T$.

Overall, we find that there is a large range of DM masses consistent with observations. For a given $m$ within this range, there exists an initial DM density which leads to the correct relic abundance. Thermodynamical considerations apply to this system as long as the self–coupling is large enough to bring it into kinetic equilibrium. For smaller couplings, the correct relic abundance can still be obtained, however the dark temperature is ill–defined.

6 Conclusion

We have performed a comprehensive study of real scalar dark matter decoupled from the Standard Model fields. We have considered two regimes where the dark sector can be assigned a temperature:

• •

DM in thermal equilibrium (larger self–coupling)

• •

DM in kinetic equilibrium (smaller self–coupling)

In the latter case, the relic abundance is fixed by the initial number density which corresponds to a non–zero effective chemical potential. In the former case, it is determined by the freeze–out temperature below which the number changing interactions are suppressed.

We have developed a relativistic approach to the dark matter evolution. In particular, we use fully relativistic expressions for the number changing and number conserving reaction rates. This allows us to explore the relativistic freeze–out regime, which occurs commonly in the allowed parameter space. When the dark temperature is much smaller than the observable one, the correct DM relic abundance can be obtained for a wide range of DM masses, 10 keV to 100 TeV. The required self–coupling is above $10^{-5}$.

If dark matter reaches only kinetic equilibrium, the correct relic density can be obtained both for $T<T_{\rm SM}$ and $T>T_{\rm SM}$ as long as DM is sufficiently dilute, $-\mu\gtrsim T$. The allowed DM mass is then above 100 eV. The presence of chemical potential also suppresses the effect of relativistic DM on nucleosynthesis.

Altogether, there is vast parameter space consistent with the thermal history of the Universe, while dark matter can be warm or cold.

Acknowledgements

O.L. is indebted to Aleksi Vuorinen for numerous enlightening discussions, to Keijo Kajantie, Kirill Boguslavski and Howie Haber for useful communications, and to Eliza Dickie for verifying some of the analytical results. G.A. thanks Stefan Vogl for fruitful discussions. The work of S.P. is partially supported by the National Science Centre, Poland, under research grants DEC-2015/18/M/ST2/00054 and DEC-2016/23/G/ST2/04301. T.T. acknowledges funding from the Natural Sciences and Engineering Research Council of Canada (NSERC). This research was enabled in part by support provided by Compute Ontario, WestGrid, Compute Canada, the T30/CIP-cluster at Technical University of Munich and the Yukawa Institute Computer Facility.

Appendix A

Here we provide some of the explicit expressions for the conversion from a general reference frame to the center–of–mass frame. Given 2 momenta $p_{1},p_{2}$, we define $p=(p_{1}+p_{2})/2$ and $k=(p_{1}-p_{2})/2$. The center–of–mass frame is defined by the relation

[TABLE]

where $E$ is the particle energy in the center–of–mass frame. In the convention $p=(p^{0},p^{3},p^{2},p^{1})^{T}$, the explicit form of $\Lambda$ and its inverse in terms of the rapidity $\eta$ and angular variables $\theta,\phi$ is given by

[TABLE]

In this frame, the 6 degrees of freedom of $p_{1},p_{2}$ become $E,\eta$ and 4 angles $\theta,\phi,\theta_{k},\phi_{k}$, where $\theta_{k}$ and $\phi_{k}$ are the spherical coordinate angles parameterizing $k$. Note that $k^{0}=0$ and $|{\bf{k}}|=\sqrt{E^{2}-m^{2}}$ due to the on–shell condition for the initial particles.

As a result, for the rest frame velocity $u=(1,0,0,0)^{T}$ and final state momenta $k_{i}$ in the center–of–mass frame, we have

[TABLE]

which appears in the Bose–Einstein enhancement factors for the final state. In the initial state thermal averaging, we encounter

[TABLE]

where we have used $k^{3}=|{\bf k}|\cos\theta_{k}$.

Appendix B

Here we present a non-relativistic analog of Fig. 2. The reaction rate at $\mu\sim 0$ is

[TABLE]

with a non–relativistic number density $n_{\rm eq}$ and $\langle\sigma_{4\rightarrow 2}v^{3}\rangle$ given by (55). We require $T<m/5$ for DM to be sufficiently non–relativistic. The condition $3nH<2\Gamma_{42}$ translates into the lower bound on $\lambda$ shown in Fig. 10. The result depends on how $T_{\rm SM}/T$ is fixed: the left panel shows the bounds for $T_{\rm SM}/T=1..100$ at $T=m/5$, whereas the right panel displays the bounds for $T_{\rm SM}/T=1..100$ fixed in the relativistic regime, $T=5m$, and continued to lower $T<m/5$ using entropy conservation,

[TABLE]

The reaction rates involve at least $n^{2}$ and thus drop sharply with the temperature, as illustrated in Fig. 3. Therefore, the equilibrium condition for a given $(\lambda,m)$ should be tested at the highest $T$ consistent with the non–relativistic regime. If it is not satisfied at this point, dark matter will not thermalize at lower $T$ either. Similar considerations apply to the kinetic equilibrium condition.

The resulting bounds on the coupling are stronger than those in the relativistic case. In particular, comparison of the right panel of Fig. 10 with Fig. 2 shows that if DM is relativistic initially and does not thermalize in the relativistic regime, it will not reach thermal equilibrium later either.

Appendix C

In this appendix, we provide details of the $4\rightarrow 2$ cross section calculation in the non–relativistic limit. We believe this is useful since there are significant discrepancies with the results in the literature [13, 16].

There are two distinct contributions to the amplitude shown in Fig. 1. The initial state momenta are all $(m,\vec{0})^{T}$, while the final state momenta can be chosen as $(2m,\sqrt{3}m,0,0)^{T}$ and $(2m,-\sqrt{3}m,0,0)^{T}$. Consider the diagram on the right. The momentum flow through the propagator is $p^{2}=-3m^{2}$, while the symmetry factor is $1/(2!)\times 12^{2}\times 2!4!$, where $1/(2!)$ is due to the 2d order in the coupling and $2!4!$ is due to the permutations among the initial and final state legs. For the left diagram, the momentum flow through the propagator is $p^{2}=9m^{2}$ and the symmetry factor is $1/(2!)\times 96\times 2!4!$ Thus, the QFT amplitude defined in the $standard$ way is given by

[TABLE]

In our convention, we include the phase space symmetry factor for the initial and final state, $1/(2!4!)$, directly in the cross section. Thus, according to (42) we have

[TABLE]

which reproduces (55).

This result can be verified numerically. The $2\rightarrow 4$ cross section is calculated with CalcHEP, while in equilibrium the $2\rightarrow 4$ and $4\rightarrow 2$ rates are related through $\langle\sigma v^{3}\rangle=\langle\sigma v\rangle/n_{\rm eq}^{2}$, according to our cross section convention. We find excellent agreement with our analytical formula.