Boltzmann Code Docs

Basic documentation for the python version of Boltzmann code.

Last update: Sep 20, 2017

Detailed Contents:

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 Eq. (1), 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\]

Installation

The code is written in Python and requires Python version 2.6 or later (but not Python 3) with the following external Python libraries:

• numpy
• scipy
• Assimulo >=2.9

In addition, some of the plotting functions may require:

• matplotlib
• pyROOT

For a successful installation of Assilumo, SUNDIALS must also be installed.

Running the code

A basic example of how to run the code is provided by Example.py. If all dependencies have been successfully installed, this code can be run as:

Example.py -p parameters.ini -o output.dat -P

The resulting output will be written to the output file output.dat. If the flag -P is set, a simple plot showing the evolution of the energy densities will also be displayed. The main arguments are described below:

usage: nonThermalDM.py [-h] [-p PARAMETERFILE] [-o OUTPUTFILE] [-P]

optional arguments:

-h, --help	show this help message and exit
-p PARAMETERFILE, --parameterFile PARAMETERFILE
	name of parameter file, where most options are defined
-o OUTPUTFILE, --outputFile OUTPUTFILE
	name of output file (optional argument). If not define, no output will be saved
-P, --plotResult
	show simple plot for the evolution of densities

Model Definitions

For each new BSM species, the user must supply the main functions required by the Boltzmann equations. These functions are supplied when defining the new state in the main file (using the Component class ) and correspond to:

• decays : specify the BRs and total width. Can be temperature-dependent.
• mass : specify the mass. Can be temperature-dependent.
• sigv : specify the thermally averaged annihilation cross-section. Can be temperature-dependent.

For the example provided these functions are given in modelDefinitions.py.

Furthermore, for each species, the user must specify if its spin set by the dof property \(dof = \pm 2*spin + 1\) (where the plus sign is for bosons and minus for fermions).

The Parameters File

The basic models parameters are set in the parameters file. These are model dependent and should be adapted for each individual input model. For the simple models assumed in nonThermalDM.py, the following parameters are defined:

• model parameters:
  • yCoupling (float): Mediator-DM-SM coupling
  • mMediator (float): Mediator mass (in GeV)
  • mDM (float): DM mass (in GeV)
  • TRH (float): reheat temperature (in GeV)
  • TF (float): Final temperature for evolution (should be after mediator decay)

The Output

If an output file is defined, all the output will be written to the file. This output includes the input parameters used, a summary of the main cosmological observables (relic densities, effective number of additional neutrinos,...). Furthermore, a table containing the values for the number and energy densities of each species defined as a function of the temperature (or scale factor) is also printed.

If no output file is defined, only the summary information will be written to the file. An example of a typical output is the following:

#-------------
# Parameters:
# mdm = 1.
# mmediator = 1000.
# tf = 1e-5
# trh = 100.
# ycoupling = 1e-10
#-------------
#-------------
# Summary
# TF=1e-05
# DM: T(osc)= None | T(decouple)= 100.0 | T(decay)= None | Omega h^2 (@TF) = 0.113472811107
# 
# Mediator: T(osc)= None | T(decouple)= None | T(decay)= 0.0475383532448 | Omega h^2 (@decay) = 5.53551558304e-29
# 
# Delta Neff (@TF) = 0.0
#-------------
#-------------
# Header:
#             R                      T (GeV)               n_{DM} (GeV^{3})        #rho_{DM} (GeV^{2})      n_{Mediator} (GeV^{3})  #rho_{Mediator} (GeV^{2})
        1.0000E+00                1.0000E+02                1.0901E-17                2.9445E-15                2.1794E+02                2.5433E+05       
        1.0000E+00                1.0000E+02                1.0954E-17                2.9757E-15                2.1794E+02                2.5433E+05       
        1.0000E+00                1.0000E+02                1.1008E-17                3.0069E-15                2.1794E+02                2.5433E+05       
        1.0000E+00                1.0000E+02                1.1062E-17                3.0381E-15                2.1794E+02                2.5433E+05       

...

Code Documentation

Example File

Example.main(parameterFile, outputFile, showPlot=True)

Main code to define the BSM contents and properties and solve the Boltzmann equations

Parameters:

• parameterFile – Path to the file defining the main model parameters
• outputFile – Path to the output file. If None, no results will be written.
• showPlot – If True, will show a simple plot for the evolution of the energy densities

modelDefinitions File

component File

boltzSolver File

boltzEqs File

AuxFuncs File

AuxDecays File

synopsis:This module provides decay classes and methods author:Andre Lessa <lessa.a.p@gmail.com>

class AuxDecays.Decay(instate=None, fstates=None, br=0.0)

Bases: object

Main class to store information about a specific decay.

Parameters:

• motherID – ID for the decaying mother
• fstateIDs – list of IDs for the daughters
• br – branching ratio value

compress(ID, IDecays)

Generates a new DecayList() from compressing the ID decays

class AuxDecays.DecayList

Bases: object

Main class to store information about all the decays of a specific species.

Parameters:

• dlist – List of Decay objects
• width – Total width
• Xfraction – fraction of energy injected in the thermal bath

addDecay(decay)

getAllFinalStates()

Search