#  A precision calculation of relic neutrino decoupling

Kensuke Akita1*** E-mail address<EMAIL_ADDRESS>and Masahide
Yamaguchi1††† E-mail address<EMAIL_ADDRESS>
1Department of Physics, Tokyo Institute of Technology, Tokyo 152-8551, Japan

( Abstract
We study the distortions of equilibrium spectra of relic neutrinos due to the
interactions with electrons, positrons, and neutrinos in the early Universe.
We solve the integro-differential kinetic equations for the neutrino density
matrix, including three-flavor oscillations and finite temperature corrections
from QED up to the next-to-leading order $\mathcal{O}(e^{3})$ for the first
time. In addition, the equivalent kinetic equations in the mass basis of
neutrinos are directly solved, and we numerically evaluate the distortions of
the neutrino spectra in the mass basis as well, which can be easily
extrapolated into those for non-relativistic neutrinos in the current
Universe. In both bases, we find the same value of the effective number of
neutrinos, $N_{\rm eff}=3.044$, which parameterizes the total neutrino energy
density. The estimated error for the value of $N_{\rm eff}$ due to the
numerical calculations and the choice of neutrino mixing parameters would be
at most $0.0005$. )

###### Contents

  1. 1 Introduction
  2. 2 Neutrino decoupling
    1. 2.1 Boltzmann equations in the flavor basis
    2. 2.2 Boltzmann equations in the mass basis
    3. 2.3 Finite temperature corrections from QED
    4. 2.4 Computational method and initial conditions
  3. 3 Results
    1. 3.1 The flavor basis
    2. 3.2 The mass basis
    3. 3.3 Transformation of distributions in the flavor and mass bases
  4. 4 Conclusions
  5. A Kinetic equations for neutrinos in comoving variables
  6. B Analytic estimation of the collision integral

##  1 Introduction

The successful hot big bang model after inflation predicts that neutrinos
produced in the early Universe still exist in the current Universe. These
relic neutrinos are confirmed indirectly by the observations of primordial
abundances of light elements from Big Bang Nucleosynthesis (BBN), the
anisotropies of the Cosmic Microwave Background (CMB) and the distribution of
Large Scale Structure (LSS) of the Universe.

The cosmic neutrino background was generated at high temperature and kept in
thermal equilibrium through weak interactions. When the temperature of the
Universe decreased, weak interactions became ineffective and cosmic neutrinos
were decoupled with other particles at the decoupling temperature $T_{\rm
dec}\sim 2\ {\rm MeV}$. In the instantaneous decoupling limit, the energy
spectrum of neutrinos takes a form of Fermi-Dirac distribution function and
receives only the effect of redshift of physical momentum after the
decoupling. Soon after the decoupling of neutrinos, electrons and positrons
start to annihilate and to heat photons when the temperature of the Universe
is almost equal to the electron mass $m_{e}=0.511\ {\rm MeV}$. If we assume
that electrons and positrons annihilate only into photons, we can
approximately estimate the ratio of the temperatures of cosmic photons and
neutrinos, $T_{\gamma}/T_{\nu}\simeq 1.40102$, using the entropy conservation
of the Universe.

However, the decoupling temperature of neutrinos $T_{\rm dec}$ and the
temperature of annihilation of electrons and positrons are so close that some
neutrinos keep interacting with electrons and positrons. These interaction
processes become more efficient for neutrinos with higher energies because the
interaction rates of relativistic particles with higher energies are larger[1,
2]. Due to these processes, non-thermal distortions in the neutrino spectra
are produced and the photon temperature increases less than that in the
instantaneous decoupling limit. In particular, the non-thermal distortions
increase the total energy density of neutrinos, which can be parameterized by
the effective number of neutrinos $N_{\rm eff}$. This parameter can be
constrained by cosmological observations such as the measurement of the CMB
anisotropies.

The non-instantaneous decoupling of neutrinos was studied numerically for a
long time. In particular, this numerical study requires solving integro-
differential kinetic equations, which correspond to the Boltzmann equations
for neutrino momentum distributions. 111Recently, various approximate
evaluations on $N_{\rm eff}$ are also developing. In refs. [19, 20], a simple
evaluation model of $N_{\rm eff}$ is proposed, including finite temperature
QED corrections and neutrino oscillations. In ref. [21], the authors
approximately estimate $N_{\rm eff}$ with including QED corrections up to
$\mathcal{O}(e^{4})$. First, several studies [3, 4, 5] solved these Boltzmann
equations under some approximations such as Maxwell-Boltzmann statistics
approximation for neutrinos. A few years later, the Boltzmann equations for
the distorted Fermi-Dirac statistics of neutrinos were solved in refs. [6, 7,
8, 9]. Finally, the kinetic equations were solved with including finite
temperature radiative corrections at leading order [10, 11, 12, 13, 14, 15]
and three-flavor neutrino oscillations [16, 17, 18]. The kinetic equations
including neutrino oscillations correspond to the Boltzmann equations for the
neutrino density matrix. In refs. [16, 18], the authors solved the Boltzmann
equations under the damping approximation, where the off-diagonal parts of the
collision terms are treated as the damping factors. In ref. [17], the
Boltzmann equations for the neutrino density matrix with the full collision
terms were solved.

In the present Universe, since the average magnitude of momenta for neutrinos
is $\langle p\rangle\sim 0.53~{}{\rm meV}\ll\sqrt{\Delta
m^{2}_{21}},\sqrt{|\Delta m^{2}_{31}|}$, two massive neutrinos at least are
non-relativistic. In the non-relativistic epoch for neutrinos, we cannot
quantize flavor neutrinos, which are flavor eigenstates of neutrinos, in the
conventional way as we quantize fields whose masses are diagonalized, and
hence the flavor neutrino spectra in this epoch do not make sense. In order to
investigate and to detect the neutrino spectra in the present epoch, we need
to consider massive neutrino spectra. Then, in this paper, by solving the
kinetic equations for massive neutrinos in the early Universe, we study the
distortions for neutrino spectra in the mass basis too, which could be easily
extrapolated to those in the current epoch, and compare the results both in
the flavor and mass bases. In addition, since we can expect the better
accuracy on future measurements of $N_{\rm eff}$, we solve the Boltzmann
equations for the neutrino density matrix with full collision terms, including
finite temperature corrections from QED up to the next-to-leading order
$\mathcal{O}(e^{3})$ for the first time.

This paper is organized as follows. In Sec. 2, we give the Boltzmann equations
both in the flavor and mass bases in order to analyze the decoupling process
of neutrinos. In this section, we also discuss finite temperature corrections
from QED and comment on the computational method and initial conditions. In
Sec. 3, we present our results of neutrino spectra and the value of $N_{\rm
eff}$ for each basis. We also discuss the relation of distribution functions
in the flavor and mass bases. Finally, we give our conclusions in Sec. 4. In
appendices, the kinetic equations for neutrinos in comoving variables and
analytic estimation of the collision integral are given.

##  2 Neutrino decoupling

###  2.1 Boltzmann equations in the flavor basis

In order to describe the process of neutrino decoupling in the early Universe,
in particular, to estimate the spectral distortion with good precision, we
first consider field operators of flavor neutrinos and their density matrices
in a homogeneous system. In ultra-relativistic limit, the field operators of
left-handed neutrinos are expanded as

$\displaystyle\bm{\nu}_{\alpha}(x)=\int\frac{d^{3}{\bm{p}}}{(2\pi)^{3}\sqrt{2p_{0}}}\left(a_{\alpha}({\bm{p}},t)u_{\bm{p}}e^{i\bm{p\cdot
x}}+b_{\alpha}^{{\dagger}}({\bm{p}},t)v_{\bm{p}}e^{-i\bm{p\cdot x}}\right),$
(1)

where $a_{\alpha}(\bm{p},t)=a_{\alpha}(\bm{p})e^{-ip_{0}t}$ and
$b_{\alpha}(\bm{p},t)=b_{\alpha}(\bm{p})e^{-ip_{0}t}$ are annihilation
operators for negative-helicity neutrinos and positive-helicity anti-
neutrinos, respectively. $\alpha$ and $\bm{p}$ are a flavor index and a three
dimensional momentum with $p_{0}\simeq|\bm{p}|$, respectively. $u_{\bm{p}}\
(v_{\bm{p}})$ denotes the Dirac spinor for a massless negative-helicity
particle (positive-helicity anti-particle), which is normalized to be
$u_{\bm{p}}^{\dagger}u_{\bm{p}}=v_{\bm{p}}^{\dagger}v_{\bm{p}}=2p_{0}$, and
satisfies

$\displaystyle/pu_{\bm{p}}=0,\ \ \ \ \ \ \ \ /pv_{\bm{p}}=0.$ (2)

These expansions of the field operators make sense only in the ultra-
relativistic limit. The annihilation and creation operators satisfy the anti-
commutation relations,

$\displaystyle\\{a_{\alpha}(\bm{p}),a_{\beta}^{{\dagger}}(\bm{p}^{\prime})\\}=\\{b_{\alpha}(\bm{p}),b_{\beta}^{{\dagger}}(\bm{p}^{\prime})\\}=\delta_{\alpha\beta}(2\pi)^{3}\delta^{(3)}(\bm{p}-\bm{p}^{\prime}).$
(3)

The density matrices for neutrinos and anti-neutrinos are defined through the
following expectation values of these operators with regard to the initial
thermal equilibrium states,

$\displaystyle\langle
a^{{\dagger}}_{\beta}(\bm{p},t)a_{\alpha}(\bm{p}^{\prime},t)\rangle$
$\displaystyle=(2\pi)^{3}\delta^{(3)}(\bm{p}-\bm{p^{\prime}})\left(\rho_{p}\right)_{\alpha\beta},$
$\displaystyle\langle
b^{{\dagger}}_{\alpha}(\bm{p},t)b_{\beta}(\bm{p^{\prime}},t)\rangle$
$\displaystyle=(2\pi)^{3}\delta^{(3)}(\bm{p}-\bm{p^{\prime}})\left(\bar{\rho}_{p}\right)_{\alpha\beta},$
(4)

where $p=|\bm{p}|$. Due to the reversed order of flavor indices in
$\bar{\rho}_{p}(t)$, both density matrices transform in the same way under a
unitary transformation of flavor space. Here the diagonal parts are the usual
distribution functions of flavor neutrinos and the off-diagonal parts are non-
zero in the presence of flavor mixing. Note again that since the off-diagonal
parts of the neutrino mass matrix are not zero, we cannot define annihilation
operators, creation operators and density matrices for non-relativistic flavor
neutrinos in the conventional way. In Sec. 2.2, we write field operators of
neutrinos including their masses in the mass basis so that we can define
density matrices for non-relativistic neutrinos.

The neutrino density matrix takes the following form,

$\displaystyle\rho_{p}(t)=\begin{pmatrix}\rho_{ee}&\rho_{e\mu}&\rho_{e\tau}\\\
\rho_{\mu e}&\rho_{\mu\mu}&\rho_{\mu\tau}\\\ \rho_{\tau
e}&\rho_{\tau\mu}&\rho_{\tau\tau}\end{pmatrix}=\begin{pmatrix}f_{\nu_{e}}&a_{1}+ia_{2}&b_{1}+ib_{2}\\\
a_{1}-ia_{2}&f_{\nu_{\mu}}&c_{1}+ic_{2}\\\
b_{1}-ib_{2}&c_{1}-ic_{2}&f_{\nu_{\tau}}\end{pmatrix},$ (5)

where $f_{\nu_{\alpha}}$ is the distribution function for flavor neutrinos and
the off-diagonal parts are characterized by the real parameters $a_{i},b_{i}$
and $c_{i}\ (i=1,2)$. Hereafter we neglect a neutrino asymmetry since neutrino
oscillations leading to flavor equilibrium before BBN impose a stringent
constraint on this asymmetry [22, 23, 24, 25, 26]. Under this assumption,
neutrinos and anti-neutrinos satisfy the same density matrices and the same
evolutions in the Universe, $\rho_{p}(t)=\bar{\rho}_{p}(t)^{\mathrm{T}}$.

The equations of motion for the neutrino density matrix in the expanding
Universe are [27, 28]

$\displaystyle(\partial_{t}-Hp\partial_{p})\rho_{p}(t)=-i\left[\left(\frac{M^{2}}{2p}-\frac{8\sqrt{2}G_{F}p}{3m_{W}^{2}}E\right),\
\rho_{p}(t)\right]+C[\rho_{p}(t)],$ (6)

where $H$ is the Hubble parameter, $G_{F}$ is the Fermi coupling constant,
$m_{W}$ is the W boson mass, and $[\cdot,\cdot]$ represents the commutator of
matrices with a flavor index. The first term222Only when we derive the first
term which comes from the free neutrino Hamiltonian including the mass matrix,
we replace the operators $a_{\alpha}(\bm{p},t)$ and $b_{\alpha}(\bm{p},t)$ in
Eq. (1) with $a_{\alpha}^{\rm
osc}(\bm{p},t)=(\exp(-i\Omega_{\bm{p}}t))_{\alpha\beta}a_{\beta}(\bm{p})$ and
$b_{\alpha}^{\rm
osc}(\bm{p},t)=(\exp(-i\Omega_{\bm{p}}t))_{\alpha\beta}b_{\beta}(\bm{p})$ as
in [27], where $\Omega_{\bm{p}}=\sqrt{\bm{p}^{2}+M^{2}}$.  in the commutator
is the vacuum oscillation term proportional to the mass-squared matrix in the
flavor basis $M^{2}$. The mass-squared matrix is related to the diagonal mass-
squared matrix in the mass basis $M_{\rm diag}^{2}$ through the Pontecorvo-
Maki-Nakagawa-Sakata matrix, assuming the CP conservation,

$\displaystyle U_{\rm PMNS}$
$\displaystyle\equiv\begin{pmatrix}U_{e1}&U_{e2}&U_{e3}\\\ U_{\mu 1}&U_{\mu
2}&U_{\mu 3}\\\ U_{\tau 1}&U_{\tau 2}&U_{\tau 3}\end{pmatrix},$
$\displaystyle=\begin{pmatrix}c_{12}c_{13}&s_{12}c_{13}&s_{13}\\\
-s_{12}c_{23}-c_{12}s_{23}s_{13}&c_{12}c_{23}-s_{12}s_{23}s_{13}&s_{23}c_{13}\\\
s_{12}s_{23}-c_{12}c_{23}s_{13}&-c_{12}s_{23}-s_{12}c_{23}s_{13}&c_{23}c_{13}\end{pmatrix},$
(7)

where $c_{ij}=\cos\theta_{ij}$ and $s_{ij}=\sin\theta_{ij}$ for $ij=12,\ 13$,
or $23$. The relation between the mass-squared matrices of neutrinos in the
two bases is

$\displaystyle M_{\rm diag}^{2}$ $\displaystyle=U_{\rm
PMNS}^{{\dagger}}M^{2}U_{\rm PMNS},$ $\displaystyle={\rm diag}(m_{1}^{2},\
m_{2}^{2},\ m_{3}^{2}).$ (8)

From the global analysis of neutrino oscillation experiments in [30], we
consider the following best-fit values of neutrino masses and mixing
parameters,

$\displaystyle\left(\frac{\Delta m^{2}_{21}}{10^{-5}\ {\rm eV^{2}}},\
\frac{\Delta m_{31}^{2}}{10^{-3}\ {\rm eV^{2}}},\ s_{12}^{2},\ s_{23}^{2},\
s_{13}^{2}\right)_{\rm NH}$ $\displaystyle=(7.39,\ 2.525,\ 0.310,\ 0.582,\
0.0224),$ $\displaystyle\left(\frac{\Delta m^{2}_{21}}{10^{-5}\ {\rm
eV^{2}}},\ \frac{\Delta m_{31}^{2}}{10^{-3}\ {\rm eV^{2}}},\ s_{12}^{2},\
s_{23}^{2},\ s_{13}^{2}\right)_{\rm IH}$ $\displaystyle=(7.39,\ -2.512,\
0.310,\ 0.582,\ 0.02263),$ (9)

where $\Delta m^{2}_{ij}=m^{2}_{i}-m^{2}_{j}$. The first (second) equation in
Eq. (9) corresponds to the normal (inverted) hierarchy ordering of neutrino
masses.

The second term in the commutator in Eq. (6) represents the refractive effect
in the medium which comes from one-loop thermal contributions to the neutrino
self-energy. The diagonal matrix $E$ is the energy density of the charged
leptons and in the temperature of ${\rm MeV}$ scale, $E$ takes the following
form approximately,

$\displaystyle E={\rm diag}(\rho_{ee},\ 0,\ 0),$ (10)

where $\rho_{ee}=\rho_{e^{-}}+\rho_{e^{+}}$ is the energy density of electrons
and positrons. We neglect other refractive terms coming from the charged
lepton asymmetries and neutrino self-interactions, which are significantly
suppressed [27, 29].

The final term in Eq. (6) represents the collisions of neutrinos with
electrons, positrons, and themselves, which are dominated by two-body
reactions $1+2\rightarrow 3+4$. As done in the previously most accurate
calculation of neutrino decoupling in the early Universe [17], we also deal
with both diagonal and off-diagonal collision terms for the processes which
involve electrons and positrons. On the other hand, we do not treat the off-
diagonal terms for the self-interactions of neutrinos, such as
$\nu\nu\leftrightarrow\nu\nu$ or $\nu\bar{\nu}\leftrightarrow\nu\bar{\nu}$
since the annihilations of electrons and positrons are important for the
heating process of neutrinos while the self-interactions of neutrinos less
contribute to this heating process.

The diagonal collision term from the self-interaction processes
$\nu(p_{1})\nu(p_{2})\leftrightarrow\nu(p_{3})\nu(p_{4})$ and
$\nu(p_{1})\bar{\nu}(p_{2})\leftrightarrow\nu(p_{3})\bar{\nu}(p_{4})$ is

$\displaystyle C_{S}[\nu_{\alpha}(p_{1})]$
$\displaystyle=\frac{2^{5}G_{F}^{2}}{2\left|\bm{p}_{1}\right|}\int\frac{d^{3}\bm{p}_{2}}{(2\pi)^{3}2\left|\bm{p}_{2}\right|}\frac{d^{3}\bm{p}_{3}}{(2\pi)^{3}2\left|\bm{p}_{3}\right|}\frac{d^{3}\bm{p}_{4}}{(2\pi)^{3}2\left|\bm{p}_{4}\right|}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle\times\Bigl{[}\left\\{4(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})+2(p_{1}\cdot p_{2})(p_{3}\cdot
p_{4})\right\\}F(\nu_{\alpha}^{(1)},\nu_{\alpha}^{(2)},\nu_{\alpha}^{(3)},\nu_{\alpha}^{(4)})$
$\displaystyle+\left\\{(p_{1}\cdot p_{4})(p_{2}\cdot p_{3})+(p_{1}\cdot
p_{2})(p_{3}\cdot
p_{4})\right\\}F(\nu_{\alpha}^{(1)},\nu_{\beta}^{(2)},\nu_{\alpha}^{(3)},\nu_{\beta}^{(4)})$
$\displaystyle+(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})F\left(\nu_{\alpha}^{(1)},\nu_{\alpha}^{(2)},\nu_{\beta}^{(3)},\nu_{\beta}^{(4)}\right)$
$\displaystyle+\left\\{(p_{1}\cdot p_{4})(p_{2}\cdot p_{3})+(p_{1}\cdot
p_{2})(p_{3}\cdot
p_{4})\right\\}F(\nu_{\alpha}^{(1)},\nu_{\gamma}^{(2)},\nu_{\alpha}^{(3)},\nu_{\gamma}^{(4)})$
$\displaystyle+(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})F(\nu_{\alpha}^{(1)},\nu_{\alpha}^{(2)},\nu_{\gamma}^{(3)},\nu_{\gamma}^{(4)})\Bigl{]},$
(11)

where $\alpha,\beta,\gamma=e,\mu,\tau$ and
$\alpha\neq\beta,\alpha\neq\gamma,\beta\neq\gamma$. We define
$F(\nu_{\alpha}^{(1)},\nu_{\beta}^{(2)},\nu_{\gamma}^{(3)},\nu_{\delta}^{(4)})$
as

$\displaystyle
F(\nu_{\alpha}^{(1)},\nu_{\beta}^{(2)},\nu_{\gamma}^{(3)},\nu_{\delta}^{(4)})$
$\displaystyle=f_{\nu_{\gamma}}(p_{3})f_{\nu_{\delta}}(p_{4})\left(1-f_{\nu_{\alpha}}(p_{1})\right)\left(1-f_{\nu_{\beta}}(p_{2})\right)$
$\displaystyle-
f_{\nu_{\alpha}}(p_{1})f_{\nu_{\beta}}(p_{2})\left(1-f_{\nu_{\gamma}}(p_{3})\right)\left(1-f_{\nu_{\delta}}(p_{4})\right),$
(12)

where $\alpha,\beta,\delta,\gamma=e,\mu,\tau$.

The collision term from the annihilation processes
$\nu(p_{1})\bar{\nu}(p_{2})\leftrightarrow e^{-}(p_{3})e^{+}(p_{4})$ is

$\displaystyle C_{A}$
$\displaystyle=\frac{1}{2}\frac{2^{5}G_{F}^{2}}{2\left|\bm{p}_{1}\right|}\int\frac{d^{3}\bm{p}_{2}}{(2\pi)^{3}2\left|\bm{p}_{2}\right|}\frac{d^{3}\bm{p}_{3}}{(2\pi)^{3}2E_{3}}\frac{d^{3}\bm{p}_{4}}{(2\pi)^{3}2E_{4}}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle\times\Bigl{[}4(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})F^{LL}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)$
$\displaystyle+4(p_{1}\cdot p_{3})(p_{2}\cdot
p_{4})F^{RR}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)$
$\displaystyle+2(p_{1}\cdot
p_{2})m_{e}^{2}\Bigl{(}F^{LR}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)+F^{RL}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)\Bigl{)}\Bigl{]},$
(13)

where

$\displaystyle
F^{ab}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)$
$\displaystyle=f_{e}(p_{3})f_{e}(p_{4})\Bigl{(}Y^{a}(1-\bar{\rho}_{2})Y^{b}(1-\rho_{1})+(1-\rho_{1})Y^{b}(1-\bar{\rho}_{2})Y^{a}\Bigl{)}$
$\displaystyle-(1-f_{e}(p_{3}))(1-f_{e}(p_{4}))\Bigl{(}Y^{a}\bar{\rho}_{2}Y^{b}\rho_{1}+\rho_{1}Y^{b}\bar{\rho}_{2}Y^{a}\Bigl{)}.$
(14)

Here $f_{e}(p)$ is the distribution functions of electrons and positrons, and
$\bar{\rho}_{2}=\rho_{2}^{\mathrm{T}}$. We assume that electrons and positrons
are always in thermal equilibrium since electrons, positrons and photons
interact with each other through rapid electromagnetic interactions. Under
this assumption, the distribution functions of electrons and positrons take
the following form,

$\displaystyle f_{e}(p)=\frac{1}{\exp(\sqrt{p^{2}+m_{e}^{2}}/T_{\gamma})+1}.$
(15)

$Y^{a}(a=L,R)$ is a $3\times 3$ matrix of couplings and becomes in the flavor
basis

$\displaystyle Y^{L}$ $\displaystyle={\rm
diag}(g_{L},\tilde{g}_{L},\tilde{g}_{L}),$ $\displaystyle Y^{R}$
$\displaystyle={\rm diag}(g_{R},g_{R},g_{R}),$ (16)

where

$\displaystyle g_{L}=\frac{1}{2}+\sin^{2}\theta_{W},\ \ \ \
\tilde{g}_{L}=-\frac{1}{2}+\sin^{2}\theta_{W},\ \ \ \
g_{R}=\sin^{2}\theta_{W}.$ (17)

Here $\sin^{2}\theta_{W}\simeq 0.231$ and $\theta_{W}$ is the weak mixing
angle.

The collision term from the scattering processes
$\nu(p_{1})e^{-}(p_{2})\leftrightarrow\nu(p_{3})e^{-}(p_{4})$ and
$\nu(p_{1})e^{+}(p_{2})\leftrightarrow\nu(p_{3})e^{+}(p_{4})$ is

$\displaystyle C_{SC}$
$\displaystyle=\frac{1}{2}\frac{2^{5}G_{F}^{2}}{2\left|\bm{p}_{1}\right|}\int\frac{d^{3}\bm{p}_{2}}{(2\pi)^{3}2E_{2}}\frac{d^{3}\bm{p}_{3}}{(2\pi)^{3}2\left|\bm{p}_{3}\right|}\frac{d^{3}\bm{p}_{4}}{(2\pi)^{3}2E_{4}}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle\times\Bigl{[}4\left\\{(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})+(p_{1}\cdot p_{2})(p_{3}\cdot p_{4})\right\\}$
$\displaystyle\times\Bigl{(}F^{LL}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)+F^{RR}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)\Bigl{)}$
$\displaystyle-4(p_{1}\cdot
p_{3})m_{e}^{2}\Bigl{(}F^{LR}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)+F^{RL}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)\Bigl{)}\Bigl{]},$
(18)

where

$\displaystyle F^{ab}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)$
$\displaystyle=f_{e}(p_{4})(1-f_{e}(p_{2}))\Bigl{(}Y^{a}\rho_{3}Y^{b}(1-\rho_{1})+(1-\rho_{1})Y^{b}\rho_{3}Y^{a}\Bigl{)}$
$\displaystyle-
f_{e}(p_{2})(1-f_{e}(p_{4}))\Bigl{(}\rho_{1}Y^{b}(1-\rho_{3})Y^{a}+Y^{a}(1-\rho_{3})Y^{b}\rho_{1}\Bigl{)}.$
(19)

These collision terms are described in detail in appendix A.

In addition to the Boltzmann equations for the neutrino density matrix, the
energy conservation law must be satisfied,

$\displaystyle\frac{d\rho}{dt}=-3H(\rho+P),$ (20)

where $\rho$ and $P$ are the total energy density and pressure of the standard
model particles ($\gamma,e^{\pm},\nu_{i}$) respectively. Though we will
discuss finite temperature corrections from QED to $\rho,P$ and $m_{e}$ later,
in the ideal gas limit, they are given as follows, which are denoted by
$\rho_{(0)}$ and $P_{(0)}$ respectively,

$\displaystyle\rho_{(0)}$
$\displaystyle=\frac{\pi^{2}T_{\gamma}^{4}}{15}+\frac{2}{\pi^{2}}\int\frac{dpp^{2}\sqrt{p^{2}+m_{e}^{2}}}{\exp(\sqrt{p^{2}+m_{e}^{2}}/T_{\gamma})+1}+\sum_{\alpha=e,\mu,\tau}\frac{1}{\pi^{2}}\int
dp~{}p^{3}f_{\nu_{\alpha}}(p),$ $\displaystyle P_{(0)}$
$\displaystyle=\frac{\pi^{2}T_{\gamma}^{4}}{45}+\frac{2}{\pi^{2}}\int\frac{dpp^{4}}{3\sqrt{p^{2}+m_{e}^{2}}[\exp(\sqrt{p^{2}+m_{e}^{2}}/T_{\gamma})+1]}+\sum_{\alpha=e,\mu,\tau}\frac{1}{3\pi^{2}}\int
dp~{}p^{3}f_{\nu_{\alpha}}(p).$ (21)

The energy conservation law governs the evolution of the photon temperature
$T_{\gamma}$. The Hubble parameter in Eqs. (6) and (20) is calculated using
the usual relation, $3H^{2}m_{\rm Pl}^{2}=8\pi\rho$ with $m_{\rm Pl}$ being
the Planck mass, where we ignore the curvature term and the cosmological
constant because they are negligible in the radiation dominated epoch.

###  2.2 Boltzmann equations in the mass basis

In this section, we formulate the Boltzmann equations for the density matrices
of massive neutrinos at the early Universe in the mass basis. If we would like
to observe the distortions of neutrinos in the current Universe in future, it
is easier to follow the evolution of negative-helicity neutrinos in the
diagonal mass basis since the helicity states of neutrinos are conserved while
non-relativistic neutrinos are freely streaming. Thus, it is quite useful to
formulate the Boltzmann equations for negative-helicity neutrinos in the mass
basis though we concentrate on neutrino decoupling processes in this paper,
where ultra-relativistic limit is a good approximation and there are no much
difference between the two bases. This approach is also complementary to that
in the flavor basis given in the previous subsection 2.1 and is useful for the
cross-check of the results.

Since the negative-helicity neutrinos in the mass basis satisfy the free Dirac
equation, they are expanded as

$\displaystyle\bm{\nu}_{i}(x)=\int\frac{d^{3}{\bm{p}}}{(2\pi)^{3}\sqrt{2E_{i}}}\left(a_{i}({\bm{p}},t)u^{(i)}_{\bm{p}}e^{i\bm{p\cdot
x}}+b_{i}^{{\dagger}}({\bm{p}},t)v^{(i)}_{\bm{p}}e^{-i\bm{p\cdot x}}\right),$
(22)

where $i(=1,2,3)$ represents a mass eigenstate,
$a_{i}({\bm{p}},t)=a_{i}({\bm{p}})e^{-iE_{i}t},b_{i}({\bm{p}},t)=b_{i}({\bm{p}})e^{-iE_{i}t}$,
$E_{i}=\sqrt{\bm{p}^{2}+m_{i}^{2}}$ and $m_{i}$ is the neutrino mass in the
mass basis. Since $u^{(i)}_{\bm{p}}\ (v^{(i)}_{\bm{p}})$ denotes the Dirac
spinor for massive negative-helicity particles (positive-helicity anti-
particles), which is also normalized to be
$u^{(i)}_{\bm{p}}{}^{\dagger}u^{(i)}_{\bm{p}}=v^{(i)}_{\bm{p}}{}^{\dagger}v^{(i)}_{\bm{p}}=2E_{i}$,
the Dirac spinors satisfy

$\displaystyle(/p-m_{i})u^{(i)}_{\bm{p}}=0,\ \ \ \ \ \ \ \
(/p+m_{i})v^{(i)}_{\bm{p}}=0.$ (23)

As in the flavor basis, $a_{i}(\bm{p})$ and $b_{i}(\bm{p})$ are annihilation
operators for negative-helicity neutrinos and for positive-helicity anti-
neutrinos in the mass basis, respectively, which satisfy

$\displaystyle\\{a_{i}(\bm{p}),a_{j}^{{\dagger}}(\bm{p}^{\prime})\\}=\\{b_{i}(\bm{p}),b_{j}^{{\dagger}}(\bm{p}^{\prime})\\}=\delta_{ij}(2\pi)^{3}\delta^{(3)}(\bm{p}-\bm{p}^{\prime}).$
(24)

The density matrices for neutrinos and anti-neutrinos in the mass basis are
given by

$\displaystyle\langle
a^{{\dagger}}_{j}(\bm{p},t)a_{i}(\bm{p}^{\prime},t)\rangle$
$\displaystyle=(2\pi)^{3}\delta^{(3)}(\bm{p}-\bm{p^{\prime}})\left(\rho_{p}\right)_{ij},$
$\displaystyle\langle
b^{{\dagger}}_{i}(\bm{p},t)b_{j}(\bm{p^{\prime}},t)\rangle$
$\displaystyle=(2\pi)^{3}\delta^{(3)}(\bm{p}-\bm{p^{\prime}})\left(\bar{\rho}_{p}\right)_{ij},$
(25)

where the diagonal parts are the distribution functions for massive neutrinos.

In the following, we study the kinetic equations for the neutrino density
matrix in ultra-relativistic limit. In this limit, negative-helicity neutrinos
coincide with left-handed neutrinos due to no distinction between helicity and
chirality. The diagonalization of mass matrix for left-handed neutrinos in the
flavor basis is achieved through the transformations,

$\displaystyle\bm{\nu}_{\alpha}(x)=\sum_{i=1}^{3}U_{\alpha i}\bm{\nu}_{i}(x),$
(26)

with $\alpha=e,\mu,\tau$. Here $U_{\alpha i}$ represents a component of the
unitary matrix $U_{\rm PMNS}$ given in Eq. (2.1).

In order to specify the collision processes for massive neutrinos, we discuss
weak neutral currents and charged currents in the mass basis. The weak neutral
currents of electrons, positrons, and neutrinos are

$\displaystyle J_{\rm NC}^{\mu}=J_{ee}^{L\mu}+J_{ee}^{R\mu}+J_{\nu\nu}^{\mu},$
(27)

where $J_{ee}^{L\mu}$ and $J_{ee}^{R\mu}$ are the neutral currents of left-
handed electrons and right-handed electrons respectively and given by

$\displaystyle J_{ee}^{L\mu}$
$\displaystyle=\tilde{g}_{L}\bar{\bm{e}}\gamma^{\mu}(1-\gamma_{5})\bm{e},$
$\displaystyle J_{ee}^{R\mu}$
$\displaystyle=g_{R}\bar{\bm{e}}\gamma^{\mu}(1+\gamma_{5})\bm{e},$ (28)

where $\bm{e}$ is the field operator of electron. The neutral currents for the
left-handed neutrinos in both bases are given by the following form and they
are related through the unitary $U_{\rm PMNS}$ matrix,

$\displaystyle
J_{\nu\nu}^{\mu}=\sum_{\alpha=e,\mu,\tau}\bar{\bm{\nu}}_{\alpha}\gamma^{\mu}(1-\gamma_{5})\bm{\nu}_{\alpha}=\sum_{i=1,2,3}\bar{\bm{\nu}}_{i}\gamma^{\mu}(1-\gamma_{5})\bm{\nu}_{i}.$
(29)

The charged currents for left-handed electrons and left-handed electron
neutrinos are given by

$\displaystyle
J_{e\nu_{e}}^{\mu}=\bar{\bm{\nu}}_{e}\gamma^{\mu}(1-\gamma_{5})\bm{e}.$ (30)

Using the Fierz transformations for fermionic fields, we can describe the
Hamiltonian density including the charged currents at the neutrino decoupling
process as

$\displaystyle\mathcal{H}_{CC}=\frac{G_{F}}{\sqrt{2}}J^{{\dagger}\mu}_{e\nu_{e}}(J_{e\nu_{e}})_{\mu}$
$\displaystyle=\frac{G_{F}}{\sqrt{2}}J_{ee}^{L\mu}(J_{\nu_{e}\nu_{e}})_{\mu},$
$\displaystyle J_{\nu_{e}\nu_{e}}^{\mu}$
$\displaystyle=\bar{\bm{\nu}}_{e}\gamma^{\mu}(1-\gamma_{5})\bm{\nu}_{e}.$ (31)

Since Eq. (31) implies that the charged currents in the Hamiltonian density
$\mathcal{H}_{CC}$ can be replaced by the equivalent neutral currents, we
implicitly take the charged currents into account by making the following
replacement of the coefficient of neutral currents for electron neutrinos in
the collision terms,

$\displaystyle\tilde{g}_{L}\rightarrow g_{L}=\tilde{g}_{L}+1.$ (32)

In the mass basis, the charged currents for neutrinos and electrons are,
through the $U_{\rm PMNS}$ matrix,

$\displaystyle J_{e\nu_{e}}^{\mu}$
$\displaystyle=\bar{\bm{\nu}}_{e}\gamma^{\mu}(1-\gamma_{5})\bm{e},$
$\displaystyle=\sum_{i=1}^{3}U^{*}_{ei}\bar{\bm{\nu}}_{i}\gamma^{\mu}(1-\gamma_{5})\bm{e}.$
(33)

The replacement in Eq. (31) corresponds to the following relation in the mass
basis,

$\displaystyle J^{{\dagger}\mu}_{e\nu_{e}}(J_{e\nu_{e}})_{\mu}$
$\displaystyle=\sum_{i=1}^{3}\sum_{j=1}^{3}U_{ei}^{*}U_{ej}J_{ee}^{L\mu}(J_{{\nu}_{i}{\nu}_{j}})_{\mu},$
$\displaystyle J_{\nu_{i}\nu_{j}}^{\mu}$
$\displaystyle=\bar{\bm{\nu}}_{i}\gamma^{\mu}(1-\gamma_{5})\bm{\nu}_{j}.$ (34)

From the above equation, the corresponding couplings to $Y^{a}$ in Eqs. (14)
and (19) are changed into

$\displaystyle Y^{L}$ $\displaystyle\rightarrow
Z^{L}=\begin{pmatrix}\tilde{g}_{L}+U_{e1}^{*}U_{e1}&U_{e1}^{*}U_{e2}&U_{e1}^{*}U_{e3}\\\
U_{e2}^{*}U_{e1}&\tilde{g}_{L}+U_{e2}^{*}U_{e2}&U_{e2}^{*}U_{e3}\\\
U_{e3}^{*}U_{e1}&U_{e3}^{*}U_{e2}&\tilde{g}_{L}+U_{e3}^{*}U_{e3}\end{pmatrix},$
$\displaystyle Y^{R}$ $\displaystyle\rightarrow Z^{R}=Y^{R}={\rm
diag}(g_{R},g_{R},g_{R}).$ (35)

The collision term from the self-interaction processes in the mass basis takes
the same form as that in the flavor basis except for the subscripts,

$\displaystyle C_{S}[\nu_{i}(p_{1})]$
$\displaystyle=\frac{2^{5}G_{F}^{2}}{2\left|\bm{p}_{1}\right|}\int\frac{d^{3}\bm{p}_{2}}{(2\pi)^{3}2\left|\bm{p}_{2}\right|}\frac{d^{3}\bm{p}_{3}}{(2\pi)^{3}2\left|\bm{p}_{3}\right|}\frac{d^{3}\bm{p}_{4}}{(2\pi)^{3}2\left|\bm{p}_{4}\right|}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle\times\Bigl{[}\left\\{4(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})+2(p_{1}\cdot p_{2})(p_{3}\cdot
p_{4})\right\\}F(\nu_{i}^{(1)},\nu_{i}^{(2)},\nu_{i}^{(3)},\nu_{i}^{(4)})$
$\displaystyle+\left\\{(p_{1}\cdot p_{4})(p_{2}\cdot p_{3})+(p_{1}\cdot
p_{2})(p_{3}\cdot
p_{4})\right\\}F(\nu_{i}^{(1)},\nu_{j}^{(2)},\nu_{i}^{(3)},\nu_{j}^{(4)})$
$\displaystyle+(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})F\left(\nu_{i}^{(1)},\nu_{i}^{(2)},\nu_{j}^{(3)},\nu_{j}^{(4)}\right)$
$\displaystyle+\left\\{(p_{1}\cdot p_{4})(p_{2}\cdot p_{3})+(p_{1}\cdot
p_{2})(p_{3}\cdot
p_{4})\right\\}F(\nu_{i}^{(1)},\nu_{k}^{(2)},\nu_{i}^{(3)},\nu_{k}^{(4)})$
$\displaystyle+(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})F(\nu_{i}^{(1)},\nu_{i}^{(2)},\nu_{k}^{(3)},\nu_{k}^{(4)})\Bigl{]},$
(36)

where $i,j,k=1,2,3$ and $i\neq j,i\neq k,j\neq k$. The collision term from the
annihilation processes is

$\displaystyle C_{A}$
$\displaystyle=\frac{1}{2}\frac{2^{5}G_{F}^{2}}{2\left|\bm{p}_{1}\right|}\int\frac{d^{3}\bm{p}_{2}}{(2\pi)^{3}2\left|\bm{p}_{2}\right|}\frac{d^{3}\bm{p}_{3}}{(2\pi)^{3}2E_{3}}\frac{d^{3}\bm{p}_{4}}{(2\pi)^{3}2E_{4}}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle\times\Bigl{[}4(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})G^{LL}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)$
$\displaystyle+4(p_{1}\cdot p_{3})(p_{2}\cdot
p_{4})G^{RR}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)$
$\displaystyle+2(p_{1}\cdot
p_{2})m_{e}^{2}\Bigl{(}G^{LR}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)+G^{RL}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)\Bigl{)}\Bigl{]},$
(37)

where

$\displaystyle
G^{ab}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)$
$\displaystyle=f_{e}(p_{3})f_{e}(p_{4})\Bigl{(}Z^{a}(1-\bar{\rho}_{2})Z^{b}(1-\rho_{1})+(1-\rho_{1})Z^{b}(1-\bar{\rho}_{2})Z^{a}\Bigl{)}$
$\displaystyle-(1-f_{e}(p_{3}))(1-f_{e}(p_{4}))\Bigl{(}Z^{a}\bar{\rho}_{2}Z^{b}\rho_{1}+\rho_{1}Z^{b}\bar{\rho}_{2}Z^{a}\Bigl{)}.$
(38)

The collision term from the scattering processes is

$\displaystyle C_{SC}$
$\displaystyle=\frac{1}{2}\frac{2^{5}G_{F}^{2}}{2\left|\bm{p}_{1}\right|}\int\frac{d^{3}\bm{p}_{2}}{(2\pi)^{3}2E_{2}}\frac{d^{3}\bm{p}_{3}}{(2\pi)^{3}2\left|\bm{p}_{3}\right|}\frac{d^{3}\bm{p}_{4}}{(2\pi)^{3}2E_{4}}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle\times\Bigl{[}4\left\\{(p_{1}\cdot p_{4})(p_{2}\cdot
p_{3})+(p_{1}\cdot p_{2})(p_{3}\cdot p_{4})\right\\}$
$\displaystyle\times\Bigl{(}G^{LL}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)+G^{RR}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)\Bigl{)}$
$\displaystyle-4(p_{1}\cdot
p_{3})m_{e}^{2}\Bigl{(}G^{LR}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)+G^{RL}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)\Bigl{)}\Bigl{]},$
(39)

where

$\displaystyle G^{ab}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)$
$\displaystyle=f_{e}(p_{4})(1-f_{e}(p_{2}))\Bigl{(}Z^{a}\rho_{3}Z^{b}(1-\rho_{1})+(1-\rho_{1})Z^{b}\rho_{3}Z^{a}\Bigl{)}$
$\displaystyle-
f_{e}(p_{2})(1-f_{e}(p_{4}))\Bigl{(}\rho_{1}Z^{b}(1-\rho_{3})Z^{a}+Z^{a}(1-\rho_{3})Z^{b}\rho_{1}\Bigl{)}.$
(40)

Finally, the equations of motion for the neutrino density matrix in the mass
basis are given by

$\displaystyle(\partial_{t}-Hp\partial_{p})\rho_{p}(t)=-i\left[\left(\frac{M_{\rm
diag}^{2}}{2p}-\frac{8\sqrt{2}G_{F}p}{3m_{W}^{2}}\tilde{E}\right),\
\rho_{p}(t)\right]+C[\rho_{p}(t)].$ (41)

Since $\tilde{E}$ and $E$ are the thermal contributions to the self-energies
of left-handed neutrinos, these have the same relation as that of $M^{2}_{\rm
diag}$ and $M^{2}$, which is given by,

$\displaystyle\tilde{E}=U_{\rm PMNS}^{{\dagger}}EU_{\rm PMNS}.$ (42)

###  2.3 Finite temperature corrections from QED

In this section, we discuss finite temperature corrections from QED up to the
next-to-leading order $\mathcal{O}(e^{3})$, which modify electron, positron
and photon masses. Through these corrections, several points in the former
sections are changed. First, the modification of masses affects the energy
density and the pressure of the electromagnetic plasma, and the collision
rates involving electrons and positrons. In addition, the expansion rate $H$
in the Boltzmann equations changes through the total energy density of the
plasma.

The corrections to electron, positron and photon masses can be obtained
perturbatively by calculating the loop corrections to the self-energies of
these particles. The corrections to the electron and positron masses from
finite temperature effects at $\mathcal{O}(e^{2})$ are given by [31],

$\displaystyle\delta m_{e(2)}^{2}(p,T_{\gamma})$
$\displaystyle=\frac{2\pi\alpha
T_{\gamma}^{2}}{3}+\frac{4\alpha}{\pi}\int_{0}^{\infty}dk\frac{k^{2}}{E_{k}}\frac{1}{\exp(E_{k}/T_{\gamma})+1}$
$\displaystyle-\frac{2m_{e}^{2}\alpha}{\pi
p}\int^{\infty}_{0}dk\frac{k}{E_{k}}\log\left|\frac{p+k}{p-k}\right|\frac{1}{\exp(E_{k}/T_{\gamma})+1},$
(43)

where $\alpha=e^{2}/4\pi$ and $E_{k}=\sqrt{k^{2}+m_{e}^{2}}$. The last term
gives less than a $10\%$ correction to $\delta m_{e(2)}^{2}$ around the
decoupling temperature and the average momentum of electron [32], and
contributes about $-0.00005$ to $N_{\rm eff}$ in the instantaneous decoupling
limit [21]. Due to this smallness, we neglect the last term and consider only
the first two terms, which depend only on $T_{\gamma}$. On the other hand, the
thermal corrections to the photon mass at $\mathcal{O}(e^{2})$ are given by
[10],

$\displaystyle\delta
m_{\gamma(2)}^{2}(T_{\gamma})=\frac{8\alpha}{\pi}\int^{\infty}_{0}dk\frac{k^{2}}{E_{k}}\frac{1}{\exp(E_{k}/T_{\gamma})+1}.$
(44)

The total pressure and the total energy density of the electromagnetic plasma
are given by, including thermal mass corrections of electrons, positrons and
photons,

$\displaystyle P$
$\displaystyle=\frac{T_{\gamma}}{\pi^{2}}\int^{\infty}_{0}dk~{}k^{2}\log\left[\frac{(1+e^{-E_{e}/T_{\gamma}})^{2}}{(1-e^{-E_{\gamma}/T_{\gamma}})}\right],$
$\displaystyle\rho$ $\displaystyle=-P+T_{\gamma}\frac{dP}{dT_{\gamma}},$ (45)

where $E_{\gamma}=\sqrt{k^{2}+\delta m_{\gamma}^{2}}$ and
$E_{e}=\sqrt{k^{2}+m_{e}^{2}+\delta m_{e}^{2}}$. $\delta m_{\gamma}^{2}$ and
$\delta m_{e}^{2}$ denote the thermal mass corrections of photons and
electrons, respectively. We expand $P$ in terms of $\delta m_{e}^{2}$ and
$\delta m_{\gamma}^{2}$ at $\mathcal{O}(e^{2})$ and get the leading order
correction to the pressure [11],

$\displaystyle
P_{(2)}=-\int^{\infty}_{0}\frac{dk}{2\pi^{2}}\left[\frac{k^{2}}{E_{k}}\frac{\delta
m_{e(2)}^{2}}{\exp(E_{k}/T_{\gamma})+1}+\frac{k}{2}\frac{\delta
m_{\gamma(2)}^{2}}{\exp(k/T_{\gamma})-1}\right].$ (46)

Here we need to introduce the symmetry factor $1/2$ in Eq. (46) in order to
avoid the double counting of the thermal corrections to the total pressure.
Then the leading order correction to the energy density is obtained as

$\displaystyle\rho_{(2)}=-P_{(2)}+T_{\gamma}\frac{dP_{(2)}}{dT_{\gamma}}.$
(47)

The next-to-leading order of thermal corrections to the photon mass is
$\mathcal{O}(e^{3})$. These nontrivial corrections to the photon mass come
from the resummation of ring diagrams at all orders. Through this mass
correction, the thermal corrections to the pressure and energy density are
given by [21],

$\displaystyle P_{(3)}=\frac{e^{3}T_{\gamma}}{12\pi^{4}}I^{3/2}(T_{\gamma}),$
$\displaystyle\rho_{(3)}=\frac{e^{3}T_{\gamma}^{2}}{8\pi^{4}}I^{1/2}\partial_{T_{\gamma}}I,$
(48)

where

$\displaystyle
I(T_{\gamma})=2\int^{\infty}_{0}dk\left(\frac{k^{2}+E_{k}^{2}}{E_{k}}\right)\frac{1}{\exp(E_{k}/T_{\gamma})+1}.$
(49)

Note that the thermal corrections at $\mathcal{O}(e^{3})$ do not modify the
collision terms since these corrections change only the photon mass while the
next corrections to the electron mass would appear at $\mathcal{O}(e^{4})$.
Finally, we read the total energy density and the total pressure of the
electromagnetic plasma as

$\displaystyle P$ $\displaystyle=P_{(0)}+P_{(2)}+P_{(3)},$ $\displaystyle\rho$
$\displaystyle=\rho_{(0)}+\rho_{(2)}+\rho_{(3)}.$ (50)

###  2.4 Computational method and initial conditions

We solve a set of Eqs. (6) and (20) with the following comoving variables
instead of the cosmic time $t$, a momentum $p$, and the photon temperature
$T_{\gamma}$,

$\displaystyle x=m_{e}a,\ \ \ \ \ \ \ y=pa,\ \ \ \ \ \ \ \ z=T_{\gamma}a,$
(51)

where we take an arbitrary mass scale to be the electron mass $m_{e}$ and $a$
is the scale factor of the Universe, normalized as $z\rightarrow 1$ in the
high temperature limit.

Since the Boltzmann equations (6) are integro-differential equations due to
the collision terms, these equations were solved by a discretization in a grid
of comoving momenta $y_{i}$ in refs. [6, 7, 8, 13, 16, 17, 18], by an
expansion of the distortions of neutrinos in moments in refs. [9, 11, 12], or
by a hybrid method combining the former two methods in ref. [15]. In this
study, we use the former discretization method and take 100 grid points for
the comoving momentum, equally spaced in the region $y_{i}\in[0.02,\ 20]$.

We have numerically calculated the evolution of the density matrix and the
photon temperature in the interval $x_{\rm in}\leq x\leq x_{f}$. We have
chosen $x_{\rm in}=m_{e}/10\ {\rm MeV}$ as an initial time. Since neutrinos
keep in thermal equilibrium with the electromagnetic plasma at $x_{\rm in}$,
the initial values of density matrix $\rho_{y_{i}}^{\rm in}(x)$ are,

$\displaystyle\rho_{y_{i}}^{\rm in}(x)={\rm
diag}\left(\frac{1}{e^{y_{i}/z_{\rm in}}+1},\frac{1}{e^{y_{i}/z_{\rm
in}}+1},\frac{1}{e^{y_{i}/z_{\rm in}}+1}\right).$ (52)

The initial value of the dimensionless photon temperature at $x_{\rm in}$,
$z_{\rm in}$, slightly differs from unity because of the finite electron and
positron masses. Due to the entropy conservation of electromagnetic plasma,
neutrinos and anti-neutrinos, $z_{\rm in}$ is set as in [8],

$\displaystyle z_{\rm in}=1.00003.$ (53)

We set $x_{f}=30$ as a final time, when the neutrino density matrix and $z$
can be regarded as frozen.

##  3 Results

###  3.1 The flavor basis

First, we have numerically solved a set of Eqs. (6) and (20) in the flavor
basis, during the process of neutrino decoupling. In order to compare with
previous results, we discuss the cases with and without neutrino mixing, and
those with and without finite temperature corrections from QED up to
$\mathcal{O}(e^{2})$ and $\mathcal{O}(e^{3})$. In the case with neutrino
mixing, we also consider the normal neutrino mass hierarchy with the latest
best-fit values of neutrino mixing parameters. Although we have considered the
inverted mass hierarchy too, the results are almost the same in the case of
normal mass hierarchy as in ref. [17]. Hereafter we only show the case of the
normal mass hierarchy.

In Figs. 1 and 2, we show the evolution of the photon temperature and the
distortions of the flavor neutrino spectra for a comoving momentum $(y=5)$ in
the case with QED corrections up to $\mathcal{O}(e^{3})$ respectively, where
we plot the comoving photon temperature $z(x)$ and the neutrino spectra
$f_{\nu_{\alpha}}/f_{\rm eq}\ (f_{\rm eq}=[\exp(y)+1]^{-1})$ as a function of
the normalized cosmic scale factor $x$. At high temperature with $(x\lesssim
0.2)$, the temperature differences between photons and neutrinos are
negligible and neutrinos are in thermal equilibrium with electrons and
positrons. In the intermediate regime with $(0.2\lesssim x\lesssim 4)$, weak
interactions gradually become ineffective with shifting from small to large
momenta. In this period, the neutrino spectra are distorted since the energies
of electrons and positrons partially convert into those of neutrinos coupled
with electromagnetic plasma. Finally, at low temperature with $(x>4)$, the
collision term $C[\rho_{p}(t)]$ becomes ineffective and the distortions are
frozen.

In Fig. 2, we show the results of the two cases with and without neutrino
mixing. We find that the final values of $f_{\nu_{\alpha}}[y=5]$ without
flavor mixing are $1.17\%$ for $\nu_{e}$ and $0.500\%$ for $\nu_{\mu,\tau}$
larger than those in the instantaneous decoupling limit, that is, $f_{\rm
eq}[y=5]$. This difference between electron-type neutrinos and mu(tau)-type
neutrinos arises from the fact that only electron-type neutrinos interact with
electrons and positrons through the weak charged-currents. On the other hand,
in the cases with neutrino mixing, neutrino oscillations mix the distortions
of the flavor neutrinos too. Though the flavor oscillation effects are
subdominant to the refractive effects in the period with $(x\lesssim 0.2)$,
the refractive term gets ineffective in the lower temperature due to the
annihilations of electrons and positrons. Thus, the oscillation terms finally
mix the flavor neutrino distortions. We also find that the final values of
$f_{\nu_{\alpha}}[y=5]$ with flavor mixing are $0.895\%$ for $\nu_{e}$,
$0.648\%$ for $\nu_{\mu}$ and $0.663\%$ for $\nu_{\tau}$ larger than those in
the instantaneous decoupling limit, that is, $f_{\rm eq}[y=5]$.

In Fig. 3, we show the frozen values of the flavor neutrino spectra
$f_{\nu_{\alpha}}/f_{\rm eq}$ as a function of a comoving momentum $y$ for
both cases with and without neutrino mixing, including QED corrections up to
$\mathcal{O}(e^{3})$. This figure shows the fact that neutrinos with higher
energies interact with electrons and positrons until a later epoch. In
addition, we see that neutrino oscillations tend to equilibrate the flavor
neutrino distortions. Although the neutrino spectra $f_{\nu_{\alpha}}/f_{\rm
eq}$ with low energies are very slightly less than unity, these extractions of
low energy neutrinos stem from an energy boost through the scattering by
electrons, positrons, (and neutrinos) with sufficiently high energies, which
are not yet annihilated and hence still effective at the neutrino decoupling
process.

We also give several important quantities charactering the decoupling process
of neutrinos. In Tables. 1 and 2, we give final values (at $x_{f}=30$) of the
dimensionless photon temperature $z_{\rm fin}$, the difference of energy
densities and number densities of flavor neutrinos from those in the
instantaneous decoupling limit denoted by $\rho_{\nu_{0}}$ and $n_{\nu_{0}}$,
and the effective number of neutrinos $N_{\rm eff}$ defined as

$\displaystyle\rho_{r}=\left[1+\frac{7}{8}\left(\frac{4}{11}\right)^{4/3}N_{\rm
eff}\right]\rho_{\gamma}.$ (54)

Here $\rho_{r}$ and $\rho_{\gamma}$ are the energy densities of the total
radiations and photons, respectively. The effective number of neutrinos
$N_{\rm eff}$ can be rewritten,

$\displaystyle N_{\rm eff}=\left(\frac{z_{0}}{z_{\rm
fin}}\right)^{4}\left(3+\frac{\delta\rho_{\nu_{e}}}{\rho_{\nu_{0}}}+\frac{\delta\rho_{\nu_{\mu}}}{\rho_{\nu_{0}}}+\frac{\delta\rho_{\nu_{\tau}}}{\rho_{\nu_{0}}}\right),$
(55)

where $z_{0}=(11/4)^{1/3}\simeq 1.40102$ is the final value of the
dimensionless photon temperature in the instantaneous decoupling limit and
$\delta\rho_{\nu_{\alpha}}=\rho_{\nu_{\alpha}}-\rho_{\nu_{0}}$.

Without neutrino mixing, we find the final values of $N_{\rm eff}$ are 3.03404
for the case without QED corrections and 3.04430 for the case with those up to
$\mathcal{O}(e^{2})$, which agree with recent previous works [12, 14, 15]. In
addition, our results for the case without neutrino mixing but with QED
corrections up to $\mathcal{O}(e^{3})$ show that the final value of $N_{\rm
eff}$ is slightly modified to 3.04335. Thus, the difference with QED
corrections up to $\mathcal{O}(e^{2})$ and $\mathcal{O}(e^{3})$ (but without
neutrino mixing) is 0.00095 in terms of $N_{\rm eff}$, which is very close to
the value estimated in the instantaneous decoupling limit [21].

In the cases with neutrino mixing, Table. 2 shows that the energy densities of
$\mu,\tau$-type neutrinos increase more while those of electron-type neutrinos
increase less, compared to the cases without neutrino mixing. This
modification leads to the enhancement of the total energy density for
neutrinos with final values of $N_{\rm eff}=3.04486$ with QED corrections up
to $\mathcal{O}(e^{2})$ and $N_{\rm eff}=3.04391$ with QED corrections up to
$\mathcal{O}(e^{3})$. Since the blocking factor for electron neutrinos,
$(1-f_{\nu_{e}}),$ is decreased by neutrino mixing, the annihilation of
electrons and positrons into electron neutrinos increases. Although the
annihilation into the other neutrinos decreases, electron neutrinos contribute
to the neutrino heating most efficiently, and neutrino oscillations enhance
the annihilation of electrons and positrons into neutrinos. From these
processes, we conclude that neutrino oscillations slightly promote neutrino
heating and the difference of $N_{\rm eff}$ is $0.00056$, which agrees with
the results of previous works[17, 33, 20]. In this case, we also find that the
difference of $N_{\rm eff}$ between the cases including QED corrections up to
$\mathcal{O}(e^{2})$ and $\mathcal{O}(e^{3})$ is 0.00095.

Finally, we comment on errors of the results due to the numerical calculations
and the choice of physical parameters. Our numerical calculations converge
very well, which is confirmed by the same values of QED corrections at
$\mathcal{O}(e^{3})$ in the two cases with and without neutrino mixing.
However, the results could be affected by neglecting some QED corrections such
as the second line in Eq. (43), neglecting the off-diagonal terms for the
self-interactions of neutrinos and dependence on the values of mixing
parameters, which would change $N_{\rm eff}$ by less than $0.0005$. Taking
into account these errors, we conclude $N_{\rm eff}=3.044$ for the most
accurate case.

Case | $z_{\rm fin}$ | $N_{\rm eff}$
---|---|---
Instantaneous decoupling | 1.40102 | 3.000
No mixing + No QED | 1.39910 | 3.034
No mixing + QED up to $\mathcal{O}(e^{2})$ | 1.39789 | 3.044
No mixing + QED up to $\mathcal{O}(e^{3})$ | 1.39800 | 3.043
mixing + QED up to $\mathcal{O}(e^{2})$ | 1.39786 | 3.045
mixing + QED up to $\mathcal{O}(e^{3})$ | 1.39797 | 3.044
Table 1: The final values of comoving photon temperature and the effective number of neutrinos for flavor neutrinos in several cases. Case | $\delta\bar{\rho}_{\nu_{e}}(\%)$ | $\delta\bar{\rho}_{\nu_{\mu}}(\%)$ | $\delta\bar{\rho}_{\nu_{\tau}}(\%)$ | $\delta\bar{n}_{\nu_{e}}(\%)$ | $\delta\bar{n}_{\nu_{\mu}}(\%)$ | $\delta\bar{n}_{\nu_{\tau}}(\%)$
---|---|---|---|---|---|---
Instantaneous decoupling | 0 | 0 | 0 | 0 | 0 | 0
No mixing + No QED | 0.949 | 0.397 | 0.397 | 0.583 | 0.240 | 0.240
No mixing + QED up to $\mathcal{O}(e^{2})$ | 0.937 | 0.391 | 0.391 | 0.575 | 0.236 | 0.236
No mixing + QED up to $\mathcal{O}(e^{3})$ | 0.937 | 0.391 | 0.391 | 0.575 | 0.236 | 0.236
mixing + QED up to $\mathcal{O}(e^{2})$ | 0.712 | 0.511 | 0.523 | 0.435 | 0.311 | 0.319
mixing + QED up to $\mathcal{O}(e^{3})$ | 0.712 | 0.511 | 0.523 | 0.436 | 0.312 | 0.319
Table 2: The final values of the distortions of energy densities
$\delta\bar{\rho}_{\nu_{\alpha}}\equiv\delta\rho_{\nu_{\alpha}}/\rho_{\nu_{0}}$
and number densities
$\delta\bar{n}_{\nu_{\alpha}}\equiv(n_{\nu_{\alpha}}-n_{\nu_{0}})/n_{\nu_{0}}$
for flavor neutrinos in several cases. Figure 1: The time evolution of the
comoving photon temperature $z(x)$ as a function of the normalized scale
factor $x=m_{e}a$ in the case both with neutrino mixing and QED finite
temperature corrections up to $\mathcal{O}(e^{3})$ . Figure 2: The time
evolution of the distortions of flavor neutrinos for a fixed momentum ($y=5$)
as a function of the normalized scale factor $x=m_{e}a$ in the case with QED
finite temperature corrections up to $\mathcal{O}(e^{3})$. Upper (lower)
dotted line is for $\nu_{e}\ (\nu_{\mu,\tau})$ without neutrino oscillations.
Inner three lines are for flavor neutrinos with neutrino oscillations. Figure
3: The final distortions of flavor neutrino spectra as a function of the
comoving momentum $y$ in the case with QED finite temperature corrections up
to $\mathcal{O}(e^{3})$. Dotted lines represent those for $\nu_{e}\
(\nu_{\mu,\tau})$ without neutrino oscillations, while solid and dashed lines
represent those for flavor neutrinos with neutrino oscillations.

###  3.2 The mass basis

In this section, we present the results of numerical calculations with a set
of Eqs. (20) and (41). In the mass basis, we have also solved these equations
with and without QED finite temperature corrections up to $\mathcal{O}(e^{2})$
and $\mathcal{O}(e^{3})$. Also in the mass basis, we show only the case of the
normal mass hierarchy with the best-fit values of mixing parameters of
neutrinos since we have checked that the results in the normal and inverted
hierarchies are almost the same.

In Fig. 4, we show that the evolution of the massive neutrino spectra,
$f_{\nu_{i}}/f_{\rm eq}$, for a comoving momentum $(y=5)$ as a function of the
normalized scale factor $x$ with QED corrections up to $\mathcal{O}(e^{3})$.
The final values of $f_{\nu_{i}}[y=5]$ are found to be $0.958\%$ for
$\nu_{1}$, $0.724\%$ for $\nu_{2}$, and $0.522\%$ for $\nu_{3}$ larger than
those in the instantaneous decoupling limit, $f_{eq}[y=5]$. These differences
of distortions arise since each massive neutrino interacts with electrons and
positrons through the different coupling $Z^{L}$ in Eq. (35) and the
refractive effects in the mass basis generate the off-diagonal parts of the
mass matrix for massive neutrinos effectively. In Fig. 5, we show the
asymptotic values of the massive neutrino spectra $f_{\nu_{i}}/f_{\rm eq}$ as
a function of $y$ with QED corrections up to $\mathcal{O}(e^{3})$.

In Tables. 3 and 4, we also show that the final values of the dimensionless
photon temperature $z_{\rm fin}$, the energy densities
$\rho_{\nu_{i}}/\rho_{\nu_{0}}$ and number densities $n_{\nu_{i}}/n_{\nu_{0}}$
of massive neutrinos , and the effective number of neutrinos $N_{\rm eff}$.
For the cases with finite temperature corrections up to $\mathcal{O}(e^{2})$
and $\mathcal{O}(e^{3})$, we find very good agreement for $z_{\rm fin}$ and
$N_{\rm eff}$ both in the mass basis and in the flavor basis (with neutrino
mixing). The final values of $N_{\rm eff}$ are 3.04483 for the case with QED
corrections up to $\mathcal{O}(e^{2})$ and 3.04388 for the case with those up
to $\mathcal{O}(e^{3})$. We also find that the difference of $N_{\rm eff}$
between the cases including QED corrections up to $\mathcal{O}(e^{2})$ and
$\mathcal{O}(e^{3})$ is 0.00095, which is the same in the flavor basis. The
small difference for $N_{\rm eff}$ in both bases may come from the fact that
we neglect the off-diagonal parts for self-interaction processes in the
collision terms of the Boltzmann equations.

Case | $z_{\rm fin}$ | $N_{\rm eff}$
---|---|---
QED up to $\mathcal{O}(e^{2})$ | 1.39786 | 3.045
QED up to $\mathcal{O}(e^{3})$ | 1.39797 | 3.044
Table 3: The final values of comoving photon temperature and the effective number of neutrinos for massive neutrinos in several cases. Case | $\delta\bar{\rho}_{\nu_{1}}(\%)$ | $\delta\bar{\rho}_{\nu_{2}}(\%)$ | $\delta\bar{\rho}_{\nu_{3}}(\%)$ | $\delta\bar{n}_{\nu_{1}}(\%)$ | $\delta\bar{n}_{\nu_{2}}(\%)$ | $\delta\bar{n}_{\nu_{3}}(\%)$
---|---|---|---|---|---|---
QED up to $\mathcal{O}(e^{2})$ | 0.764 | 0.573 | 0.409 | 0.468 | 0.350 | 0.248
QED up to $\mathcal{O}(e^{3})$ | 0.764 | 0.574 | 0.409 | 0.468 | 0.350 | 0.248
Table 4: The final values of the distortions of energy densities
$\delta\bar{\rho}_{\nu_{i}}\equiv\delta\rho_{\nu_{i}}/\rho_{\nu_{0}}$ and
number densities
$\delta\bar{n}_{\nu_{i}}\equiv(n_{\nu_{i}}-n_{\nu_{0}})/n_{\nu_{0}}$ for
massive neutrinos in several cases. Figure 4: The time evolution of the
distortions of massive neutrinos for a fixed momentum ($y=5$) in the case with
QED finite temperature corrections up to $\mathcal{O}(e^{3})$. Figure 5: The
final distortions of massive neutrino spectra as a function of the comoving
momentum $y$ in the case with QED finite temperature corrections up to
$\mathcal{O}(e^{3})$.

###  3.3 Transformation of distributions in the flavor and mass bases

In this section, we derive the relation between the distribution functions in
the flavor and mass bases in ultra-relativistic limit and check the numerical
results given in previous sections. The relation of annihilation operators for
negative-helicity neutrinos between flavor and mass eigenstates is given by,

$\displaystyle a_{\alpha}(\bm{p},t)=\sum_{i=1,2,3}U_{\alpha
i}a_{i}(\bm{p},t),$ (56)

with $\alpha=e,\mu,\tau$. Using the above relation, we describe the density
operators for flavor neutrinos $a_{\beta}^{{\dagger}}a_{\alpha}$ through the
operators for massive neutrinos,

$\displaystyle
a^{{\dagger}}_{\beta}(\bm{p},t)a_{\alpha}(\bm{p}^{\prime},t)=\sum_{i,j=1,2,3}U^{\ast}_{\beta
j}U_{\alpha i}a^{{\dagger}}_{j}(\bm{p},t)a_{i}(\bm{p}^{\prime},t),$ (57)

and we get the distribution functions for flavor neutrinos as the density
matrix for massive neutrinos easily,

$\displaystyle f_{\nu_{\alpha}}(\bm{p},t)=\sum_{i,j=1,2,3}U^{\ast}_{\alpha
j}U_{\alpha i}(\rho_{p})_{ij}.$ (58)

In particular, after the decoupling process of neutrinos, the off-diagonal
parts of neutrino density matrix in the mass basis are expected to vanish
since all interactions involving neutrinos are ineffective in this period and
neutrinos in the mass basis do not oscillate. Then the relation between
distribution functions in the flavor and mass bases after the decoupling of
neutrinos is given by,

$\displaystyle f_{\nu_{\alpha}}(\bm{p},t)=\sum_{i}|U_{\alpha
i}|^{2}f_{\nu_{i}}(\bm{p},t).$ (59)

We have numerically confirmed Eq.(59) and $(\rho_{p}(t))_{ij}\simeq 0\ (i\neq
j)$ after the decoupling of neutrinos. Using Eq. (59) and Tables in the mass
basis, we can estimate $N_{\rm eff}=3.04389$ in the flavor basis, which is
very close to the numerical value of $N_{\rm eff}=3.04391$ in the flavor
basis.

##  4 Conclusions

We have studied the neutrino decoupling process in the early Universe by
solving the kinetic equations for neutrinos numerically. We have calculated
the evolution of the neutrino spectral distortions not only in the flavor
basis but also in the mass basis. The latter approach enables us to easily
reveal the neutrino momentum spectra at the current Universe in future work.
The calculations in both bases are also useful for the cross-check of the
results. In addition, in preparation for precision measurements of the
effective number of neutrinos $N_{\rm eff}$, we have also considered the
effects due to QED finite temperature corrections up to $\mathcal{O}(e^{3})$
on the relevant kinetic equations for the first time.

In both bases, we have solved the momentum-dependent kinetic equations for the
neutrino density matrix, where we have considered the full collision terms for
the processes including neutrinos, electrons and their anti-particles while we
neglect the off-diagonal parts of the collision terms for neutrino self-
interaction processes. We find in both bases that the effective number of
neutrinos is $N_{\rm eff}=3.044$. The effects of neutrino oscillations
increase $N_{\rm eff}$ by about 0.0005 compared to $N_{\rm eff}$ without
neutrino oscillations since neutrino oscillations promote the annihilation of
electron-positron pairs into neutrinos. On the other hand, the impacts of QED
corrections up to $\mathcal{O}(e^{3})$ decrease $N_{\rm eff}$ by about 0.001
compared to $N_{\rm eff}$ with QED corrections up to $\mathcal{O}(e^{2})$. The
estimated error for the value of $N_{\rm eff}$ would be at most $0.0005$,
which mainly stems from neglecting some finite temperature corrections and the
off-diagonal terms for self-interactions of neutrinos, and from the dependence
on the values of mixing parameters. We also find the final values of the
number density and energy density for each neutrino. In particular, these
values in the mass basis may be important for detection processes of relic
neutrinos in the current Universe.

The current constraint on $N_{\rm eff}$ from the Planck data analyses [34] in
$\Lambda$CDM model is $N_{\rm eff}=2.99^{+0.34}_{-0.33}$ at $95\%$ CL [35],
consistent with our prediction, $N_{\rm eff}=3.044$. Upcoming CMB and LSS
experiments are expected to improve neutrino masses and energy density bounds
over the next years (see e.g. [36, 37]) and determine $N_{\rm eff}$ with $1\%$
precision in the near future (see e.g. [38, 39, 40, 41, 42]).

Finally, in ultra-relativistic limit and after the decoupling for neutrinos,
we also find the simple transformation formula between the distribution
functions in the flavor and mass bases, which is confirmed by our numerical
calculation. Using this relation, we can easily switch the distribution
functions in the flavor and mass bases for neutrinos without direct numerical
calculations in the two bases.

## Acknowledgments

We are grateful to Shoichi Yamada for useful discussions. KA is supported by
JSPS Grant-in-Aid for Research Fellows No. 19J14449. KA and MY are supported
in part by JSPS Bilateral Open Partnership Joint Research Projects. MY is
supported in part by JSPS Grant-in-Aid for Scientific Research Numbers
18K18764 and Mitsubishi Foundation.

##  Appendix A Kinetic equations for neutrinos in comoving variables

In this appendix, we write the Boltzmann equations and the energy conservation
law in terms of the comoving variables, $x=m_{e}a,\ y=pa,\ z=T_{\gamma}a$. In
terms of these variables, we can write the Boltzmann equations (6) as in [17],

$\displaystyle\frac{d\rho_{y}(x)}{dx}=m_{Pl}\sqrt{\frac{3}{8\pi\bar{\rho}}}\left\\{-i\frac{x^{2}}{m_{e}^{3}}\left[\left(\frac{M^{2}}{2y}-\frac{8\sqrt{2}G_{F}ym_{e}^{6}}{3m_{W}^{2}x^{6}}\bar{E}\right),\rho_{y}(x)\right]+\frac{m_{e}^{3}}{x^{4}}\bar{C}[\rho_{y}(x)]\right\\}.$
(60)

where $\bar{\rho}$, $\bar{E}$, and $\bar{C}[\rho_{y}(x)]$ are quantities
written in the comoving variables, $x,\ y,\ z$. We can write $\bar{\rho}$ and
$\bar{E}$ as

$\displaystyle\bar{\rho}$
$\displaystyle=\rho\left(\frac{x}{m_{e}}\right)^{4},$ $\displaystyle\bar{E}$
$\displaystyle={\rm
diag}\left(\rho_{ee}\left(\frac{x}{m_{e}}\right)^{4},0,0\right).$ (61)

We also give the diagonal collision term from the self-interaction processes
in the comoving variables $\bar{C}_{S}[\nu_{\alpha}(y_{1})]$, where nine-
dimensional collision integrals are reduced to two-dimensional collision
integrals as in appendix B,

$\displaystyle\bar{C}_{S}[\nu_{\alpha}(y_{1})]$
$\displaystyle=\frac{G_{F}^{2}}{2\pi^{3}y_{1}}\int dy_{2}dy_{3}\
y_{2}y_{3}y_{4}\biggl{[}\Pi_{S}^{1}F(\nu_{\alpha}^{(1)},\nu_{\alpha}^{(2)},\nu_{\alpha}^{(3)},\nu_{\alpha}^{(4)})$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
+\Pi_{S}^{2}F(\nu_{\alpha}^{(1)},\nu_{\beta}^{(2)},\nu_{\alpha}^{(3)},\nu_{\beta}^{(4)})+\Pi_{S}^{3}F(\nu_{\alpha}^{(1)},\nu_{\alpha}^{(2)},\nu_{\beta}^{(3)},\nu_{\beta}^{(4)})$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
+\Pi_{S}^{2}F(\nu_{\alpha}^{(1)},\nu_{\gamma}^{(2)},\nu_{\alpha}^{(3)},\nu_{\gamma}^{(4)})+\Pi_{S}^{3}F(\nu_{\alpha}^{(1)},\nu_{\alpha}^{(2)},\nu_{\gamma}^{(3)},\nu_{\gamma}^{(4)})\biggl{]}.$
(62)

Similarly, the collision terms from the annihilation processes and scattering
processes are

$\displaystyle\bar{C}_{A}$ $\displaystyle=\frac{G_{F}^{2}}{2\pi^{3}y_{1}}\int
dy_{2}dy_{3}\ y_{2}y_{3}\bar{E}_{4}$ $\displaystyle\ \ \ \
\times\biggl{[}\Pi_{A}^{1}F^{LL}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)+\Pi_{A}^{2}F^{RR}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)$
$\displaystyle\ \ \ \ \ \ \ \
+\Pi_{A}^{3}\Bigl{(}F^{RL}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)+F^{LR}_{A}\left(\nu^{(1)},\bar{\nu}^{(2)},e^{(3)},\bar{e}^{(4)}\right)\Bigl{)}\biggl{]},$
(63) $\displaystyle\bar{C}_{SC}$
$\displaystyle=\frac{G_{F}^{2}}{2\pi^{3}y_{1}}\int dy_{2}dy_{3}\
y_{2}y_{3}\bar{E}_{4}$ $\displaystyle\ \ \ \
\times\biggl{[}\Pi_{SC}^{1}\Bigl{(}F^{LL}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)+F^{RR}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)\Bigl{)}$
$\displaystyle\ \ \ \ \ \ \ \
-\Pi_{SC}^{2}\Bigl{(}F^{LR}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)+F^{RL}_{SC}\left(\nu^{(1)},e^{(2)},\nu^{(3)},e^{(4)}\right)\Bigl{)}\biggl{]},$
(64)

where $E_{i}=\sqrt{y_{i}^{2}+x^{2}+\delta\bar{m}_{e}^{2}}$ and
$\delta\bar{m}_{e}$ is the finite temperature correction to the electron mass
up to $\mathcal{O}(e^{2})$ in the comoving variables,

$\displaystyle\delta\bar{m}_{e}^{2}=\frac{2\pi\alpha
z}{3}+\frac{4\alpha}{\pi}\int
dy\frac{y^{2}}{\sqrt{y^{2}+x^{2}}}\frac{1}{\exp(\sqrt{y^{2}+x^{2}}/z)+1}.$
(65)

The functions $\Pi_{S,A,SC}^{1,2,3}$ in Eqs. (62), (63) and (64) take the
following forms,

$\displaystyle\Pi_{S}^{1}$
$\displaystyle=6D_{1}-\frac{4D_{2}(y_{1},y_{4})}{y_{1}y_{4}}-\frac{4D_{2}(y_{2},y_{3})}{y_{2}y_{3}}+\frac{2D_{2}(y_{1},y_{2})}{y_{1}y_{2}}+\frac{2D_{2}(y_{3},y_{4})}{y_{3}y_{4}}+\frac{6D_{3}}{y_{1}y_{2}y_{3}y_{4}},$
$\displaystyle\Pi_{S}^{2}$
$\displaystyle=2D_{1}+\frac{D_{2}(y_{1},y_{2})}{y_{1}y_{2}}+\frac{D_{2}(y_{3},y_{4})}{y_{3}y_{4}}-\frac{D_{2}(y_{1},y_{4})}{y_{1}y_{4}}-\frac{D_{2}(y_{2},y_{3})}{y_{2}y_{3}}+\frac{2D_{3}}{y_{1}y_{2}y_{3}y_{4}},$
$\displaystyle\Pi_{S}^{3}$
$\displaystyle=D_{1}-\frac{D_{2}(y_{2},y_{3})}{y_{2}y_{3}}-\frac{D_{2}(y_{1},y_{4})}{y_{1}y_{4}}+\frac{D_{3}}{y_{1}y_{2}y_{3}y_{4}},$
$\displaystyle\Pi_{A}^{1}$
$\displaystyle=2D_{1}-\frac{2D_{2}(y_{2},y_{3})}{y_{2}\bar{E}_{3}}-\frac{2D_{2}(y_{1},y_{4})}{y_{1}\bar{E}_{4}}+\frac{2D_{3}}{y_{1}y_{2}\bar{E}_{3}\bar{E}_{4}},$
$\displaystyle\Pi_{A}^{2}$
$\displaystyle=2D_{1}-\frac{2D_{2}(y_{2},y_{4})}{y_{2}\bar{E}_{4}}-\frac{2D_{2}(y_{1},y_{3})}{y_{1}\bar{E}_{3}}+\frac{D_{3}}{y_{1}y_{2}\bar{E}_{3}\bar{E}_{4}},$
$\displaystyle\Pi_{A}^{3}$
$\displaystyle=(x^{2}+\delta\bar{m}_{e}^{2})\left(D_{1}+\frac{D_{2}(y_{1},y_{2})}{y_{1}y_{2}}\right)\frac{1}{\bar{E}_{3}\bar{E}_{4}},$
$\displaystyle\Pi_{SC}^{1}$
$\displaystyle=4D_{1}-\frac{2D_{2}(y_{2},y_{3})}{\bar{E}_{2}y_{3}}-\frac{2D_{2}(y_{1},y_{4})}{y_{1}\bar{E}_{4}}+\frac{2D_{2}(y_{3},y_{4})}{y_{3}\bar{E}_{4}}+\frac{2D_{2}(y_{1},y_{2})}{y_{1}\bar{E}_{2}}+\frac{4D_{3}}{y_{1}\bar{E}_{2}y_{3}\bar{E}_{4}},$
$\displaystyle\Pi_{SC}^{2}$
$\displaystyle=2(x^{2}+\delta\bar{m}_{e}^{2})\left(D_{1}-\frac{D_{2}(y_{1},y_{3})}{y_{1}y_{3}}\right)\frac{1}{\bar{E}_{2}\bar{E}_{4}}.$
(66)

The functions of $D_{1,2,3}$ are written as,

$\displaystyle D_{1}$
$\displaystyle=\frac{4}{\pi}\int^{\infty}_{0}\frac{d\lambda}{\lambda^{2}}\sin(\lambda
y_{1})\sin(\lambda y_{2})\sin(\lambda y_{3})\sin(\lambda y_{4}),$
$\displaystyle D_{2}(y_{3},y_{4})$
$\displaystyle=\frac{4y_{3}y_{4}}{\pi}\int^{\infty}_{0}\frac{d\lambda}{\lambda^{2}}\sin(\lambda
y_{1})\sin(\lambda y_{2})\left[\cos(\lambda y_{3})-\frac{\sin(\lambda
y_{3})}{\lambda y_{3}}\right]\left[\cos(\lambda y_{4})-\frac{\sin(\lambda
y_{4})}{\lambda y_{4}}\right],$ $\displaystyle D_{3}$
$\displaystyle=\frac{4y_{1}y_{2}y_{3}y_{4}}{\pi}\int^{\infty}_{0}\frac{d\lambda}{\lambda^{2}}\left[\cos(\lambda
y_{1})-\frac{\sin(\lambda y_{1})}{\lambda y_{1}}\right]\left[\cos(\lambda
y_{2})-\frac{\sin(\lambda y_{2})}{\lambda y_{2}}\right]$ $\displaystyle\ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \times\left[\cos(\lambda
y_{3})-\frac{\sin(\lambda y_{3})}{\lambda y_{3}}\right]\left[\cos(\lambda
y_{4})-\frac{\sin(\lambda y_{4})}{\lambda y_{4}}\right],$ (67)

which can be integrated analytically as in appendix B.

Finally, the energy conservation law (20) is translated into the evolution
equation for $z$, including the finite temperature corrections from QED up to
$\mathcal{O}(e^{3})$ [11, 21],

$\displaystyle\frac{dz}{dx}=\frac{\frac{x}{z}J(x/z)-\frac{1}{2\pi^{2}z^{3}}\int^{\infty}_{0}dy\
y^{3}\left(\frac{df_{\nu_{e}}}{dx}+\frac{df_{\nu_{\mu}}}{dx}+\frac{df_{\nu_{\tau}}}{dx}\right)+G^{(2)}_{1}(x/z)+G^{(3)}_{1}(x/z)}{\frac{x^{2}}{z^{2}}J(x/z)+Y(x/z)+\frac{2\pi^{2}}{15}+G^{(2)}_{2}(x/z)+G^{(3)}_{2}(x/z)},$
(68)

where

$\displaystyle G_{1}^{(2)}(\omega)$
$\displaystyle=2\pi\alpha\left[\frac{1}{\omega}\biggl{(}\frac{K(\omega)}{3}+2K(\omega)^{2}-\frac{J(\omega)}{6}-K(\omega)J(\omega)\right)$
$\displaystyle\ \ \ \ \ \ \ \
+\left(\frac{K^{\prime}(\omega)}{6}-K(\omega)K^{\prime}(\omega)+\frac{J^{\prime}(\omega)}{6}+J^{\prime}(\omega)K(\omega)+J(\omega)K^{\prime}(\omega)\right)\biggl{]},$
$\displaystyle G_{2}^{(2)}(\omega)$
$\displaystyle=-8\pi\alpha\left(\frac{K(\omega)}{6}+\frac{J(\omega)}{6}-\frac{1}{2}K(\omega)^{2}+K(\omega)J(\omega)\right)$
$\displaystyle\ \ \ \ \ \ \ \
+2\pi\alpha\omega\left(\frac{K^{\prime}(\omega)}{6}-K(\omega)K^{\prime}(\omega)+\frac{J^{\prime}(\omega)}{6}+J^{\prime}(\omega)K(\omega)+J(\omega)K^{\prime}(\omega)\right)$
$\displaystyle G_{1}^{(3)}(\omega)$
$\displaystyle=\frac{e^{3}}{4\pi}\left(K(\omega)+\frac{\omega^{2}}{2}k(\omega)\right)^{1/2}\biggl{[}\frac{1}{\omega}\left(2J(\omega)-4K(\omega)\right)-2J^{\prime}(\omega)-\omega^{2}j^{\prime}(\omega)$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
-\omega\left(2k(\omega)+j(\omega)\right)-\frac{\left(2J(\omega)+\omega^{2}j(\omega)\right)\left(\omega\left(k(\omega)-j(\omega)\right)+K^{\prime}(\omega)\right)}{2\left(2K+\omega^{2}k(\omega)\right)}\biggl{]}$
$\displaystyle G_{2}^{(3)}(\omega)$
$\displaystyle=\frac{e^{3}}{4\pi}\left(K(\omega)+\frac{\omega^{2}}{2}k(\omega)\right)^{1/2}\left[\frac{(2J(\omega)+\omega^{2}j(\omega))^{2}}{2(2K(\omega)+\omega^{2}k(\omega))}-\frac{2}{\omega}Y^{\prime}(\omega)-\omega\left(3J^{\prime}(\omega)+\omega^{2}j^{\prime}(\omega)\right)\right]$

with

$\displaystyle K(\omega)$
$\displaystyle=\frac{1}{\pi^{2}}\int^{\infty}_{0}du~{}\frac{u^{2}}{\sqrt{u^{2}+\omega^{2}}}\frac{1}{\exp\left(\sqrt{u^{2}+\omega^{2}}\right)+1},$
$\displaystyle J(\omega)$
$\displaystyle=\frac{1}{\pi^{2}}\int^{\infty}_{0}du~{}u^{2}\frac{\exp\left(\sqrt{u^{2}+\omega^{2}}\right)}{\left(\exp\left(\sqrt{u^{2}+\omega^{2}}\right)+1\right)^{2}},$
$\displaystyle Y(\omega)$
$\displaystyle=\frac{1}{\pi^{2}}\int^{\infty}_{0}du~{}u^{4}\frac{\exp\left(\sqrt{u^{2}+\omega^{2}}\right)}{\left(\exp\left(\sqrt{u^{2}+\omega^{2}}\right)+1\right)^{2}},$
$\displaystyle k(\omega)$
$\displaystyle=\frac{1}{\pi^{2}}\int^{\infty}_{0}du\frac{1}{\sqrt{u^{2}+\omega^{2}}}\frac{1}{\exp\left(\sqrt{u^{2}+\omega^{2}}\right)+1},$
$\displaystyle j(\omega)$
$\displaystyle=\frac{1}{\pi^{2}}\int^{\infty}_{0}du\frac{\exp\left(\sqrt{u^{2}+\omega^{2}}\right)}{\left(\exp\left(\sqrt{u^{2}+\omega^{2}}\right)+1\right)^{2}}.$
(70)

The prime represents the derivative with respect to $\omega$. Note that
$G^{(2)}$ and $G^{(3)}$ indicate the finite temperature corrections at
$\mathcal{O}(e^{2})$ and $\mathcal{O}(e^{3})$ respectively.

##  Appendix B Analytic estimation of the collision integral

In this appendix, we analytically perform seven out of nine integrations in
the collision terms for four-Fermi interaction processes in the isotropic
Universe, following refs. [7, 28]. We consider the general form of the
collision term in this case,

$\displaystyle C_{\rm
coll}=\frac{1}{2E_{1}}\int(2\pi)^{4}\delta^{4}(\sum_{i}p_{i})\mathcal{|M|}^{2}F\left(\rho_{p}\right)\prod_{i=2}^{4}\frac{d^{3}\bm{p}_{i}}{(2\pi)^{3}2E_{i}},$
(71)

where $E_{i}$ is the energy of $i$-th particle. The matrix
$F\left(\rho_{p}\right)$ is a function of neutrino density matrix and
$|\mathcal{M}|^{2}$ is a part of the possible squared amplitudes summed over
spin degrees of freedom of all particles except for the first particle. We use
the following relation:

$\displaystyle\delta^{(3)}(\sum_{i}\bm{p}_{i})=\int
e^{\bm{\lambda}\cdot(\bm{p}_{1}+\bm{p}_{2}-\bm{p}_{3}-\bm{p}_{4})}\frac{d^{3}\bm{\lambda}}{(2\pi)^{3}},$
(72)

and decompose momentum integrations into the radial integration and the angle
integrations,

$\displaystyle
d^{3}\bm{p}_{i}=p_{i}^{2}dp_{i}\sin\theta_{i}d\theta_{i}d\phi_{i}\equiv
p_{i}^{2}dp_{i}d\Omega_{i}.$ (73)

Using Eqs. (72) and (73), we write the general collision term (71) as

$\displaystyle C_{\rm
coll}=\frac{1}{64\pi^{3}E_{1}p_{1}}\int\delta(E_{1}+E_{2}-E_{3}-E_{4})F(\rho_{p}(t))D(p_{1},p_{2},p_{3},p_{4})\frac{p_{2}dp_{2}}{E_{2}}\frac{p_{3}dp_{3}}{E_{3}}\frac{p_{4}dp_{4}}{E_{4}},$
(74)

where

$\displaystyle D(p_{1},p_{2},p_{3},p_{4})$
$\displaystyle=\frac{p_{1}p_{2}p_{3}p_{4}}{64\pi^{5}}\int^{\infty}_{0}\lambda^{2}d\lambda\int
e^{i\bm{\lambda}\cdot\bm{p}_{1}}d\Omega_{\lambda}\int
e^{i\bm{\lambda}\cdot\bm{p}_{2}}d\Omega_{p_{2}}$ $\displaystyle\ \ \ \ \ \ \ \
\ \ \ \times\int e^{-i\bm{\lambda}\cdot\bm{p}_{3}}d\Omega_{p_{3}}\int
e^{-i\bm{\lambda}\cdot\bm{p}_{4}}d\Omega_{p_{4}}|\mathcal{M}|^{2}.$ (75)

In the cases of four-Fermi interaction processes, all of $|\mathcal{M}|^{2}$
have two kinds of forms,

$\displaystyle K_{1}({q_{1}}_{\mu}q_{2}^{\mu})({q_{3}}_{\nu}q_{4}^{\nu})$
$\displaystyle=K_{1}(E_{1}E_{2}-\bm{q}_{1}\cdot\bm{q}_{2})(E_{3}E_{4}-\bm{q}_{3}\cdot\bm{q}_{4}),$
(76) $\displaystyle K_{2}m^{2}({q_{3}}_{\mu}q_{4}^{\mu})$
$\displaystyle=K_{2}m^{2}(E_{3}E_{4}-\bm{q}_{3}\cdot\bm{q}_{4}),$ (77)

where $q_{i}$ corresponds to one of $p_{j}$ and the angle between $\bm{q}_{i}$
and $\bm{q}_{j}$ is written in terms of the integration variables of angle,

$\displaystyle\cos\psi_{ij}=\sin\theta_{i}\sin\theta_{j}\cos(\phi_{i}-\phi_{j})+\cos\theta_{i}\cos\theta_{j}.$
(78)

In both cases of Eqs. (76) and (77), we can perform all angle integrals in Eq.
(75) so that we can write $D(q_{1},q_{2},q_{3},q_{4})$ in the case of Eq. (76)
as

$\displaystyle
D=K_{1}[E_{1}E_{2}E_{3}E_{4}D_{1}+E_{1}E_{2}D_{2}(q_{3},q_{4})+E_{3}E_{4}D_{2}(q_{1},q_{2})+D_{3}],$
(79)

while in the case of Eq. (77), $D(q_{1},q_{2},q_{3},q_{4})$ is expressed as

$\displaystyle D=K_{2}E_{1}E_{2}[E_{3}E_{4}D_{1}+D_{2}(q_{3},q_{4})].$ (80)

Here $D_{1,2,3}$ are defined in Eq. (67) and hereafter we only consider
$D_{1},D_{2}(q_{3},q_{4}),D_{3}$.

Although we can perform the integrals in $D_{1,2,3}$ and get the exact
expressions given in ref. [28], we assume for simplicity that $q_{1}>q_{2}$
and $q_{3}>q_{4}$ without loss of generality. Then we get the simplified
expressions of $D_{1,2,3}$ in four different cases:

$(1)\ q_{1}+q_{2}>q_{3}+q_{4}$, $q_{1}+q_{4}>q_{2}+q_{3}$ and $q_{1}\leq
q_{2}+q_{3}+q_{4}$

$\displaystyle D_{1}$ $\displaystyle=\frac{1}{2}(q_{2}+q_{3}+q_{4}-q_{1}),$
$\displaystyle D_{2}(q_{3},q_{4})$
$\displaystyle=\frac{1}{12}\left((q_{1}-q_{2})^{3}+2(q_{3}^{3}+q_{4}^{3})-3(q_{1}-q_{2})(q_{3}^{2}+q_{4}^{2})\right),$
$\displaystyle D_{3}$
$\displaystyle=\frac{1}{60}\bigl{(}q_{1}^{5}-5q_{1}^{3}q_{2}^{2}+5q_{1}^{2}q_{2}^{3}-q_{2}^{5}$
$\displaystyle-5q_{1}^{3}q_{3}^{2}+5q_{2}^{3}q_{3}^{2}+5q_{1}^{2}q_{3}^{3}+5q_{2}^{2}q_{3}^{3}-q_{3}^{5}$
$\displaystyle-5q_{1}^{3}q_{4}^{2}+5q_{2}^{3}q_{4}^{2}+5q_{3}^{3}q_{4}^{2}+5q_{1}^{2}q_{4}^{3}+5q_{2}^{2}q_{4}^{3}+5q_{3}^{2}q_{4}^{3}-q_{4}^{5}\bigl{)}.$
(81)

Note that the case $q_{1}>q_{2}+q_{3}+q_{4}$ is unphysical so that
$D_{1}=D_{2}=D_{3}=0$ in this case.

$(2)\ q_{1}+q_{2}>q_{3}+q_{4}$ and $q_{1}+q_{4}<q_{2}+q_{3}$

$\displaystyle D_{1}$ $\displaystyle=q_{4},$ $\displaystyle
D_{2}(q_{3},q_{4})$ $\displaystyle=\frac{1}{3}q_{4}^{3},$ $\displaystyle
D_{3}$
$\displaystyle=\frac{1}{30}q_{4}^{3}\left(5q_{1}^{2}+5q_{2}^{2}+5q_{3}^{2}-q_{4}^{2}\right).$
(82)

$(3)\ q_{1}+q_{2}<q_{3}+q_{4}$, $q_{1}+q_{4}<q_{2}+q_{3}$ and $q_{3}\leq
q_{1}+q_{2}+q_{4}$

$\displaystyle D_{1}$ $\displaystyle=\frac{1}{2}(q_{1}+q_{2}+q_{4}-q_{3}),$
$\displaystyle D_{2}(q_{3},q_{4})$
$\displaystyle=\frac{1}{12}\left(-(q_{1}+q_{2})^{3}-2q_{3}^{3}+2q_{4}^{3}+3(q_{1}+q_{2})(q_{3}^{3}+q_{4}^{3})\right).$
(83)

$D_{3}$ is equal to that in Eq. (81) with the replacement of variables
$q_{1}\leftrightarrow q_{3}$ and $q_{2}\leftrightarrow q_{4}$ and the case
$q_{3}>q_{1}+q_{2}+q_{4}$ is unphysical so that $D_{1}=D_{2}=D_{3}=0$ in this
case.

$(4)\ q_{1}+q_{2}<q_{3}+q_{4}$ and $q_{1}+q_{4}>q_{2}+q_{3}$

$\displaystyle D_{1}$ $\displaystyle=q_{2},$ $\displaystyle
D_{2}(q_{3},q_{4})$
$\displaystyle=\frac{1}{6}q_{2}\left(3q_{3}^{2}+3q_{4}^{2}-3q_{1}^{2}-q_{2}^{2}\right),$
$\displaystyle D_{3}$
$\displaystyle=\frac{1}{30}q_{2}^{3}\left(5q_{1}^{2}+5q_{3}^{2}+5q_{4}^{2}-q_{2}^{2}\right).$
(84)

After we have integrated the $\delta$-function, we get the simplified
expression of the collision term, leaving two integrals,

$\displaystyle C_{\rm coll}=\frac{1}{64\pi^{3}E_{1}p_{1}}\int\int
F\left(\rho_{p}(t)\right)D(p_{1},p_{2},p_{3},p_{4})\frac{p_{2}dp_{2}}{E_{2}}\frac{p_{3}dp_{3}}{E_{3}},$
(85)

where $E_{4}=E_{1}+E_{2}-E_{3}$ and $p_{4}=\sqrt{E_{4}^{2}-m_{4}^{2}}$.

## References

  * [1] D. A. Dicus, E. W. Kolb, A. M. Gleeson, E. C. G. Sudarshan, V. L. Teplitz and M. S. Turner, Phys. Rev. D 26, 2694 (1982). doi:10.1103/PhysRevD.26.2694
  * [2] A. D. Dolgov, Phys. Rept. 370, 333 (2002) doi:10.1016/S0370-1573(02)00139-4 [hep-ph/0202122].
  * [3] S. Dodelson and M. S. Turner, Phys. Rev. D 46, 3372 (1992). doi:10.1103/PhysRevD.46.3372
  * [4] A. D. Dolgov and M. Fukugita, Phys. Rev. D 46, 5378 (1992). doi:10.1103/PhysRevD.46.5378
  * [5] B. D. Fields, S. Dodelson and M. S. Turner, Phys. Rev. D 47, 4309 (1993) doi:10.1103/PhysRevD.47.4309 [astro-ph/9210007].
  * [6] S. Hannestad and J. Madsen, Phys. Rev. D 52, 1764 (1995) doi:10.1103/PhysRevD.52.1764 [astro-ph/9506015].
  * [7] A. D. Dolgov, S. H. Hansen and D. V. Semikoz, Nucl. Phys. B 503, 426 (1997) doi:10.1016/S0550-3213(97)00479-3 [hep-ph/9703315].
  * [8] A. D. Dolgov, S. H. Hansen and D. V. Semikoz, Nucl. Phys. B 543, 269 (1999) doi:10.1016/S0550-3213(98)00818-9 [hep-ph/9805467].
  * [9] S. Esposito, G. Miele, S. Pastor, M. Peloso and O. Pisanti, Nucl. Phys. B 590, 539 (2000) doi:10.1016/S0550-3213(00)00554-X [astro-ph/0005573].
  * [10] N. Fornengo, C. W. Kim and J. Song, Phys. Rev. D 56, 5123 (1997) doi:10.1103/PhysRevD.56.5123 [hep-ph/9702324].
  * [11] G. Mangano, G. Miele, S. Pastor and M. Peloso, Phys. Lett. B 534, 8 (2002) doi:10.1016/S0370-2693(02)01622-2 [astro-ph/0111408].
  * [12] J. Birrell, C. T. Yang and J. Rafelski, Nucl. Phys. B 890, 481 (2014) doi:10.1016/j.nuclphysb.2014.11.020 [arXiv:1406.1759 [nucl-th]].
  * [13] E. Grohs, G. M. Fuller, C. T. Kishimoto, M. W. Paris and A. Vlasenko, Phys. Rev. D 93, no. 8, 083522 (2016) doi:10.1103/PhysRevD.93.083522 [arXiv:1512.02205 [astro-ph.CO]].
  * [14] E. Grohs and G. M. Fuller, Nucl. Phys. B 923, 222 (2017) doi:10.1016/j.nuclphysb.2017.07.019 [arXiv:1706.03391 [astro-ph.CO]].
  * [15] J. Froustey and C. Pitrou, Phys. Rev. D 101, no. 4, 043524 (2020) doi:10.1103/PhysRevD.101.043524 [arXiv:1912.09378 [astro-ph.CO]].
  * [16] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti and P. D. Serpico, Nucl. Phys. B 729, 221 (2005) doi:10.1016/j.nuclphysb.2005.09.041 [hep-ph/0506164].
  * [17] P. F. de Salas and S. Pastor, JCAP 1607, 051 (2016) doi:10.1088/1475-7516/2016/07/051 [arXiv:1606.06986 [hep-ph]].
  * [18] S. Gariazzo, P. F. de Salas and S. Pastor, JCAP 1907, 014 (2019) doi:10.1088/1475-7516/2019/07/014 [arXiv:1905.11290 [astro-ph.CO]].
  * [19] M. Escudero, JCAP 02, 007 (2019) doi:10.1088/1475-7516/2019/02/007 [arXiv:1812.05605 [hep-ph]].
  * [20] M. Escudero Abenza, JCAP 05, 048 (2020) doi:10.1088/1475-7516/2020/05/048 [arXiv:2001.04466 [hep-ph]].
  * [21] J. J. Bennett, G. Buldgen, M. Drewes and Y. Y. Y. Wong, JCAP 2003, no. 03, 003 (2020) doi:10.1088/1475-7516/2020/03/003 [arXiv:1911.04504 [hep-ph]].
  * [22] A. D. Dolgov, S. H. Hansen, S. Pastor, S. T. Petcov, G. G. Raffelt and D. V. Semikoz, Nucl. Phys. B 632, 363 (2002) doi:10.1016/S0550-3213(02)00274-2 [hep-ph/0201287].
  * [23] Y. Y. Wong, Phys. Rev. D 66, 025015 (2002) doi:10.1103/PhysRevD.66.025015 [arXiv:hep-ph/0203180 [hep-ph]].
  * [24] K. N. Abazajian, J. F. Beacom and N. F. Bell, Phys. Rev. D 66, 013008 (2002) doi:10.1103/PhysRevD.66.013008 [arXiv:astro-ph/0203442 [astro-ph]].
  * [25] G. Mangano, G. Miele, S. Pastor, O. Pisanti and S. Sarikas, Phys. Lett. B 708, 1 (2012) doi:10.1016/j.physletb.2012.01.015 [arXiv:1110.4335 [hep-ph]].
  * [26] E. Castorina, U. Franca, M. Lattanzi, J. Lesgourgues, G. Mangano, A. Melchiorri and S. Pastor, Phys. Rev. D 86, 023517 (2012) doi:10.1103/PhysRevD.86.023517 [arXiv:1204.2510 [astro-ph.CO]].
  * [27] G. Sigl and G. Raffelt, Nucl. Phys. B 406, 423 (1993). doi:10.1016/0550-3213(93)90175-O
  * [28] D. N. Blaschke and V. Cirigliano, Phys. Rev. D 94, no. 3, 033009 (2016) doi:10.1103/PhysRevD.94.033009 [arXiv:1605.09383 [hep-ph]].
  * [29] D. Notzold and G. Raffelt, Nucl. Phys. B 307, 924-936 (1988) doi:10.1016/0550-3213(88)90113-7
  * [30] I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, JHEP 1901, 106 (2019) doi:10.1007/JHEP01(2019)106 [arXiv:1811.05487 [hep-ph]].
  * [31] A. F. Heckler, Phys. Rev. D 49, 611 (1994). doi:10.1103/PhysRevD.49.611
  * [32] R. E. Lopez and M. S. Turner, Phys. Rev. D 59, 103502 (1999) doi:10.1103/PhysRevD.59.103502 [astro-ph/9807279].
  * [33] S. Hannestad, Phys. Rev. D 65, 083006 (2002) doi:10.1103/PhysRevD.65.083006 [astro-ph/0111423].
  * [34] Y. Akrami et al. [Planck], [arXiv:1807.06205 [astro-ph.CO]].
  * [35] N. Aghanim et al. [Planck], [arXiv:1807.06209 [astro-ph.CO]].
  * [36] B. Benson et al. [SPT-3G], Proc. SPIE Int. Soc. Opt. Eng. 9153, 91531P (2014) doi:10.1117/12.2057305 [arXiv:1407.2973 [astro-ph.IM]].
  * [37] P. Ade et al. [Simons Observatory], JCAP 02, 056 (2019) doi:10.1088/1475-7516/2019/02/056 [arXiv:1808.07445 [astro-ph.CO]].
  * [38] K. Abazajian et al. [Topical Conveners: K.N. Abazajian, J.E. Carlstrom, A.T. Lee], Astropart. Phys. 63, 66-80 (2015) doi:10.1016/j.astropartphys.2014.05.014 [arXiv:1309.5383 [astro-ph.CO]].
  * [39] E. Di Valentino et al. [CORE], JCAP 04, 017 (2018) doi:10.1088/1475-7516/2018/04/017 [arXiv:1612.00021 [astro-ph.CO]].
  * [40] S. Hanany et al. [NASA PICO], [arXiv:1902.10541 [astro-ph.IM]].
  * [41] N. Sehgal et al. [arXiv:1906.10134 [astro-ph.CO]].
  * [42] K. Abazajian et al. [arXiv:1907.04473 [astro-ph.IM]].