Boltzmann Equations

General Formalism and Approximations

The general Boltzmann equation for the number distribution of a particle species can be written as (assuming isotropy):

{\frac{\partial F_{i}}{\partial t}} -H p {\frac{\partial F_{i}}{\partial p}} = C_{i}[F_{i},F_{j},p]

where \(F_{i}(p)\) is the number distribution of particle \(i\) as function of momentum \(p\), \(C\) represents a source/sink term and \(H\) is the Hubble constant:

\[H = \sqrt{\frac{8 \pi}{3} \frac{\rho_T}{M_P^2}} \label{H}\]

with \(M_{P} = 1.22\times 10^{19}\) GeV and \(\rho_T = \sum_{i} \rho_i\). The number, energy and pressure densities are given in terms of \(F_{i}\) as:

\[\begin{split}\begin{align} n_{i}(t) & = \int \frac{dp}{2 \pi^2} p^2 F_i(p) \nonumber \\ \rho_{i}(t) & = \int \frac{dp}{2 \pi^2} p^2 E_i F_i(p) \label{beqs}\\ P_{i}(t) & = \frac{1}{3} \int \frac{dp}{2 \pi^2} \frac{p^4}{E_i} F_i(p) \nonumber\end{align}\end{split}\]

where \(m_i\) is the mass of particle \(i\) and \(E_i = \sqrt{p_i^2 + m_i^2}\). Using the Boltzmann equation above, we obtain the following equations for the number and energy densities:

\[\begin{split}\begin{aligned} {\frac{d n_i}{d t}} + 3H n_i & = \int \frac{dp}{2 \pi^2} p^2 C_i \nonumber \\ {\frac{d \rho_i}{d t}} + 3H (\rho_i + P_i) & = \int \frac{dp}{2 \pi^2} p^2 E_i C_i \label{eq:meqs}\end{aligned}\end{split}\]

The collision term, \(C_i\), for the process \(i + f + \ldots \leftrightarrow a + b + c + \ldots\) is given by:

\[\begin{split}\begin{aligned} C_i & = \frac{1}{E_i} \int \prod_{j,a} \frac{d^3 p_j}{2 E_j (2 \pi)^3} \frac{d^3 p_a}{2 E_a (2 \pi)^3} (2 \pi)^4 \delta^{4}\left(p_i + p_j + \ldots - p_a - p_b \ldots\right) |\mathcal{M}|^2 \nonumber \\ &\times \left[(1 \pm f_a) (1 \pm f_b)\ldots f_i f_j\ldots - f_a f_b \ldots (1 \pm f_i)(1 \pm f_j)\ldots \right]\end{aligned}\end{split}\]

where the plus (minus) sign is for bosons (fermions). Below we always assume \(f_{i,j,a,..} \ll 1\), so:

\[\begin{split}\begin{aligned} C_i & \simeq \frac{1}{E_i} \int \prod_{j,a} \frac{d^3 p_j}{2 E_j (2 \pi)^3} \frac{d^3 p_a}{2 E_a (2 \pi)^3} (2 \pi)^4 \delta^{4}\left(p_i + p_j + \ldots - p_a - p_b \ldots\right) |\mathcal{M}|^2 \nonumber \\ &\times \left[f_i f_j\ldots - f_a f_b \ldots \right]\end{aligned}\end{split}\]

We will assume that \(C\) is given by:

\[C = C_{dec} + C_{prod} + C_{ann}\]

where \(C_{dec}\) contains the contributions from decays and inverse decays (\(i \leftrightarrow a + b + \ldots\)), \(C_{prod}\) contains the contributions from decay injection and inverse decay injection (\(a \leftrightarrow i + b + \ldots\)) and \(C_{ann}\) from annihilations with the thermal plasma (\(i + i \leftrightarrow a + b\)). Below we compute each term separately, under some assumptions.

Annihilation Term

The annihilation term \(C_{ann}\) for the \(i + j \leftrightarrow a + b\) process is given by:

\[\int \frac{dp}{2 \pi^2} p^2 C_{ann} = \int d\Pi_{i} d\Pi_{j} d\Pi_{a} d\Pi_{b} (2 \pi)^4 \delta^{(4)}(p_i + p_j - p_a - p_b) |M|^2 \left[ f_a f_b - f_i f_j \right]\]

where \(d\Pi_{i} = d^{3} p_i/((2\pi)^3 2 E_i)\). Since we are ultimately interested in Eqs.([eq:meqs]) for the number and energy densities, we will consider the following integral:

\[\begin{split}\int \frac{dp}{2 \pi^2} p^2 C_{ann} E_i^{\alpha} = \int d\Pi_{i} d\Pi_{j} d\Pi_{a} d\Pi_{b} (2 \pi)^4 \delta^{(4)}(p_i + p_j - p_a - p_b) |M|^2 \\ \times \left[ f_a f_b - f_i f_j \right] E_i^{\alpha}\end{split}\]

where \(\alpha = 0 (1)\) for the number (energy) density. Here we assume that the distributions can be approximated by:

\[f_i \simeq \exp(-(E_i - \mu_i)/T)\]

so the annihilation term can then be written as:

\[\begin{split}\begin{aligned} & \int \frac{dp}{2 \pi^2} p^2 C_{ann} E_i^{\alpha} = -\left( \exp((\mu_i + \mu_j)/T) -\exp((\mu_a + \mu_b)/T)\right) \nonumber \\ & \times \int d\Pi_{i} d\Pi_{j} d\Pi_{a} d\Pi_{b} (2 \pi)^4 \delta^{(4)}(p_i + p_j - p_a - p_b) |M|^2 \exp(-(E_i + E_j)/T) \times E_i^{\alpha} \nonumber\end{aligned}\end{split}\]

where above we have used conservation of energy (\(E_i + E_j = E_a + E_b\)). Since for the cases of interest the equilibrium distributions have zero chemical potential, we have:

\[\frac{n_i}{\bar{n}_i} = \exp(\mu_i/T)\]

so:

\[\begin{split}\begin{aligned} & \int \frac{dp}{2 \pi^2} p^2 C_{ann} E_i^{\alpha} = -\left( \frac{n_i n_j}{\bar{n}_i \bar{n}_j} - \frac{n_a n_b}{\bar{n}_a \bar{n}_b}\right) \nonumber \\ & \times \int d\Pi_{i} d\Pi_{j} d\Pi_{a} d\Pi_{b} (2 \pi)^4 \delta^{(4)}(p_i + p_j - p_a - p_b) |M|^2 \exp(-(E_i + E_j)/T) \times E_i^{\alpha} \nonumber\end{aligned}\end{split}\]

In particular, for the process \(i + i \leftrightarrow a + b\), where \(a\) and \(b\) are in thermal equilibrium (\(\mu_a = \mu_b = 0\)):

\[\begin{split}\begin{aligned} & \int \frac{dp}{2 \pi^2} p^2 C_{ann} E_i^{\alpha} = -\left( \frac{n_i^2}{\bar{n}_i^2} - 1 \right) \nonumber \\ & \times \int d\Pi_{i} d\Pi_{j} d\Pi_{a} d\Pi_{b} (2 \pi)^4 \delta^{(4)}(p_i + p_j - p_a - p_b) |M|^2 \exp(-(E_i + E_j)/T) \times E_i^{\alpha} \nonumber \\ & = -\left( n_i^2 - \bar{n}_i^2 \right) \langle \sigma v E_i^{\alpha} \rangle\end{aligned}\end{split}\]

For \(\alpha = 0\), the above equation is the well known contribution from thermal scatterings to the annihilation term. To estimate its value for \(\alpha = 1\), we assume:

\[\langle \sigma v E \rangle \simeq \langle \sigma v \rangle \langle E_i \rangle = \langle \sigma v \rangle \frac{\rho_i}{n_i} \label{eq:app}\]

where \(\langle \;\; \rangle\) represents thermal average. Thus:

\[\begin{split}\int \frac{dp}{2 \pi^2} p^2 C_{ann} E_i^{\alpha} = \left( \bar{n}_i^2 - n_i^2 \right) \left\{ \begin{array}{rl} \langle \sigma v \rangle & \mbox{, for $\alpha = 0$} \\ \langle \sigma v \rangle \frac{\rho_i}{n_i} &\mbox{, for $\alpha = 1$} \end{array} \right. \label{eq:collfin}\end{split}\]

Decay Term

Now we derive a simplified expression for the decay (and inverse decay) term, under approximations similar to the ones used in the last section. The decay term includes the contributions from particle decay and inverse decay:

\[C_{dec} \simeq \frac{1}{E_i} \int \prod_{a} \frac{d^3 p_a}{2 E_a (2 \pi)^3} (2 \pi)^4 \delta^{4}\left(p_i - p_a - p_b \ldots\right) |\mathcal{M}|^2 \left[f_i - f_a f_b \ldots \right] \label{eq:dec0}\]

As in the case of the annihilation term, we assume that the distributions for \(a,b,\ldots\) can be approximated by \(f_x \simeq \exp(-(E_x - \mu_x)/T)\), so we can write:

\[\begin{split}f_a f_b \ldots \simeq \exp\left(\frac{\mu_a + \mu_b + \ldots}{T}\right) \exp(-E_i/T) = \frac{n_a n_b \ldots}{\bar{n}_a \bar{n}_b \ldots} \exp(-E_i/T) \\ = \frac{n_a n_b \ldots}{\bar{n}_a \bar{n}_b \ldots} \bar{f}_{i}\end{split}\]

where we used conservation of energy (\(E_a + E_b + \ldots = E_i\)) and \(\bar{f}_i\) is the equilibrium distribution for the species \(i\). Hence we can write Eq.([eq:dec0]) as:

\[\begin{split}\begin{aligned} C_{dec} & \simeq \left[f_i - \frac{n_a n_b \ldots}{\bar{n}_a \bar{n}_b \ldots} \bar{f}_{i} \right] \frac{1}{E_i} \int \prod_{a} \frac{d^3 p_a}{2 E_a (2 \pi)^3} (2 \pi)^4 \delta^{4}\left(p_i - p_a - p_b \ldots\right) |\mathcal{M}|^2 \nonumber \\ & = \mathcal{B}_{ab\ldots} \frac{\Gamma_i m_i}{E_i} \left[f_i - \frac{n_a n_b \ldots}{\bar{n}_a \bar{n}_b \ldots} \bar{f}_{i} \right] \end{aligned}\end{split}\]

where \(\Gamma_i\) is the width for \(i\) and \(\mathcal{B}_{ab\ldots} \equiv BR(i \to a + b + \ldots)\)

Once again we consider the integral:

\[\begin{split}\begin{aligned} \int \frac{dp}{2 \pi^2} p^2 C_{dec}(p) E_i^{\alpha} = & - \Gamma_i \int \frac{dp}{2 \pi^2} p^2 \frac{m_i}{E_i} f_i E_i^{\alpha} \nonumber \\ & + \sum_{i \; decays} \mathcal{B}_{ab\ldots} \Gamma_i \frac{n_a n_b \ldots}{\bar{n}_a \bar{n}_b \ldots} \int \frac{dp}{2 \pi^2} p^2 \frac{m_i}{E_i} \bar{f}_{i} E_i^\alpha \label{eq:dec2}\end{aligned}\end{split}\]

where we have included the sum over all decay channels and \(\alpha = 0 (1)\) for the contribution to the number (energy) density equation. Note that both integrals are identical, except for the replacement \(f_i \to \bar{f_i}\). The first integral in Eq.([eq:dec2]) gives:

\[\begin{split}-\Gamma_i \int \frac{dp}{2 \pi^2} p^2 \frac{m_i}{E_i} f_i(p) E_i^{\alpha} = \left\{ \begin{array}{rl} -\Gamma_i m_i n_i \langle \frac{1}{E_i} \rangle & \mbox{, for $\alpha = 0$} \\ -\Gamma_i m_i n_i &\mbox{, for $\alpha = 1$} \end{array} \right. \label{eq:dec1a}\end{split}\]

where

\[\langle \frac{1}{E_i} \rangle \equiv \frac{1}{n_i} \int \frac{dp}{2 \pi^2} p^2 \frac{1}{E_i} f_i(p)\]

Hence we can write Eq.([eq:dec2]) as:

\[\begin{split}\int \frac{dp}{2 \pi^2} p^2 C_{dec}(p) E_i^{\alpha} = -\Gamma_i m_i \left\{ \begin{array}{ll} n_i \langle \frac{1}{E_i} \rangle - \bar{n}_i \langle \frac{1}{E_i} \rangle_{eq} \sum \mathcal{B}_{ab\ldots} \frac{n_a n_b\ldots}{\bar{n}_a \bar{n}_b\ldots} & \mbox{, for $\alpha = 0$} \\ n_i - \bar{n}_i \sum \mathcal{B}_{ab\ldots} \frac{n_a n_b\ldots}{\bar{n}_a \bar{n}_b\ldots} & \mbox{, for $\alpha = 1$} \end{array} \right. \label{eq:decfin}\end{split}\]

For the non-equilibrium average we assume:

\[\langle \frac{1}{E_i} \rangle \simeq \frac{1}{\langle E_i \rangle} = \frac{n_i}{\rho_i}\]

which is exact in the non-relativistic limit, but it is only an approximation for the relativistic case. Although we can compute the equilibrium average (\(\langle \frac{1}{E_i}\rangle_{eq}\)) explicitly, in order to have an exact cancellation between the decay and inverse decay terms when \(i\), \(a\) and \(b\) are all in equilibrium, we take:

\[\langle \frac{1}{E_i} \rangle_{eq} \simeq \langle \frac{1}{E_i} \rangle = \frac{n_i}{\rho_i}\]

With the above approximations we finally obtain:

\[\begin{split}\int \frac{dp}{2 \pi^2} p^2 C_{dec}(p) E_i^{\alpha} = -\Gamma_i m_i \left\{ \begin{array}{ll} \frac{n_i}{\rho_i}\left( n_i - \bar{n}_i \sum \mathcal{B}_{ab\ldots} \frac{n_a n_b \ldots}{\bar{n}_a \bar{n}_b \ldots} \right) & \mbox{, for $\alpha = 0$} \\ n_i - \bar{n}_i \sum \mathcal{B}_{ab\ldots} \frac{n_a n_b \ldots}{\bar{n}_a \bar{n}_b \ldots} & \mbox{, for $\alpha = 1$} \end{array} \right.\end{split}\]

where \(\mathcal{B}_{ab\ldots} \equiv BR(i\to a+b+\ldots)\).

Production Term

The decay and inverse decay of other particles (\(a \to i + b + \ldots\)) can also affect the species \(i\). The contribution from these terms we label \(C_{prod}\), which is given by:

\[\begin{split}C_{prod} \simeq \frac{1}{E_i} \int \frac{d^3 p_a}{2 E_a (2 \pi)^3} \prod_{b} \frac{d^3 p_b}{2 E_b (2 \pi)^3} (2 \pi)^4 \delta^{4}\left(p_a - p_i - p_b \ldots\right) |\mathcal{M}|^2\\ \times \left[f_a - f_i f_b \ldots \right]\end{split}\]

Using the same approximations of the previous section, we write:

\[f_i f_b\ldots \simeq \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} e^{-E_a/T} = \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \bar{f}_{a}\]

Hence:

\[\begin{split}C_{prod} = \frac{1}{E_i} \int \frac{d^3 p_a}{2 E_a (2 \pi)^3} \prod_{b} \frac{d^3 p_b}{2 E_b (2 \pi)^3} (2 \pi)^4 \delta^{4}\left(p_a - p_i - p_b \ldots\right) |\mathcal{M}|^2 \\ \times \left(f_a - \bar{f}_a \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \right)\end{split}\]

and

\[\begin{split}\begin{aligned} \int \frac{dp}{2 \pi^2} p^2 C_{prod}(p) E_i^\alpha & = \int \frac{d^3 p_a}{E_a (2 \pi)^3} \left(f_a - \bar{f}_a \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \right) \nonumber \\ & \times \frac{d^3 p E_i^{\alpha}}{2 E_i (2 \pi)^3} \prod_{b} \frac{d^3 p_b}{2 E_b (2 \pi)^3} (2 \pi)^4 \delta^{4}\left(p_a - p_i - p_b \ldots\right) |\mathcal{M}|^2 \label{eq:prod2}\end{aligned}\end{split}\]

with \(\alpha = 0 (1)\) for the contribution to the number (energy) density equation. For \(\alpha = 0\) we obtain:

\[\begin{split}\begin{aligned} \int \frac{dp}{2 \pi^2} p^2 C_{prod}(p) & = \Gamma_a \mathcal{B}_{i} m_a \int \frac{d^3 p_a}{E_a (2 \pi)^3} \left(f_a - \bar{f}_a \sum_b \frac{\mathcal{B}_{ib\ldots}}{\mathcal{B}_{i}}\frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \right) \nonumber \\ & = \Gamma_a \mathcal{B}_{i} m_a \frac{n_a}{\rho_a} \left( n_a - \bar{n}_a \sum_b \frac{\mathcal{B}_{ib\ldots}}{\mathcal{B}_{i}} \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \right)\end{aligned}\end{split}\]

where \(\mathcal{B}_{ib\ldots} \equiv BR(a \to i + b + \ldots)\), \(\mathcal{B}_i = \sum_{b} \mathcal{B}_{ib\ldots}\) and we have once again assumed \(\langle 1/E_a \rangle \simeq \langle 1/E_a \rangle_{eq} \simeq n_a/\rho_a\).

For \(\alpha = 1\), the integral in Eq.([eq:prod2]) does not take a simple form. In order to compute it, we assume:

\[E_i \simeq \frac{E_a}{2}\]

The above expression is only exact for 2-body decays and \(m_a \gg m_i,m_b\). For the remaining cases, it is only an estimate.

\[\begin{split}\begin{aligned} \int \frac{dp}{2 \pi^2} p^2 C_{prod}(p) E_i & \simeq \Gamma_a \mathcal{B}_{i} \frac{m_a}{2} \int \frac{d^3 p_a}{(2 \pi)^3} \left(f_a - \bar{f}_a \sum_b \frac{\mathcal{B}_{ib\ldots}}{\mathcal{B}_{i}} \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \right) \nonumber \\ & = \Gamma_a \mathcal{B}_{i} \frac{m_a}{2} \left( n_a - \bar{n}_a \sum_b \frac{\mathcal{B}_{ib\ldots}}{\mathcal{B}_{i}} \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \right)\end{aligned}\end{split}\]

Combining the results for \(\alpha = 0\) and 1, we have:

\[\begin{split}\int \frac{dp}{2 \pi^2} p^2 C_{prod}(p) E_i^{\alpha} = \Gamma_a \mathcal{B}_{i} m_a \left( n_a - \bar{n}_a \sum_b \frac{\mathcal{B}_{ib\ldots}}{\mathcal{B}_{i}} \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \right) \left\{ \begin{array}{ll} \frac{n_a}{\rho_a} & \mbox{, for $\alpha = 0$} \\ \frac{1}{2} & \mbox{, for $\alpha = 1$} \end{array} \right. \label{eq:prodfin}\end{split}\]

Number and Energy Density Equations

Using the results of Eqs.([eq:collfin]), ([eq:decfin]) and ([eq:prodfin]) in the Boltzmann equations for \(n_i\) and \(\rho_i\) (Eq.([eq:meqs])), we obtain:

\[\begin{split}\begin{aligned} {\frac{d n_i}{d t}} + 3H n_i & = \left( \bar{n}_i^2 - n_i^2 \right) \langle \sigma v \rangle - \Gamma_i m_i \frac{n_i}{\rho_i}\left(n_i - \bar{n}_i \sum_{i\to\ldots} \mathcal{B}_{ab\ldots} \frac{n_a n_b \ldots}{\bar{n}_a \bar{n}_b \ldots} \right) \nonumber \\ & + \sum_a \Gamma_a \mathcal{B}_i m_a \frac{n_a}{\rho_a} \left(n_a - \bar{n}_a \sum_{a \to i\ldots} \frac{\mathcal{B}_{ib\ldots}}{\mathcal{B}_{i}} \frac{n_i n_b \ldots}{\bar{n}_i \bar{n}_b \ldots} \right) + C_{i}(T) \label{eq:nieq} \\ {\frac{d \rho_i}{d t}} + 3H (\rho_i + P_i) & = \left( \bar{n}_i^2 - n_i^2 \right) \langle \sigma v \rangle \frac{\rho_i}{n_i} - \Gamma_i m_i \left( n_i - \bar{n}_i \sum_{i\to\ldots} \mathcal{B}_{ab\ldots} \frac{n_a n_b\ldots}{\bar{n}_a \bar{n}_b\ldots}\right) \nonumber \\ & + \sum_a \Gamma_a \mathcal{B}_i \frac{m_a}{2} \left( n_a - \bar{n}_a \sum_{a \to i\ldots} \frac{\mathcal{B}_{ib\ldots}}{\mathcal{B}_{i}} \frac{n_i n_b..}{\bar{n}_i \bar{n}_b..} \right) + \tilde{C}_{i}(T) \frac{\rho_i}{n_i}\end{aligned}\end{split}\]

where \(\mathcal{B}_{ab\ldots} = BR(i \to a + b+ \ldots)\), \(\mathcal{B}_{ib\ldots} = BR(a \to i + b + \ldots)\), \(\mathcal{B}_i = \sum_b \mathcal{B}_{ib\ldots}\) and we have included an extra term (\(C_i\) and \(\tilde{C}_i\)) to allow for other possible sources for the number and energy densities. For simplicity we assume \(C_i = \tilde{C}_{i}\) from now on.

It is also convenient to use the above results to obtain a simpler equation for \(\rho_i/n_i\):

\[\begin{split}\begin{aligned} {\frac{d \rho_i/n_i}{d t}} \equiv {\frac{d R_i}{d t}} & = -3 H \frac{P_i}{n_i} \\ & + \sum_{a} \mathcal{B}_{i} \frac{\Gamma_a m_a}{n_i} \left( \frac{1}{2} - \frac{n_a}{\rho_a} \frac{\rho_i}{n_i} \right) \left(n_a - \bar{n}_a \sum_{a \to i\ldots} \frac{\mathcal{B}_{ib\ldots}}{\mathcal{B}_{i}} \frac{n_i n_b..}{\bar{n}_i \bar{n}_b..}\right) \label{eq:Rieq} \end{aligned}\end{split}\]

Besides the above equations, it is useful to consider the evolution equation for entropy:

\[dS \equiv \frac{dQ^{dec}}{T}\]

where \(dQ^{dec}\) is the net energy injected from decays. With the above definition we have:

\[\begin{split}\begin{aligned} \dot{S} & = \frac{1}{T}\sum_i BR(i,X) \frac{d\left(R^3 \rho_i\right)^{dec}}{dt} \nonumber \\ \Rightarrow \dot{S} & = \frac{R^3}{T}\sum_i BR(i,X) \Gamma_i m_i\left(n_i - \bar{n}_i \sum_{i\to\ldots} \mathcal{B}_{ab\ldots} \frac{n_a n_b\ldots}{\bar{n}_a \bar{n}_b\ldots} \right) \label{Seq}\end{aligned}\end{split}\]

where \(R\) is the scale factor and \(BR(i,X)\) is the fraction of energy injected in the thermal bath from \(i\) decays.

The above expressions can be written in a more compact form if we define the following ”effective thermal densities” and ”effective BR”:

\[\begin{split}\begin{aligned} \mathcal{N}^{th}_{X} & \equiv \bar{n}_X \sum_{X \to \ldots} BR(X \to 1 + 2 + \ldots) \prod_{k} \frac{n_k}{\bar{n}_k} \nonumber \\ \mathcal{N}^{th}_{XY} & \equiv \frac{\bar{n}_X}{\mathcal{B}^{eff}_{XY}} \sum_{X \to Y + \ldots} g_Y BR(X \to g_Y Y + 1 + \ldots) \left(\frac{n_Y}{\bar{n}_Y}\right)^{g_Y} \prod_{k} \frac{n_k}{\bar{n}_k} \nonumber \\ \mathcal{B}^{eff}_{XY} & \equiv \sum_{X \to Y + \ldots} g_Y BR(X \to g_Y Y + 1+\ldots) \nonumber\end{aligned}\end{split}\]

where \(g_Y\) is the \(Y\) multiplicity in the final state of \(X\) decays. In addition, defining:

\[x = \ln(R/R_0),\;\; N_i = \ln(n_i/s_0),\;\; {\rm and}\;\; N_S = \ln(S/S_0)\]

we can write Eqs.([Seq]), ([eq:nieq]) and ([eq:Rieq]) as:

\[\begin{split}\begin{aligned} N_S' & = \frac{e^{(3 x - N_S)}}{HT} \sum_{i} BR(i,X) \Gamma_i m_i \left(n_i - \mathcal{N}_{i}^{th} \right) \label{Seqb} \\ N_i' & = -3 + \frac{{\langle \sigma v \rangle}_i}{H} n_i [\left(\frac{\bar{n}_i}{n_i}\right)^2 -1] - \frac{\Gamma_i}{H} \frac{m_i}{R_i}\left(1 - \frac{\mathcal{N}_{i}^{th}}{n_i} \right) \nonumber \\ & + \sum_{a} \mathcal{B}_{ai}^{eff} \frac{\Gamma_a}{H} \frac{m_a}{R_a}\left(\frac{n_a}{n_i} - \frac{\mathcal{N}_{ai}^{th}}{n_i} \right) \\ R_i' & = -3 \frac{P_i}{n_i} + \sum_{a} \mathcal{B}_{ai}^{eff} \frac{\Gamma_a}{H} m_a \left( \frac{1}{2} - \frac{R_i}{R_a} \right) \left(\frac{n_a}{n_i} - \frac{\mathcal{N}_{ai}^{th}}{n_i} \right) \label{Nieq}\end{aligned}\end{split}\]

where \('=d/dx\).

The above equation for \(N_i\) also applies for coherent oscillating fields, if we define:

\[N_i = \ln(n_i/s_0),\;\; {\rm and}\;\; n_i \equiv \rho_i/m_i\]

so

\[\begin{split}\begin{aligned} N_i' & = -3 - \frac{\Gamma_i}{H} \nonumber \\ R_i'& = 0 \label{Nico}\end{aligned}\end{split}\]

where we assume that the coherent oscillating component does not couple to any of the other fields.

Collecting Eqs.([Seqb])-([Nieq]) and ([Nico]) we have a closed set of first order differential equations:

• Entropy:

  \[N_S' = \frac{e^{(3 x - N_S)}}{HT} \sum_{i} BR(i,X) \Gamma_i m_i \left(n_i - \mathcal{N}_{i}^{th} \right) \label{eq:Sfin}\]

• Thermal fields:

  \[\begin{split}\begin{aligned} N_i'& = -3 + \frac{{\langle \sigma v \rangle}_i}{H} n_i [\left(\frac{\bar{n}_i}{n_i}\right)^2 -1] - \frac{\Gamma_i}{H} \frac{m_i}{R_i}\left(1 - \frac{\mathcal{N}_{i}^{th}}{n_i} \right) \\ & + \sum_{a} \mathcal{B}_{ai}^{eff} \frac{\Gamma_a}{H} \frac{m_a}{R_a}\left(\frac{n_a}{n_i} - \frac{\mathcal{N}_{ai}^{th}}{n_i} \right) \nonumber \\ R_i' & = -3 \frac{P_i}{n_i} + \sum_{a} \mathcal{B}_{ai}^{eff} \frac{\Gamma_a}{H} m_a \left( \frac{1}{2} - \frac{R_i}{R_a} \right) \left(\frac{n_a}{n_i} - \frac{\mathcal{N}_{ai}^{th}}{n_i} \right)\end{aligned}\end{split}\]

• Coherent Oscillating fields:

  \[\begin{split}\begin{aligned} N_i' & = -3 - \frac{\Gamma_i}{H} \nonumber \\ R_i' & = 0 \label{eq:COeq}\end{aligned}\end{split}\]

As seen above, the equation for \(R_i = \rho_i/n_i\) depends on \(P_i/n_i\). A proper evaluation of this quantity requires knowledge of the distribution \(F_i(p,t)\). However, for relativistic (or massless) particles we have \(P_i = \rho_i/3\), as seen from Eq.([beqs]), while for particles at rest we have \(P_i = 0\). Hence \(F_i(p,t)\) is only required to evaluate the relativistic/non-relativistic transition, which corresponds to a relatively small part of the evolution history of particle \(i\). Nonetheless, to model this transition we approximate \(F_i\) by a thermal distribution and take \(T_i, \mu_i \ll m_i\), where \(T_i\) is the temperature of the particle (which can be different from the thermal bath’s). Under these approximations we have:

\[\begin{split}\begin{aligned} \frac{P_i}{n_i} & = T_i \nonumber \\ \frac{\rho_i}{n_i} & = T_i \left[ \frac{K_1(m_i/T_i)}{K_2(m_i/T_i)} \frac{m_i}{T_i} + 3 \right] \label{eq:p1}\end{aligned}\end{split}\]

where \(K_{1,2}\) are the modified Bessel functions. In particular, if \(m_i/T_i \gg 1\):

\[\frac{\rho_i}{n_i} \simeq T_i \left[\frac{3}{2} + \frac{m_i}{T_i} + 3 \right] \Rightarrow \frac{P_i}{n_i} = T_i = \frac{2 m_i}{3}\left( \frac{R_i}{m_i} -1 \right)\]

As shown above, for a given value of \(R_i = \rho_i/n_i\), Eq.([eq:p1]) can be inverted to compute \(T_i\) (\(=P_i/n_i\)):

\[\frac{P_i}{n_i} = T_i(R_i)\]

Since we are interested in the non-relativistic/relativistic transition, we can expand the above expression around \(R_i/m_i = 1\), so \(P_i/n_i\) can be written as:

\[\frac{P_i}{n_i} = \frac{2 m_i}{3}\left( \frac{R_i}{m_i} -1 \right) + m_i \sum_{n >1} a_n \left(\frac{R_i}{m_i} -1 \right)^n\]

where the coefficients \(a_n\) can be numerically computed from Eq.([eq:p1]). The above approximation should be valid for \(m_i/T_i \gtrsim 1\) (or \(R_i \gtrsim m_i\)). On the other hand, for \(m_i/T_i \ll 1\) (or \(R_i \gg m_i\)), we have the relativistic regime, with \(P_i/n_i = R_i/3\). Therefore we can approximate the \(P_i/n_i\) function for all values of \(R_i\) by:

\[\begin{split}\frac{P_i}{n_i} = \left\{ \begin{array}{rl} & \frac{2 m_i}{3}\left( \frac{R_i}{m_i} -1 \right) + m_i \sum_{n >1} a_n \left(\frac{R_i}{m_i} -1 \right)^n \mbox{ , for $R_i < \tilde{R}$} \\ & \frac{R_i}{3} \mbox{ , for $R_i > \tilde{R}$} \end{array} \right. \label{Pfin}\end{split}\]

where the coefficients \(a_n\) are given by the numerical fit of Eq.([eq:p1]) and \(\tilde{R}\) is given by the matching of the two solutions.

Finally, to solve Eqs.([eq:Sfin])-([eq:COeq]) we need to compute \(H\) according to Eq.([H]), which requires knowledge of the energy densities for all particles (\(\rho_i\)) and for the thermal bath (\(\rho_R\)). The former are directly obtained from \(N_i\) and \(R_i\), while the latter can be computed from \(N_S\):

\[T = \left(\frac{g_{*S}(T_R)}{g_{*S}(T)}\right)^{1/3} T_R \exp[N_S/3 -x] \Rightarrow \rho_R = \frac{\pi^2}{30} g_{*}(T) T^4\]

Eqs.([eq:Sfin])-([eq:COeq]), with the auxiliary equations for \(H\) (Eq.([H])) and \(P_i/n_i\) (Eq.([Pfin])) form a set of closed equations, which can be solved once the initial conditions for the number density (\(n_i\)), energy density (\(\rho_i\)) and entropy (\(S\)) are given. For thermal fluids we assume:

\[\begin{split}\begin{aligned} n_i(T_R) & = \left\{ \begin{array}{ll} 0 & , \mbox{ if ${\langle \sigma v \rangle}_i \bar{n}_i/H|_{T=T_R} < 10$} \\ \bar{n}_i(T_R) & , \mbox{ if ${\langle \sigma v \rangle}_i \bar{n}_i/H|_{T=T_R} > 10$} \end{array} \right. \label{ni0TP} \\ \frac{\rho_i}{n_i}(T_R) & = \frac{\bar{\rho}_i}{\bar{n}_i}(T_R)\end{aligned}\end{split}\]

where \(\bar{\rho}_i\) is the equilibrium energy density (with zero chemical potential) for the particle \(i\). While for coherent oscillating fluids the initial condition is set at the beginning of oscillations:

\[\begin{split}\begin{aligned} n_i(T^{osc}_i) & =\frac{\rho_i^{0}}{m_i(T^{osc}_i)} \\ \frac{\rho_i}{n_i}(T^{osc}_i) & = m_i\end{aligned}\end{split}\]

where \(T^{osc}_i\) is the oscillation temperature, given by \(3H(T^{osc}_i) = m_i(T^{osc}_i)\) and \(\rho_i^{0}\) the initial energy density for oscillations.

Finally, the initial condition for the entropy \(S\) is trivially obtained, once we assume a radiation dominated universe at \(T=T_R\):

\[S(T_R) = \frac{2 \pi^2}{45} g_*(T_R) T_R^3 R_0^3\]