pyBBN v0.9.0: particles/IntermediateParticle.py

#

For intermediate regime equilibrium particles, density, energy density and pressure are obtained through integration of distribution function

from common import integrators
import numpy


name = 'intermediate'

#

Particle density

$\begin{equation} g \int \frac{p^2 dp}{2 \pi^2} f\left( \frac{p}{T} \right) \end{equation}$

def density(particle):

    density, _ = integrators.integrate_1D(
        lambda p: (
            particle.equilibrium_distribution_function(
                particle.energy(p) / particle.T
            ) * p**2 * particle.dof / 2. / numpy.pi**2
        ), (particle.grid.MIN_MOMENTUM / particle.params.a,
            particle.grid.MAX_MOMENTUM / particle.params.a)
    )
    return density

#

Particle energy density

#

$\begin{equation} I_\rho = \frac{g}{2 \pi^2} p^2 E(p) \, f\left( \frac{E(p)}{T} \right) \end{equation}$

def energy_density_integrand(p, particle):

    E = particle.energy(p)
    return (
        particle.equilibrium_distribution_function(E / particle.T)
        * p**2 * E
        * particle.dof / 2. / numpy.pi**2
    )

#

$\begin{equation} \rho = \int dp I_\rho \end{equation}$

def energy_density(particle):

    energy_density, _ = integrators.integrate_1D(
        lambda p: energy_density_integrand(p, particle),
        (particle.grid.MIN_MOMENTUM / particle.params.a,
         particle.grid.MAX_MOMENTUM / particle.params.a)
    )
    return energy_density

#

Particle pressure

#

$\begin{equation} I_P = \frac{g}{6 \pi^2} \frac{p^4}{E(p)} \, f\left(\frac{E(p)}{T}\right) \end{equation}$

def pressure_integrand(p, particle):

    E = particle.energy(p)
    return (
        particle.equilibrium_distribution_function(E / particle.T)
        * p**4 / E
        * particle.dof / 6. / numpy.pi**2
    )

#

$\begin{equation} P = \int dp I_P \end{equation}$

def pressure(particle):

    pressure, _ = integrators.integrate_1D(
        lambda p: pressure_integrand(p, particle),
        (particle.grid.MIN_MOMENTUM / particle.params.a,
         particle.grid.MAX_MOMENTUM / particle.params.a)
    )
    return pressure

#

Master equation terms

#

$\begin{equation} \frac{M^2 x}{m^2 a T} * I(2) \end{equation}$

def numerator(particle):

    return (particle.mass**2 * particle.params.x
            / particle.params.m**2 / particle.aT
            * I(particle, 2))

#

$\begin{equation} \frac{I(4) + M_N^2 I(2)}{(a T)^2} \end{equation}$

def denominator(particle):

    return (I(particle, 4) + particle.conformal_mass**2 * I(particle)) / particle.aT**2

#

$\begin{equation} I(n) = \frac{g}{2 \pi^2} \int y^n dy \frac{ e^{-\frac{E_N(y)}{a T}} } { \left(e^{-\frac{E_N(y)}{a T}} + \eta \right)^2 } \end{equation}$

def I(particle, y_power=2):

    return particle.dof / 2. / numpy.pi**2 * integrators.integrate_1D(
        lambda y: (
            y**y_power * numpy.exp(-particle.conformal_energy(y) / particle.aT)
            / (numpy.exp(-particle.conformal_energy(y) / particle.aT) + particle.eta) ** 2
        ),
        (particle.grid.MIN_MOMENTUM, particle.grid.MAX_MOMENTUM)
    )[0]