Next Article in Journal
Gravitating Electron Based on Overrotating Kerr-Newman Solution
Next Article in Special Issue
Editorial for the Special Issue “Recent Advances in Neutrino Physics: From Theory to Experiments”
Previous Article in Journal
Physical Properties of Three Eclipsing Binaries of V Crt, WY Cnc and CG Cyg with Radio Radiation
Previous Article in Special Issue
Shore Shadow Effect in Baikal
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

A Review of Neutrino Decoupling from the Early Universe to the Current Universe

by
Kensuke Akita
1,* and
Masahide Yamaguchi
2
1
Center for Theoretical Physics of the Universe, Institute for Basic Science, Daejeon 34126, Korea
2
Department of Physics, Tokyo Institute of Technology, Tokyo 152-8551, Japan
*
Author to whom correspondence should be addressed.
Universe 2022, 8(11), 552; https://doi.org/10.3390/universe8110552
Submission received: 6 September 2022 / Revised: 19 October 2022 / Accepted: 20 October 2022 / Published: 25 October 2022
(This article belongs to the Special Issue Recent Advances in Neutrino Physics: From Theory to Experiments)

Abstract

:
We review the distortions of spectra of relic neutrinos due to the interactions with electrons, positrons, and neutrinos in the early universe. We solve integro-differential kinetic equations for the neutrino density matrix, including vacuum three-flavor neutrino oscillations, oscillations in electron and positron background, a collision term and finite temperature corrections to electron mass and electromagnetic plasma up to the next-to-leading order O ( e 3 ) . After that, we estimate the effects of the spectral distortions in neutrino decoupling on the number density and energy density of the Cosmic Neutrino Background (C ν B) in the current universe, and discuss the implications of these effects on the capture rates in direct detection of the C ν B on tritium, with emphasis on the PTOLEMY-type experiment. In addition, we find a precise value of the effective number of neutrinos, N eff = 3.044 . However, QED corrections to weak interaction rates at order O ( e 2 G F 2 ) and forward scattering of neutrinos via their self-interactions have not been precisely taken into account in the whole literature so far. Recent studies suggest that these neglections might induce uncertainties of ± ( 10 3 10 4 ) in N eff .

1. Introduction

The successful hot big bang model after inflation predicts that neutrinos produced in the early universe still exist in the current universe. After the temperature of the universe dropped below T 2 MeV , weak interactions became ineffective and neutrinos would have decoupled from thermal plasma. Analogous to photons that make up the Cosmic Microwave Background (CMB), these decoupled neutrinos are called the Cosmic Neutrino Background (C ν B). The existence of these relic neutrinos is confirmed indirectly by the observations of primordial abundances of light elements from the Big Bang Nucleosynthesis (BBN), the anisotropies of the CMB and the distribution of Large Scale Structure (LSS) of the universe. In particular, observations from the Planck satellite impose the severe constraint on the effective number of relativistic species N eff , which describes the total neutrino energy in the Standard Model (SM), and the sum of the neutrino masses at 95 % CL as [1]
N eff 8 7 11 4 4 / 3 ρ r ρ γ 1 = 2 . 99 0.33 + 0.34 and m ν < 0.12 eV ,
where ρ γ and ρ r are the energy densities of photons and radiation, which is composed of photons and neutrinos in the SM, respectively.
Future observations of the C ν B will be developed both indirectly and directly. In fact, CMB-S4 observations are expected to determine N eff with a very good precision of ∼0.03 at 68% C.L. [2]. Thus, an estimation of N eff in the SM with 10 3 precision will be important towards the future CMB-S4 observation. In addition, although it is still very difficult to observe the C ν B in a direct way at present, it is inconceivable that the C ν B will never be directly observed. Among the various discussions on the direct observations, the most promising method of direct detection of the C ν B is neutrino capture on β -decaying nuclei [3,4], ν + n p + e , where there is no threshold energy for relic cosmic neutrinos. In both cases, the theoretical prediction of the relic neutrino spectrum is a crucial ingredient since the radiation energy density in N eff and the direct detection rates depend on the spectrum, and their deviations from the SM suggest physics beyond the SM.
Soon after the decoupling of neutrinos, e ± -pairs start to annihilate and heat photons when the temperature of the universe is T m e = 0.511 MeV . If neutrinos decoupled instantaneously and all electrons and positrons annihilated into photons, the ratio for the temperatures of cosmic photons and neutrinos would be T γ / T ν = ( 11 / 4 ) 1 / 3 1.40102 , due to entropy conservation of the universe. However, the temperatures of neutrino decoupling and e ± -annihilations are so close that e ± -pairs slightly annihilate into neutrinos, which leads to non-thermal distortions in neutrino spectra and a less increase in the photon temperature. These modifications are also parametrized by an increase in N eff from 3.
The non-thermal distortions of relic neutrino spectra and the precise value of N eff have long been studied by solving kinetic equations for neutrinos, which are the Boltzmann equations and the continuity equation. First, several studies solved the Boltzmann equations for neutrino distribution functions [5,6,7,8,9,10,11,12]. Then the kinetic equations were solved with including finite temperature radiative corrections at leading order O ( e 2 ) [13,14,15,16,17,18], and then including three-flavor neutrino oscillations the Boltzmann equations for a neutrino density matrix formalism were solved [19,20,21]. A fast and precise method to calculate effective neutrino temperature for all neutrino species and N eff was also proposed [22,23]. Recently, the authors in ref. [24] pointed out that the finite temperature corrections to electromagnetic plasma at the next-to-leading order O ( e 3 ) are expected to decrease N eff by 10 3 . After that, the present authors found a precise value of N eff = 3.0439 3.044 [25] by solving the Boltzmann equations for the neutrino density matrix including the corrections to electron mass and electromagnetic plasma up to O ( e 3 ) but neglecting off-diagonal parts derived from self-interactions of neutrinos. Later, the authors in refs. [26,27] estimate N eff = 3.0440 and 3.0440 ± 0.0002 , respectively, including off-diagonal parts of the collision term derived by neutrino self-interactions. However, QED corrections to weak interaction rates at the order O ( e 2 G F 2 ) and forward scattering of neutrinos via their self-interactions have not been precisely taken into account in the above references so far. Recent studies [23,28] suggest that these omissions might still induce uncertainties of ± ( 10 3 10 4 ) in N eff .
If we observe the C ν B in a direct way in addition to its indirect observations, we might see neutrino decoupling directly. In the current universe, since the average momentum of the C ν B is p ν 0.53 meV Δ m 21 2 , | Δ m 31 2 | , two massive neutrinos at least are non-relativistic. Under such a situation, it is quite nontrivial to quantize neutrinos in the flavor basis. To reveal the contribution of e ± -annihilation in neutrino decoupling to the spectrum of the C ν B, we calculated the spectra, number densities and energy densities for relic neutrinos in the mass-diagonal basis in the current homogeneous and isotropic universe [25,29].
In this article, we present a review of the distorted spectra of relic cosmic neutrinos from neutrino decoupling to the current universe based on refs. [25,29]. First, in Section 2, we describe the kinetic equations for cosmic neutrinos. In Section 3, we present our results of relic neutrino spectra and N eff . Here we also discuss the uncertainties in N eff . In Section 4, we calculate the number density and energy density of the C ν B in the present universe. In Section 5, the impact of the distortions of the spectra in neutrino decoupling on neutrino capture experiments is also discussed. One of such experiments, which is called the PTOLEMY-type experiment [30,31], uses 100 g of tritium [29,32,33,34,35] as a target through the reaction, ν i + 3 H e + 3 He . Tritium is an appropriate candidate for the target due to its availability, high neutrino capture cross section, low Q-value and long half lifetime of t 1 / 2 = 12.32 years. Here we also include the effects of gravitational clustering of the C ν B by our Galaxy and nearby galaxies based on the results in ref. [36]. Finally, conclusions and discussion are given in Section 6.

2. Kinetic Equations for Neutrinos in Their Decoupling

To follow relic neutrino spectra from neutrino decoupling to the current homogeneous and isotropic universe, we first discuss the field operators and the density matrix for relativistic and non-relativistic neutrinos. Then we introduce the kinetic equations for neutrinos, which are the Boltzmann equations for the evolution of the neutrino density matrix known as the quantum kinetic equations. The continuity equations for the evolution of the total energy density are also introduced.

2.1. Field Operators and Density Matrix

We consider field operators of neutrinos and their density matrices in a homogeneous and isotropic system. With neutrino masses, we cannot define annihilation and creation operators for neutrinos in flavor basis due to their off-diagonal masses in the conventional way, where we interpret these operators as operators that annihilate and create a state with eigenvalues of energy and momentum. On the other hand, in the mass-diagonal basis, we can define such annihilation and creation operators, including neutrino masses. We also compare relic cosmic neutrino spectra obtained in the two bases and confirm their match.
In the ultra-relativistic limit, the field operators for left-handed flavor neutrinos in terms of 4-component spinors, which are composed of only active states for Majorana neutrinos and both active and sterile states for Dirac neutrinos, are expanded in terms of plane wave solutions as
ν α ( x ) = d 3 p ( 2 π ) 3 2 p 0 a α ( p , t ) u p e i p · x + b α ( p , t ) v p e i p · x ,
where a α ( p , t ) = e i H t a α ( p ) e i H t and b α ( p , t ) = e i H t b α ( p ) e i H t are annihilation operators for negative-helicity neutrinos and positive-helicity anti-neutrinos, respectively, and H is the Hamiltonian. α and p are a flavor index and a three dimensional momentum with p 0 | p | , respectively. u p ( v p ) denotes the Dirac spinor for a massless negative-helicity particle (positive-helicity anti-particle), which is normalized to be u p u p = v p v p = 2 p 0 . The annihilation and creation operators satisfy the anti-commutation relations,
{ a α ( p ) , a β ( p ) } = { b α ( p ) , b β ( p ) } = δ α β ( 2 π ) 3 δ ( 3 ) ( p p ) .
For freely evolving massless neutrinos without any interactions, a α 0 ( p , t ) = a α ( p ) e i p 0 t and b α 0 ( p , t ) = b α ( p ) e i p 0 t and the Dirac spinors satisfy free Dirac equations, p u p 0 = 0 , p v p 0 = 0 . On the other hand, for free massive neutrinos in the flavor basis, a α 0 ( p , t ) and b α 0 ( p , t ) cannot be expanded in terms of a plane wave with an eigenvalue of their energy due to off-diagonal neutrino masses. Then we cannot interpret a α ( p , t ) and b α ( p , t ) as annihilation operators except in the ultra-relativistic case.
The density matrices for neutrinos and anti-neutrinos in the flavor basis are defined through the following expectation values of these operators concerning the initial states,
a β ( p , t ) a α ( p , t ) = ( 2 π ) 3 δ ( 3 ) ( p p ) ρ p α β b α ( p , t ) b β ( p , t ) = ( 2 π ) 3 δ ( 3 ) ( p p ) ρ ¯ p α β
where p = | p | . Due to the reversed order of flavor indices in ρ ¯ 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 represent non-zero in the presence of flavor mixing.
On the other hand, in the mass-diagonal basis, the field operators for the negative helicity neutrinos 1 can be expanded as, including neutrino masses,
ν i ( x ) = d 3 p ( 2 π ) 3 2 E i a i ( p , t ) u p ( i ) e i p · x + b i ( p , t ) v p ( i ) e i p · x ,
where i ( = 1 , 2 , 3 ) denotes a mass eigenstate, a i ( p , t ) = e i H t a i ( p ) e i H t , b i ( p , t ) = e i H t b i ( p ) e i H t , E i = p 2 + m i 2 and m i is the neutrino mass in the mass basis. u p ( i ) ( v p ( i ) ) denotes the Dirac spinor for negative-helicity particles (positive-helicity anti-particles), which is also normalized to be u p ( i ) u p ( i ) = v p ( i ) v p ( i ) = 2 E i . For freely evolving neutrinos, a i 0 ( p , t ) = a i ( p ) e i E i t , b i 0 ( p , t ) = b i ( p ) e i E i t and the Dirac spinors satisfy ( p m i ) u p ( i ) , 0 = 0 and ( p + m i ) v p ( i ) , 0 = 0 . As in the flavor basis, the commutation relations for a i ( p ) and b i ( p ) , and the density matrix are defined in the same way except for the exchange of the subscripts, α i .
The diagonalization of the mass matrix for left-handed neutrinos in the flavor basis is achieved through the transformations,
ν α ( x ) = i = 1 , 2 , 3 U α i ν i ( x ) ,
where U α i represents a component of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix U PMNS . Due to Equation (6), in the ultra-relativistic limit, the relation of the density matrices in the flavor and the mass bases is described as
ρ p α β = i , j = 1 , 2 , 3 U β j U α i ρ p i j
In addition, after neutrino decoupling, the off-diagonal parts of the density matrix in the mass basis are zero, ( ρ p ) i j 0 ( i j ) , since all neutrino interactions are ineffective and the oscillations do not occur after neutrino decoupling. In this case, the relations of distribution function in the two bases are simply 2
f ν α ( p , t ) = i = 1 , 2 , 3 | U α i | 2 f ν i ( p , t ) .
Note that Equation (8) is only valid when neutrinos are relativistic and decoupled with thermal plasma. Our numerical calculations also confirm Equation (8).

2.2. Boltzmann Equations

In this section, we derive the Boltzmann equations for the neutrino density matrix, known as quantum kinetic equations, including neutrino oscillations in vacuum, forward scattering with e ± , ν , ν ¯ -background, corresponding to neutrino oscillations in matter, and the collision process at tree level. The resulting Boltzmann equations for neutrinos are summarized in Section 2.5, where we will also discuss the approximations we used in our numerical calculations.

2.2.1. Boltzmann Equations in a Homogeneous and Isotropic System

The Boltzmann equations for neutrinos, including flavor conversion effects, are derived from the Heisenberg equations for the neutrino density operator,
d d t N α β ( t ) = i [ H , N α β ] ,
where [ · , · ] represents the commutator of matrices with a flavor (or mass) index and N α β is the neutrino density operator,
N α β = a β ( p , t ) a α ( p , t ) .
H is the full Hamiltonian in a system, which can be separated into
H = H free + H int ,
where H free is the free Hamiltonian and H int is the interaction Hamiltonian. We assume interactions are enough small that collisions occur individually. Then any fields can be regarded as free ones except during interactions. When the interaction Hamiltonian can be treated perturbatively, the density operator evolves at the first order of H int ,
N α β ( t ) N α β 0 ( t ) + i t 0 t d t H int 0 ( t ) , N α β 0 ( t ) ,
where t 0 is the initial time and H int 0 is the interaction Hamiltonian as a function of freely evolving fields, which are solutions of free Dirac equations, and N α β 0 ( t ) is the free density operator evolved as
N α β 0 ( t ) = e i H free ( t t 0 ) N α β ( t 0 ) e i H free ( t t 0 )
The first order solution (12) includes only neutrino oscillation in vacuum and forward (momentum conserving) scattering with a medium in the system.
To take into account momentum changing collisions, we consider the evolution equation for the density operator at second order of H int , substituting Equation (12) into Equation (9),
d d t N α β ( t ) i H free 0 ( t ) , N α β 0 ( t ) + i H int 0 ( t ) , N α β 0 ( t ) t 0 t d t H int 0 ( t ) , H int 0 ( t ) , N α β 0 ( t ) ,
and an analogous equation for anti-neutrinos [37], N ¯ α β b α ( p , t ) b β ( p , t ) , which is not solved in this article since we assume no lepton asymmetry. Here H free 0 is also the free Hamiltonian as a function of freely evolving fields, where we neglect would-be tiny corrections in the presence of interactions. We also ignore the tiny modification of oscillation and forward scattering, H free , H int , N α β 0 compared with H free 0 , N α β 0 ( t ) and H int 0 , N α β 0 ( t ) . Note that the differential Equation (14) is not closed for both N α β and N α β 0 .
To close and simplify the differential Equation (14), we impose additional approximations. We may set t 0 = 0 and t in the integral range since the time step of the change of N α β , t, may be chosen to be small enough compared to the timescale of the evolution of the universe and large enough compared to the timescale of one collision, t . In addition, at t = 0 , the free density operator coincides with the full one, N α β 0 ( 0 ) = N α β ( 0 ) . Then Equation (14) can be rewritten as
d d t N α β 0 ( 0 ) = i H free 0 ( 0 ) , N α β 0 ( 0 ) + i H int 0 ( 0 ) , N α β 0 ( 0 ) 1 2 d t H int 0 ( 0 ) , H int 0 ( t ) , N α β 0 ( 0 ) .
Thus, the time evolution of the expectation value of N α β 0 ( 0 ) concerning the initial state, ρ p ( 0 ) , is given by
( 2 π ) 3 δ ( 0 ) ( 0 ) d d t ρ p ( 0 ) = i H free 0 ( 0 ) , N α β 0 ( 0 ) + i H int 0 ( 0 ) , N α β 0 ( 0 ) 1 2 d t H int 0 ( 0 ) , H int 0 ( t ) , N α β 0 ( 0 ) .
Equation (16) will be valid at all times, even at t 0 , if in two or more collisions, the correlation of the particles in each collision is independent. This assumption is called molecular chaos in the derivation of the Boltzmann equation. In general, n-point correlation functions are produced by both forward and non-forward collisions. Under the assumption of molecular chaos, n-point correlation functions are reduced to combinations of two-point correlation functions as in ordinary scattering theory. Here two-point correlation functions correspond to distribution functions and neutrino density matrix.
The first term in the right hand side (RHS) represents neutrino oscillations in vacuum and the second term represents forward scattering of neutrinos with background in the system, which is called refractive effects and corresponds to neutrino oscillations in matter. These two terms do not change neutrino momenta but induce flavor conversions. The third term represents scattering and annihilation including both momentum conserving and changing processes, usually rewritten as
1 2 d t H int 0 ( 0 ) , H int 0 ( t ) , N α β 0 ( 0 ) ( 2 π ) 3 δ ( 3 ) ( 0 ) C ρ p ( t ) ,
where C ρ p ( t ) is called the collision term. In the following sections, we calculate the formulae of these three terms. The resulting Boltzmann equations for the neutrino density matrix are summarized in Section 2.5.

2.2.2. Neutrino Oscillation in Vacuum

The calculation of the first term in the RHS of Equation (16) is well established in the mass basis. The free Hamiltonian of neutrinos in the mass basis is given by
H free = d 3 x i = 1 3 ν ¯ i ( i γ · + m i ) ν i ,
where γ = ( γ 1 , γ 2 , γ 3 ) are the gamma matrices. After substituting the free operators for left-handed neutrinos, the free Hamiltonian becomes
H free 0 = d 3 p i = 1 3 a i ( p ) E i a i ( p ) + b i ( p ) E i b i ( p ) .
The first term in the RHS of Equation (16) in the mass basis is written as
i H free 0 , N i j 0 ( 0 ) = i ( 2 π ) 3 δ ( 3 ) ( 0 ) diag ( E 1 , E 2 , E 3 ) , ρ p , i ( 2 π ) 3 δ ( 3 ) ( 0 ) M diag 2 2 p , ρ p ,
where M diag 2 = diag ( m ν 1 2 , m ν 2 2 , m ν 3 2 ) and i , j denote mass-eigenstates. In the flavor basis, as in discussed in Section 2.1, it is quite nontrivial to quantize neutrinos in the flavor basis with non-zero masses. When we calculate the first term in the RHS of Equation (16) in the flavor basis directly, we replace the free annihilation operators a α 0 ( p , t ) and b α 0 ( p , t ) with a α osc ( p , t ) = ( exp ( i Ω p t ) ) α β a β ( p ) and b α osc ( p , t ) = ( exp ( i Ω p t ) ) α β b β ( p ) as in [37], where Ω p = p 2 + M 2 . Then we also obtain the first term of Equation (16), following the similar procedure in the mass basis,
i H free 0 , N α β 0 ( 0 ) i ( 2 π ) 3 δ ( 3 ) ( 0 ) M 2 2 p , ρ p ,
where M 2 = U PMNS M diag 2 U PMNS is the neutrino mass matrix in the flavor basis. For anti-neutrinos, the corresponding term is obtained by adding a minus sign for the reverse indices in the anti-neutrino density matrix (4), i [ H free 0 , N ¯ α β 0 ( 0 ) ] i ( 2 π ) 3 δ ( 3 ) ( 0 ) [ M 2 / 2 p , ρ ¯ p ] .

2.2.3. Forward Scattering with e ± , ν , ν ¯ -Background

In the following of Section 2, we consider the flavor basis of neutrinos. Forward scattering of neutrinos with background in the system called refractive effects modifies neutrino oscillations through the one-loop thermal interaction as given in Figure 1. Since the temperature in thermal plasma is ∼2 MeV in neutrino decoupling, particles except for photons, electrons, neutrinos and their anti-particles are already annihilated due to their heavy masses. Then we consider only e ± , ν , ν ¯ -background. The interaction Hamiltonian is described as
H int = g 2 2 d 3 x d 4 y d 4 p ( 2 π ) 4 e i p ( x y ) D μ ν Z ( p ) J N C μ ( x ) J N C ν ( y ) + 2 D μ ν W ( p ) J C C μ ( x ) J C C ν ( y ) , H N C + H C C ,
where D μ ν Z ( p ) and D μ ν W ( p ) are the full propagator of Z 0 boson and W ± boson,
D μ ν W , Z ( p ) = g μ ν p μ p ν m W , Z 2 1 m W , Z 2 p 2 , g μ ν m W , Z 2 + g μ ν p 2 p μ p ν m W , Z 4 .
Here g , m Z , m W are the electroweak coupling constant, the Z 0 boson mass and the W ± boson mass, respectively. The neutral current and the charged current are given by
J N C μ J ν ν μ + J e e L μ + J e e R μ , J C C μ J e ν e μ ,
where
J ν ν μ = 1 4 cos θ W ν ¯ γ μ ( 1 γ 5 ) ν , J e e L μ = 1 2 cos θ W 1 2 + sin 2 θ W e ¯ γ μ ( 1 γ 5 ) e , J e e R μ = 1 2 cos θ W sin 2 θ W e ¯ γ μ ( 1 + γ 5 ) e , J e ν e μ = 1 2 2 ν ¯ e γ μ ( 1 γ 5 ) e ,
with
ν = ν e ν μ ν τ .
Here θ W is the weak mixing angle, e is the field operator for electron and positron and ν α is the field operator for neutrinos and anti-neutrinos with a flavor α .
The interaction Hamiltonian is divided into the two parts corresponding to the neutral current interaction, H N C J N C μ , and to the charged current interaction, H C C J C C μ . For the charged current interactions, the second term in the RHS of Equation (16), which represents forward scattering of neutrinos with e ± -background, is given by [37,38]
i H C C 0 ( 0 ) , N α β 0 ( 0 ) = i ( 2 π ) 3 δ ( 3 ) ( 0 ) 2 G F ( N e N e + ) 2 2 G F p 3 m Z 2 E e + P e + E e + + P e + , ρ p ,
where G F is the Fermi coupling constant and N e ± , E e ± and P e ± are the number density, energy density and pressure for e ± -background, respectively, which are described in the flavor basis as
N e diag ( n e , 0 , 0 ) , N e + diag ( n e + , 0 , 0 ) , n e ± = 2 d 3 p ( 2 π ) 3 f e ± ( p ) , E e ± + P e ± diag ( ρ e ± + P e ± , 0 , 0 ) , ρ e ± + P e ± = d 3 p ( 2 π ) 3 E e + p 2 3 E e f e ± ( p ) ,
where E e = p 2 + m e 2 . In the temperature of MeV scale in neutrino decoupling, the densities for muons and tauons are enough suppressed by their heavy masses. We neglect forward scattering of neutrinos with muons and tauons, which corresponds to the second and third diagonal components in Equation (28).
For the neutral current interactions, the second term in the RHS of Equation (16), which represents forward scattering of neutrinos via neutrino self-interactions, is given by [37,38]
i H N C 0 ( 0 ) , N α β 0 ( 0 ) = i ( 2 π ) 3 δ ( 3 ) ( 0 ) 2 G F ( N ν N ν ¯ ) 8 2 G F p 3 m Z 2 E ν + E ν ¯ , ρ p ,
where N ν , N ν ¯ , E ν and E ν ¯ are the number and energy densities for the density matrices of ν , ν ¯ -background, respectively, which are described in the flavor basis as
N ν = d 3 p ( 2 π ) 3 ρ p , N ν ¯ = d 3 p ( 2 π ) 3 ρ ¯ p , E ν = d 3 p ( 2 π ) 3 p ρ p , E ν ¯ = d 3 p ( 2 π ) 3 p ρ ¯ p ,
where we neglect neutrino masses since neutrinos are relativistic in neutrino decoupling.
For anti-neutrinos, the corresponding terms, i [ H C C 0 ( 0 ) , N ¯ α β 0 ( 0 ) and i [ H C C 0 ( 0 ) , N ¯ α β 0 ( 0 ) , are obtained by adding an overall minus sign for the reverse indices in the anti-neutrino density matrix (4) and replacing N e N e + ( N e N e + ) and N ν N ν ¯ ( N ν N ν ¯ ) for an opposite evolution of anti-neutrinos due to the lepton asymmetry in Equations (27) and (29) [37].
If there is a large lepton asymmetry, the terms proportional to N e N e + and N ν N ν ¯ will be important. Note that even if there is no lepton asymmetry, the off-diagonal parts of N ν N ν ¯ have non-zero contribution since the density matrices for neutrinos and anti-neutrinos follow the same evolution, ρ p = ρ ¯ p T ρ ¯ p , in the case of no lepton asymmetry.

2.2.4. Collision Term

Finally we discuss the third term in the RHS of Equation (16) called the collision term. The temperature of ∼2 MeV in neutrino decoupling is much lower than the electroweak scale of ∼ m Z , m W . After integrating out Z 0 and W ± bosons in the instantaneous interaction limit, the interaction Hamiltonian in neutrino decoupling can be written as
H int g 2 2 d 3 x [ 1 m Z 2 J N C μ ( x ) J N C μ ( x ) + 2 m W 2 J C C μ ( x ) J C C μ ( x ) ] .
The interaction Hamiltonian can be divided into the part including both neutrinos and electrons (and their anti-particles), and the one only including neutrinos and anti-neutrinos, H int H int e ν + H int ν , while we ignore the part including only electrons and positrons,
H int e ν = G F 2 d x 3 ν ¯ γ μ ( 1 γ 5 ) Y L ν e ¯ γ μ ( 1 γ 5 ) e + ν ¯ γ μ ( 1 γ 5 ) Y R ν e ¯ γ μ ( 1 + γ 5 ) e , H int ν = G F 4 2 d x 3 ν ¯ γ μ ( 1 γ 5 ) ν ν ¯ γ μ ( 1 γ 5 ) ν ,
with
Y L = 1 2 + sin 2 θ W 0 0 0 1 2 + sin 2 θ W 0 0 0 1 2 + sin 2 θ W , Y R = sin 2 θ W × 1 .
Here we have used the following Fierz transformation in the charged currents,
ν ¯ e γ μ ( 1 γ 5 ) e e ¯ γ μ ( 1 γ 5 ) ν e = ν ¯ e γ μ ( 1 γ 5 ) ν e e ¯ γ μ ( 1 γ 5 ) e .
Due to this contribution of the charged current, only electron-type neutrinos and anti-neutrinos interact with electrons and positrons via different magnitudes of interactions, compared to other flavor neutrinos with ( Y L ) 11 = ( Y L ) 22 ( 33 ) + 1 .
The Hamiltonian of Equation (32) can be further divided as
H int e ν = H ν ν ¯ e e + + H ν e ± ν e ± + H ν ¯ e ± ν ¯ e ± , H int ν = H ν ν ν ν + H ν ν ¯ ν ν ¯ ,
where H a b c d is the term including operators of (anti-)particles, a , b , c and d. In the following, we neglect H ν ¯ e ± ν ¯ e ± since this Hamiltonian does not contribute the evolution of neutrinos. In addition, we only consider contributions proportional to the following terms as a function of freely evolving fields in the collision term in Equation (16),
[ H ν ν ¯ e e + 0 , [ H ν ν ¯ e e + 0 , N α β 0 ] , [ H ν e ± ν e ± 0 , [ H ν e ± ν e ± 0 , N α β 0 ] ] , [ H ν ν ν ν 0 , [ H ν ν ν ν 0 , N α β 0 ] ] , [ H ν ν ¯ ν ν ¯ 0 , [ H ν ν ¯ ν ν ¯ 0 , N α β 0 ] .
The other terms also denote forward scattering, which would give tiny modifications of Equations (27) and (29). The first term in Equation (36) denotes the annihilation of neutrinos and anti-neutrinos into e ± -pairs, which mainly contribute to the distortion of neutrino spectrum in their decoupling. The second term denotes the scattering between neutrinos and electrons (positrons). The third term represents the scattering process including only neutrinos while the fourth term denotes the annihilation and scattering processes of neutrinos and anti-neutrinos.
In a schematic manner, the collision term for two-body reactions 1 + 2 3 + 4 at tree level takes the following expressions,
( 2 π ) 3 δ ( 3 ) ( 0 ) C [ ρ p 1 ] = 1 2 d t [ H int 0 ( 0 ) , [ H int 0 , N α β 0 ] = ( 2 π ) 3 δ ( 3 ) ( 0 ) 1 2 E 1 d 3 p 2 ( 2 π ) 3 2 E 2 d 3 p 3 ( 2 π ) 3 2 E 3 d 3 p 4 ( 2 π ) 3 2 E 4 × ( 2 π ) 4 δ ( 4 ) ( p 1 + p 2 p 3 p 4 ) F ( ρ , f e ± , Y L , Y R ) S | M | 12 34 2 part ,
where ρ i ( i = 1 , 2 , 3 , 4 ) denote the neutrino density matrix, not the energy density, and E i | p i | for ν and ν ¯ while E i = p i 2 + m e 2 for e ± . F ( ρ , f e ± , Y L , Y R ) is a matrix depending on ρ , f e ± , Y L and/or Y R . S | M | 12 34 2 part is a part of S | M | 12 34 2 , where S is the symmetric factor and | M | 2 is the squared matrix element summed over spins of all particles except for the first one. The formulae of S | M | 2 for the relevant reaction in neutrino decoupling are shown in Table 1. Nine integrals in the collision term in Equation (37) can be reduced analytically to two integrals as in Appendix B.
In the following, we rewrite the collision terms C [ ρ p ( t ) ] including Equation (36) with neutrino density matrices and the distribution functions of electrons and positrons. The formulae of the collision terms for neutrino density matrix are originally given in refs. [37,39], and for numerical calculations of neutrino spectra, these formulae are developed in refs. [20,26].
(i)
ν ( p 1 ) + ν ¯ ( p 2 ) e ( p 3 ) + e + ( p 4 )
The collision term for the annihilation process including e ± , ν ( p 1 ) + ν ¯ ( p 2 ) e ( p 3 ) + e + ( p 4 ) , comes from the term proportional to [ H ν ν ¯ e e + 0 , [ H ν ν ¯ e e + 0 , N α β 0 ] . We can calculate the corresponding collision terms, which are denoted as ( 2 π ) 3 δ ( 3 ) ( 0 ) C ν ν ¯ e e + [ ρ p 1 ( t ) ] ,
( 2 π ) 3 δ ( 3 ) ( 0 ) C ν ν ¯ e e + [ ρ p 1 ( t ) ] = 1 2 d t [ H ν ν ¯ e e + 0 ( 0 ) , [ H ν ν ¯ e e + 0 ( t ) , N α β 0 ] = ( 2 π ) 3 δ ( 3 ) ( 0 ) 1 2 2 5 G F 2 2 | p 1 | d 3 p 2 ( 2 π ) 3 2 | p 2 | d 3 p 3 ( 2 π ) 3 2 E 3 d 3 p 4 ( 2 π ) 3 2 E 4 ( 2 π ) 4 δ ( 4 ) ( p 1 + p 2 p 3 p 4 ) × [ 4 ( p 1 · p 4 ) ( p 2 · p 3 ) F ann L L ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) + 4 ( p 1 · p 3 ) ( p 2 · p 4 ) F ann R R ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) + 2 ( p 1 · p 2 ) m e 2 ( F ann L R ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) + F ann R L ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) ) ] ,
where
F ann a b ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) = f e ( p 3 ) f e + ( p 4 ) ( Y a ( 1 ρ ¯ 2 ) ) Y b ( 1 ρ 1 ) + ( 1 ρ 1 ) Y b ( 1 ρ ¯ 2 ) Y a ) 1 f e ( p 3 ) 1 f e + ( p 4 ) ( Y a ρ ¯ 1 Y b ρ 1 + ρ 1 Y b ρ ¯ 2 Y a ) .
Here f e ± ( p ) is the distribution function for electrons and positrons, respectively.
(ii)
ν ( p 1 ) + e ± ( p 2 ) ν ( p 3 ) + e ± ( p 4 )
The collision term for the scatterings including e ± , ν ( p 1 ) + e ± ( p 2 ) ν ( p 3 ) + e ± ( p 4 ) , comes from the term proportional to [ H ν e ± ν e ± 0 , [ H ν e ± ν e ± 0 , N α β 0 ] . We can similarly calculate the corresponding collision term, which is denoted as C ν e ν e [ ρ p 1 ( t ) ] and C ν e + ν e + [ ρ p 1 ( t ) ] , respectively,
C ν e ν e [ ρ p 1 ( t ) ] = 1 2 2 5 G F 2 2 p 1 d 3 p 2 ( 2 π ) 3 2 E 2 d 3 p 3 ( 2 π ) 3 2 p 3 d 3 p 4 ( 2 π ) 3 2 E 4 ( 2 π ) 4 δ ( 4 ) ( p 1 + p 2 p 3 p 4 ) × [ 4 ( p 1 · p 2 ) ( p 3 · p 4 ) F sc L L ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) + 4 ( p 1 · p 4 ) ( p 2 · p 3 ) F sc R R ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) 2 ( p 1 · p 3 ) m e 2 ( F sc L R ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) + F sc R L ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) ) ] ,
and
C ν e + ν e + [ ρ p 1 ( t ) ] = 1 2 2 5 G F 2 2 p 1 d 3 p 2 ( 2 π ) 3 2 E 2 d 3 p 3 ( 2 π ) 3 2 p 3 d 3 p 4 ( 2 π ) 3 2 E 4 ( 2 π ) 4 δ ( 4 ) ( p 1 + p 2 p 3 p 4 ) × [ 4 ( p 1 · p 2 ) ( p 3 · p 4 ) F sc R R ν ( 1 ) , e + ( 2 ) , ν ( 3 ) , e + ( 4 ) + 4 ( p 1 · p 4 ) ( p 2 · p 3 ) F sc L L ν ( 1 ) , e + ( 2 ) , ν ( 3 ) , e + ( 4 ) 2 ( p 1 · p 3 ) m e 2 ( F sc L R ν ( 1 ) , e + ( 2 ) , ν ( 3 ) , e + ( 4 ) + F sc R L ν ( 1 ) , e + ( 2 ) , ν ( 3 ) , e + ( 4 ) ) ] ,
where
F sc a b ν ( 1 ) , e ± ( 2 ) , ν ( 3 ) , e ± ( 4 ) = f e ± ( p 4 ) ( 1 f e ± ( p 2 ) ) ( Y a ρ 3 Y b ( 1 ρ 1 ) + ( 1 ρ 1 ) Y b ρ 3 Y a ) f e ± ( p 2 ) ( 1 f e ± ( p 4 ) ) ( ρ 1 Y b ( 1 ρ 3 ) Y a + Y a ( 1 ρ 3 ) Y b ρ 1 ) .
(iii)
ν ( p 1 ) + ν ( p 2 ) ν ( p 3 ) + ν ( p 4 ) and ν ( p 1 ) + ν ¯ ( p 2 ) ν ( p 3 ) + ν ¯ ( p 4 )
The collision terms for the scatterings including only neutrinos and anti-neutrinos, ν ( p 1 ) + ν ( p 2 ) ν ( p 3 ) + ν ( p 4 ) and ν ( p 1 ) + ν ¯ ( p 2 ) ν ( p 3 ) + ν ¯ ( p 4 ) , come from the term proportional to [ H ν ν ν ν 0 , [ H ν ν ν ν 0 , N α β 0 ] ] and [ H ν ν ¯ ν ν ¯ 0 , [ H ν ν ¯ ν ν ¯ 0 , N α β 0 ] ] , respectively. The corresponding collision terms, which are denoted as C ν ν ν ν [ ρ p 1 ( t ) ] and C ν ν ¯ ν ν ¯ [ ρ p 1 ( t ) ] , respectively, are calculated as
C ν ν ν ν [ ρ p 1 ( t ) ] = 1 2 2 5 G F 2 2 p 1 d 3 p 2 ( 2 π ) 3 2 | p 2 | d 3 p 3 ( 2 π ) 3 2 p 3 d 3 p 4 ( 2 π ) 3 2 | p 4 | ( 2 π ) 4 δ ( 4 ) ( p 1 + p 2 p 3 p 4 ) × ( p 1 · p 2 ) ( p 3 · p 4 ) F sc ν ( 1 ) , ν ( 2 ) , ν ( 3 ) , ν ( 4 ) ,
C ν ν ¯ ν ν ¯ [ ρ p 1 ( t ) ] = 1 2 2 5 G F 2 2 p 1 d 3 p 2 ( 2 π ) 3 2 | p 2 | d 3 p 3 ( 2 π ) 3 2 p 3 d 3 p 4 ( 2 π ) 3 2 | p 4 | ( 2 π ) 4 δ ( 4 ) ( p 1 + p 2 p 3 p 4 ) × ( p 1 · p 4 ) ( p 2 · p 3 ) F sc ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ ( 4 ) + F ann ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ ( 4 ) ,
where F sc ν ( 1 ) , ν ( 2 ) , ν ( 3 ) , ν ( 4 ) , F sc ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ ( 4 ) and F ann ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ ( 4 ) denote contributions from scatterings for ν ν ν ν , scatterings and annihilations for ν ν ¯ ν ν ¯ , respectively,
F sc ν ( 1 ) , ν ( 2 ) , ν ( 3 ) , ν ( 4 ) = ρ 4 ( 1 ρ 2 ) + Tr ( . . . ) ρ 3 ( 1 ρ 1 ) + ( 1 ρ 1 ) ρ 3 ( 1 ρ 2 ) ρ 4 + Tr ( . . . ) ( 1 ρ 4 ) ρ 2 + Tr ( . . . ) ( 1 ρ 3 ) ρ 1 ρ 1 ( 1 ρ 3 ) ρ 2 ( 1 ρ 4 ) + Tr ( . . . ) ,
F sc ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ ( 4 ) = ( 1 ρ ¯ 2 ) ρ ¯ 4 + Tr ( . . . ) ρ 3 ( 1 ρ 1 ) + ( 1 ρ 1 ) ρ 3 ρ ¯ 4 ( 1 ρ ¯ 2 ) + Tr ( . . . ) ρ ¯ 2 ( 1 ρ ¯ 4 ) + Tr ( . . . ) ( 1 ρ 3 ) ρ 1 ρ 1 ( 1 ρ 3 ) ( 1 ρ ¯ 4 ) ρ ¯ 2 + Tr ( . . . ) ,
F ann ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ ( 4 ) = [ ρ 3 ρ ¯ 4 + Tr ( . . . ) ] ( 1 ρ ¯ 2 ) ( 1 ρ 1 ) + ( 1 ρ 1 ) ( 1 ρ ¯ 2 ) [ ρ ¯ 4 ρ 3 + Tr ( . . . ) ] [ ( 1 ρ 3 ) ( 1 ρ ¯ 4 ) + Tr ( . . . ) ] ρ ¯ 2 ρ 1 ρ 1 ρ ¯ 2 [ ( 1 ρ ¯ 4 ) ( 1 ρ 3 ) + Tr ( . . . ) ] ,
where [ α + Tr ( . . . ) ] [ α + Tr ( α ) ] .
Finally, we obtain the collision term in Equation (16), C [ ρ p ( t ) ] , combining Equations (38), (40), (41), (43) and (44),
C [ ρ p ( t ) ] = C ν ν ¯ e e + + C ν e ν e + C ν e + ν e + + C ν ν ν ν + C ν ν ¯ ν ν ¯ .
The collision terms for anti-neutrinos can be obtained by appropriately replacing the density matrices and momenta, ρ i ρ ¯ i and p i p j [37,40]. Changing the collision term for ν ( p 1 ) X ν ( p 3 ) X to ν ¯ ( p 1 ) X ν ¯ ( p 3 ) X corresponds replacing ρ 1 ρ ¯ 1 , ρ 3 ρ ¯ 3 and p 1 p 3 in this collision term while changing that for ν ( p 1 ) ν ¯ ( p 2 ) X X to ν ¯ ( p 1 ) ν ( p 2 ) X X corresponds ρ 1 ρ ¯ 1 , ρ ¯ 2 ρ 2 and p 1 p 2 . One may consider the transpose in the collision terms is necessary for the reverse indices in the anti-neutrino density matrix (4), but this is not necessary since the collision terms are invariant under the transpose.

2.3. Continuity Equation

In addition to the Boltzmann equations for the neutrino density matrix, the energy conservation law must be satisfied,
d ρ d t = 3 H ( ρ + P ) ,
where ρ and P are the total energy density and pressure of γ , e ± , ν , ν ¯ around MeV-scale temperature, respectively. The continuity equation corresponds to the evolution of the photon temperature T γ .
Though we will discuss finite temperature corrections from QED to ρ , P and m e in the next section, in the ideal gas limit, they are given as follows, which are denoted by ρ ( 0 ) and P ( 0 ) , respectively,
ρ ( 0 ) = π 2 T γ 4 15 + 2 π 2 d p p 2 p 2 + m e 2 exp ( p 2 + m e 2 / T γ ) + 1 + α = e , μ , τ 1 π 2 d p p 3 f ν α ( p ) , P ( 0 ) = π 2 T γ 4 45 + 2 π 2 d p p 4 3 p 2 + m e 2 [ exp ( p 2 + m e 2 / T γ ) + 1 ] + α = e , μ , τ 1 3 π 2 d p p 3 f ν α ( p ) .
The Hubble parameter in Equation (49) is calculated using the usual relation, 3 H 2 m Pl 2 = 8 π ρ with m 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.4. Finite Temperature QED Corrections to m e , ρ and P up to O ( e 3 )

QED interactions at finite temperature modify the energy density and pressure of electromagnetic plasma from the ideal gas limit. In addition, their interactions change the electron mass (and produce an effective photon mass). These corrections affect the kinetic equations for neutrinos discussed in the former sections. The corrections to the electron mass modify the weak interaction rates and the distribution function for e ± . Through the direct modifications of ρ and P, the expansion rate H is also changed. Note that QED interactions also modify weak interaction rates in the collision term C [ ρ p ( t ) ] and the Hamiltonian for the forward scattering (60) at order O ( e 2 G F 2 ) directly. In our numerical calculations, we consider corrections to weak interaction rates only due to the change of m e . We will discuss other QED corrections to weak interaction rates and their uncertainties in N eff in Section 3.3.1.
The corrections to the grand canonical partition function Z by interactions at finite temperature are well established perturbatively and can be calculated by the similar procedure of the functional integrals of Quantum Field Theory (QFT) at zero temperature after changing t i / T . P and ρ are described by Z as
P = T V ln Z , ρ = T 2 V ln Z T = P + T P T ,
where T and V are the temperature and volume in the system, respectively. Then we can expand ln Z in powers of the QED coupling constant e as ln Z = n = 1 ln Z ( n ) , where ln Z ( n ) e n . In the isotropic and lepton symmetric universe, the corresponding corrections to P and ρ at O ( e 2 ) , P ( 2 ) , ρ ( 2 ) , e 2 , are [41]
P ( 2 ) = e 2 T γ 2 12 π 2 0 d p p 2 E p N F ( p ) e 2 8 π 4 0 d p p 2 E p N F ( p ) 2 + e 2 m e 2 16 π 4 0 0 d p d p p p E p E p ln p + p p p N F ( p ) N F ( p ) , ρ ( 2 ) = P ( 2 ) + T γ P ( 2 ) T γ ,
where E p = p 2 + m e 2 and N F ( p ) is the sum of the distribution functions for e ± ,
N F ( p ) = 2 1 e E p / T γ + 1 .
The next-to-leading order of thermal corrections to ρ , P is O ( e 3 ) , not O ( e 4 ) . These non-trivial corrections come from the resummation of ring diagrams in the photon propagator at all orders. The thermal corrections to P , ρ at O ( e 3 ) , P ( 3 ) , ρ ( 3 ) , e 3 , are [24,41],
P ( 3 ) = e 3 T γ 12 π 4 I 3 / 2 ( T γ ) , ρ ( 3 ) = e 3 T γ 2 8 π 4 I 1 / 2 I T γ ,
where
I ( T γ ) = 0 d p p 2 + E p 2 E p N F ( p ) .
Finally, we read the total energy density and the total pressure of electromagnetic plasma up to O ( e 3 ) corrections as
P = P ( 0 ) + P ( 2 ) + P ( 3 ) , ρ = ρ ( 0 ) + ρ ( 2 ) + ρ ( 3 ) .
The thermal corrections to the e ± mass at O ( e 2 ) is given by, through modifications of the e ± self energy [42],
δ m e ( 2 ) 2 ( p , T γ ) = e 2 T γ 2 6 + e 2 2 π 2 0 d p k 2 E p N F ( p ) e 2 m e 2 4 π 2 p 0 d p p E p log p + p p p N F ( p ) .
The last logarithmic terms in Equations (52) and (57) give less than 10 % corrections to these equations around the decoupling temperature and the average momentum of electrons [43]. These terms also give contributions less than 10 4 to N eff [24,27]. In the following, we neglect the logarithmic corrections. Note that thermal corrections to m e at O ( e 3 ) do not appear because O ( e 3 ) corrections stem from ring diagrams in the photon propagator.

2.5. Summary and Approximations

In this section, we summarize the closed system of the resulting Boltzmann equations for the neutrino density matrix and the continuity equation in neutrino decoupling. We also discuss the approximations we used in our numerical calculations. The following Equations (58)–(61) have already been presented in the previous sections.
The closed system of the equations of motion for the neutrino density matrix and the continuity equation, which reads the equation of the evolution for the photon temperature, in the expanding universe are [37,39]
d ρ p ( t ) d t = ( t H p p ) ρ p ( t ) = i p , ρ p ( t ) + C [ ρ p ( t ) ] ,
d ρ d t = 3 H ( ρ + P ) ,
and analogous Boltzmann equations for anti-neutrinos [37,40], which is not solved in this article since we assume no lepton asymmetry. Here H = 1 m Pl 8 π ρ 3 is the Hubble parameter, p is the Hamiltonian which governs the neutrino oscillation in vacuum and the forward scattering of neutrinos in the e ± , ν , ν ¯ -background, C [ ρ p ( t ) ] is the collision term describing the momentum changing scatterings and annihilations, and [ · , · ] represents the commutator of matrices with a flavor (or mass) index. ρ and P in Equation (59) are the total energy density and the pressure for γ , e ± , ν , ν ¯ , respectively. Including QED finite temperature corrections up to O ( e 3 ) , ρ and P are given by Equation (56) (see also Equations (50), (52) and (54) for the detail of Equation (56)).
The effective Hamiltonian for the neutrino oscillations in vacuum and the forward scattering of neutrinos in the e ± , ν , ν ¯ -background is given by 3
p = M 2 2 p + 2 G F ( N e N e + ) + 2 G F ( N ν N ν ¯ ) 2 2 G F p m W 2 ( E e + P e + E e + + P e + ) 8 2 G F p 3 m Z 2 ( E ν + E ν ¯ ) ,
where GF is the Fermi coupling constant and mW, mZ are the W and Z boson masses, respectively.
The first term in the RHS of Equation (60) denotes neutrino oscillations in vacuum and M 2 is the mass-squared matrix. In the flavor basis, we can write M 2 = U PMNS M diag 2 U PMNS , where M diag 2 = diag ( m ν 1 2 , m ν 2 2 , m ν 3 2 ) . The other terms describe the forward scattering of neutrinos in the background of thermal plasma which comes from one-loop thermal contributions to neutrino self energy. N e ± , N ν , ν ¯ , E e ± , P e ± , E ν , ν ¯ are defined in the flavor basis around the temperature of MeV scale as
N e N e + = diag ( n e n e + , 0 , 0 ) , n e ± = 2 d 3 p ( 2 π ) 3 f e ± ( p ) , N ν N ν ¯ = d 3 p ( 2 π ) 3 ρ p ( t ) ρ ¯ p ( t ) , E e ± + P e ± = diag ( ρ e ± + P e ± , 0 , 0 ) , ρ e ± + P e ± = d 3 p ( 2 π ) 3 E e + p 2 3 E e f e ± ( p ) , E ν + E ν ¯ = d 3 p ( 2 π ) 3 p ρ p ( t ) + ρ ¯ p ( t ) ,
where E e = p 2 + m e 2 + δ m e 2 ( p , T ) and f e ± ( p ) is the distribution function of e ± . δ m e 2 ( p , T ) is the QED finite temperature correction to m e , which is given by Equation (57) up to O ( e 2 ) . Here we neglect the contributions of μ and τ since the densities of these charged particles are significantly suppressed.
In the following, we assume that electrons and positrons are always in thermal equilibrium and follow the Fermi–Dirac distributions since electrons, positrons and photons interact with each other through rapid electromagnetic interactions. In addition we neglect lepton asymmetry since neutrino oscillations leading to flavor equilibrium before the BBN imposes a stringent constraint on this asymmetry [44,45,46,47,48,49,50]. The standard baryogenesis scenarios via the sphaleron process in leptogenesis models predict that the lepton asymmetry is of the order of the current baryon asymmetry, n b / n γ 10 10 , which is much smaller than the above constraint. We also neglect any CP-violating phase in the PMNS matrix for simplicity. Note that from the recent global analysis of neutrino oscillation experiments [51,52], the CP-conserving PMNS matrix is excluded at approximately 3 σ confidence level. Strictly speaking, ignoring the CP-violating phase is inconsistent with the experimental results, but we adopt this assumption to save computational time. In fact, since effects of CP-violating phase on neutrino oscillations are sub-dominant, this ignorance will not affect the resultant neutrino spectra and N eff significantly. Under these assumptions, neutrinos and anti-neutrinos satisfy the same density matrices and the same evolutions in the Universe, ρ p ( t ) = ρ ¯ p ( t ) T , and electrons and positrons follow the same Fermi–Dirac distributions with T γ and no chemical potential.
Note that without lepton asymmetry, N ν N ν ¯ 0 due to ρ p ( t ) = ρ ¯ p ( t ) T ρ ¯ p ( t ) . However, in the following, we neglect it for reducing computational time. We will discuss this uncertainty in Section 3.3. In addition, as in refs. [19,20,21,25], we replace E e ± + P e ± as 4 / 3 E e ± for simplicity. Strictly, this replacement is valid only in the ultra-relativistic limit [38]. However, since in the non-relativistic region E e ± is suppressed by the Boltzmann factor, these difference would be quite small. Ref. [27] reported this difference in N eff is no more than 10 5 .
The final term in the RHS of Equation (58) represents both the momentum conserving and changing collisions of neutrinos with neutrinos, electrons and their anti-particles. In this term, collisions are dominated by two-body reactions 1 + 2 3 + 4 , i.e., C [ ρ p ( t ) ] G F 2 , where G F is the Fermi coupling constant. The detailed formula for C [ ρ p ( t ) ] is given by Equation (48) (see also Equations (38), (40), (41), (43) and (44) in this review and refs. [20,26]). Nine integrals in the collision term (37) can be reduced analytically to two integrals as in Appendix B. We deal with both diagonal and off-diagonal collision terms in Equations (38), (40) and (41) for the processes which involve electrons and positrons, ν e ± ν e ± and ν ν ¯ e e + . On the other hand, we do not treat the off-diagonal terms in Equations (43) and (44) for the self-interactions of neutrinos, ν ν ν ν and ν ν ¯ ν ν ¯ , 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. We treat this collision term from neutrino self-interaction in Equation (A13) of Appendix A. In refs. [26,27], the authors solve kinetic equations for neutrinos including the full collision term at tree level and reported almost the same results with very small difference in N eff , δ N eff 2 × 10 4 [27]. Here, we take into account finite temperature corrections to m e up to O ( e 2 ) in the collision term as E e = p 2 + m e 2 + δ m e 2 ( p , T ) . However, we neglect other sub-leading contributions to the collision term, i.e., other QED corrections to weak interaction rates. We also discuss these uncertainties in Section 3.3.

2.6. Computational Method, Initial Conditions and Values of Neutrino Masses and Mixing

We solve kinetic equations for neutrinos of Equations (58) and (59) with the following comoving variables instead of the cosmic time t, the momentum p, and the photon temperature T γ ,
x = m e a , y = p a , z = T γ a ,
where we choose an arbitrary mass scale in x to be the electron mass m e and a is the scale factor of the universe, normalized as z 1 ( a 1 / T γ ) in high temperature limit. The resultant kinetic equations for neutrinos in the comoving variables are described in Appendix A.
Since the Boltzmann Equations (58) are integro-differential equations due to integrations in the collision terms, their equations were solved by a discretization in a momentum grid y i in refs. [8,9,10,16,19,20,21], by an expansion of the distortions of neutrinos from the Fermi–Dirac distribution in refs. [11,14,15], or by a hybrid method combining the previous two methods in ref. [18]. In this study, we adopt the discretization method we mentioned first and take 100 grid points for y i , equally spaced in the region y i [ 0.02 , 20 ] with the Simpson method. We have used MATLAB ODE solver, in particular, ode15s with an absolute and relative tolerance of 10 6 . In these tolerances, we confirm that numerical errors for relic neutrino spectra and N eff are typically 10 4 or less.
We have numerically estimated the evolution of the density matrix for neutrinos and the photon temperature in x in x x f . We have set x in = m e / 10 MeV as an initial time. Since neutrinos are kept in thermal equilibrium with the electromagnetic plasma at x in , the initial values of density matrix ρ y i in ( x ) are regarded as
ρ y i in ( x ) = diag 1 e y i / z in + 1 , 1 e y i / z in + 1 , 1 e y i / z in + 1 .
The initial dimensionless photon temperature at x in , z in , slightly deviates from 1 because a tiny amount of e ± -pairs have already been annihilated at x in . Due to the entropy conservation of electromagnetic plasma, neutrinos and anti-neutrinos, z in is estimated as in [10],
z in = 1.00003 .
We take x f = 30 as a final time, when the neutrino density matrix and z can be regarded as frozen.
Finally we comment on values of neutrino masses and mixing we use in our numerical simulation. We use the best-fit values in the global analysis in 2019 [53], but assume CP-symmetry, δ CP = 0 . We note that in 2020 their best-fit values are updated [51,52] though their differences are very small. Their parameters include small uncertainties of about 10 % at 3 σ confidence level. Effects of their uncertainties on N eff is investigated in ref. [27] and slightly change N eff by | δ N eff | 10 4 . In our numerical simulation, we confirmed that relic neutrino spectra and the value of N eff with 10 3 precision are the same for both neutrino mass ordering. In the following, we show the results in the normal mass ordering, Δ m 31 2 > 0 , not in the inverted ordering, Δ m 31 2 < 0 , because the results do not change significantly.

3. Effective Number of Neutrino Species N eff

To describe the process of neutrino decoupling, we first numerically solve a set of Equations (58) and (59) and show relic neutrino spectra in the flavor basis. Then we present a precise value of the effective number of neutrino species, N eff = 3.044 , and discuss effects of neutrino oscillations and finite temperature corrections to m e , ρ and P up to O ( e 3 ) on N eff . We also comment on uncertainties of ingredients we ignored in estimating N eff .

3.1. Relic Neutrino Spectra in the Flavor Basis

In the left panel of Figure 2, we show the distortions of the flavor neutrino spectra for a comoving momentum ( y = 5 ) , where we plot the neutrino spectra f ν α / f eq as a function of the normalized cosmic scale factor x. f eq ( y ) is the neutrino distribution function if neutrinos decoupled instantaneously and all e ± -pairs annihilated into photons,
f eq ( y ) = 1 e y + 1 .
At high temperature with ( x 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 x 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 [ ρ p ( t ) ] becomes ineffective and the distortions are frozen.
The difference between the ν e spectrum and the ν μ , τ spectrum without flavor mixing 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.
In the right panel of Figure 2, we show the frozen values of the flavor neutrino spectra f ν α / f eq as a function of a comoving momentum y for both cases with and without neutrino mixing. This figure shows the fact that neutrinos with higher energies interact with electrons and positrons until a later epoch. In addition, we see neutrino oscillations tend to equilibrate the flavor neutrino distortions. Although the neutrino spectra f ν α / f 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 neutrino decoupling process.

3.2. Value of the Effective Number of Neutrino Species N eff

The effective number of neutrinos N eff can be rewritten,
N eff = ( 11 / 4 ) 1 / 3 z 4 3 + δ ρ ν e ρ ν eq + δ ρ ν μ ρ ν eq + δ ρ ν τ ρ ν eq ,
where δ ρ ν α = ρ ν α ρ ν eq and ρ ν eq = d 3 p ( 2 π ) 3 p f eq . In Table 2 and Table 3, we present final values (at x f = 30 ) of the dimensionless photon temperature z fin , the difference of energy densities and number densities of flavor neutrinos from those where neutrinos decoupled instantaneously denoted by ρ ν eq = d 3 p ( 2 π ) 3 p f eq and n ν eq = d 3 p ( 2 π ) 3 f eq , and the effective number of neutrinos N eff .
By comparing values of N eff in the cases without QED corrections and with QED corrections to m e , ρ and P up to O ( e 2 ) and O ( e 3 ) in Table 2, we find that the QED corrections at O ( e 2 ) and O ( e 3 ) shift N eff by + 0.01 and 0.00095 , respectively, which is very close to the value estimated in the instantaneous decoupling limit [24].
In the cases with neutrino mixing, Table 3 shows that the energy densities of μ , τ -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 eff = 3.04391 3.044 with QED corrections to m e , ρ and P up to O ( e 3 ) . Since the blocking factor for electron neutrinos, ( 1 f ν 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 eff is 0.00056 , which agrees with the results of previous works [12,20,23].
To conclude, our numerical calculation with neutrino oscillations and QED finite temperature corrections to m e , ρ and P up to O ( e 3 ) finds N eff = 3.044 . This value is in excellent agreement with later independent works [26,27].

3.3. Discussions of Uncertainties in N eff

We comment on possible errors of the results for relic neutrino spectra and N eff due to approximations in Equations (58) and (59) and the choice of physical parameters. Our numerical calculations converge very well since we have directly computed N eff in the mass basis as will be done in the next section and obtained N eff = 3.04388 3.044 .
First we neglect the off-diagonal parts for neutrino self-interactions in the collision term, ν ν ¯ ν ν ¯ and ν ν ν ν . Later, in refs. [26,27], the authors solve kinetic equations for neutrinos including their off-diagonal parts in the collision term and report the difference in N eff is δ N eff 2 × 10 4 [27]. We also neglect the O ( e 2 ) logarithmic terms and terms above O ( e 4 ) in QED finite temperature corrections to m e , ρ and P. Their corrections to ρ and P are reported to contribute δ N eff < 10 4 to N eff in refs. [24,27]. Though their corrections to m e are not taken into account, the corrections to m e even at O ( e 2 ) contribute δ N eff 10 4 to N eff [27] and we have also confirmed it.
The neutrino masses and mixing parameters contain 10–20% uncertainties at 3 σ confidence level. Since in our estimations, neutrino oscillations contribute + 0.0005 to N eff , their uncertainties are expected to be quite small. In ref. [27], the authors report that their uncertainties are δ N eff 10 4 . We also neglect the CP-violating phase δ CP in the PMNS matrix. No one has yet computed precise neutrino evolution in the decoupling including three-flavor oscillations with CP violating phase. However, since effect of the CP-violating phase on neutrino oscillations is sub-dominant, we expect neutrino and anti-neutrino spectra might not change significantly. In addition, the total energy density, i.e., N eff would change much less than 0.0005 since the changes for the energy densities of neutrinos and anti-neutrinos would be canceled out. See also discussion in Appendix F of ref. [26] and ref. [54]. Other physical parameters for electroweak interaction are measured very precisely and will not affect neutrino spectra and N eff .
However, QED corrections to weak interaction rates at order O ( e 2 G F 2 ) and forward scattering of neutrinos via their self-interactions have not been precisely taken into account in the whole literature so far.

3.3.1. QED Corrections to Weak Interaction Rates at Order e 2 G F 2

QED interactions also modify the weak interaction rates in the collision term C [ ρ p ( t ) ] and the Hamiltonian for the forward scattering of neutrinos (60) at order e 2 G F 2 in addition to the modification of the energy density and pressure for electromagnetic plasma, ρ and P. These corrections are partially taken into account by considering thermal QED corrections on m e so far. See also Section 3.1.2 in ref. [27].
QED corrections to the weak interaction rates (see also the diagrams in Figure 3) are categorized as (i) additional photon emission and absorption, (ii) corrections to the dispersion relation for external e ± , (iii) vertex corrections, and (iv) corrections mediated by photon propagator. The interference among the weak interaction at leading order G F and corrections (i)–(iv) produce modifications to the weak interaction rates at the next-to-leading order O ( e 2 G F 2 ) .
The correction (i) might be the most dominant contribution to N eff since the photon emission processes, e.g., e + e ν ν ¯ γ , would not be suppressed by the distribution function of photons in the Boltzmann equations. The photon emission processes reduce N eff . However, there are many processes in the categories (ii), (iii) and (iv). In total, these contributions to N eff might be as large as that from the correction (i).
For category (ii), corrections to the dispersion relation for e ± produce a thermal electron mass as Equation (57). One can incorporate corrections (i) in the weak interaction rates by shifting m e 2 m e 2 + δ m e ( 2 ) 2 ( p , T ) , but it is numerically difficult to take into account the momentum-dependent part of δ m e ( 2 ) 2 ( p , T ) , which corresponds to the logarithmic O ( e 2 ) corrections to m e . These logarithmic O ( e 2 ) corrections to m e are less than 10 % of corrections at O ( e 2 ) to m e around neutrino decoupling [43], and corrections even at leading O ( e 2 ) to the weak interaction rates (i.e., δ m e ( 2 ) ( T ) ) contributes N eff < 10 4 to N eff [27] and we confirmed it. Thus, we would properly be able to incorporate corrections (i) to N eff with 10 4 precision. However, we should carefully derive these corrections to the weak interaction rates and consider effects of the logarithmic O ( e 2 ) corrections and other sub-dominant neglected contributions in the collision term in the future.
For categories (i), (iii) and (iv), corrections to the weak interaction rates are typically momentum-dependent. It would be quite difficult to solve the Boltzmann equation, which is the integro-differential equation, including such momentum-dependent corrections. In ref. [55], the authors consider energy loss rate of a stellar plasma, including corrections on e e + ν ν ¯ at order O ( e 2 G F 2 ) and found such corrections modify the energy loss rate of a stellar plasma by a few percent. In ref. [23], the author suggests δ N eff 0.0007 due to correction (i) by roughly extrapolating the results in ref. [55] and using a precise and simple evaluation method of N eff proposed in ref. [23]. The contributions of (i), (iii) and (iv) to N eff should be evaluated in the future in a more precise way.

3.3.2. Forward Scattering of Neutrinos via Their Self-Interactions

In the Hamiltonian (60) in the Boltzmann Equations (58), the forward scattering terms of neutrinos via their self-interactions correspond to
p 2 G F ( N ν N ν ¯ ) 8 2 G F p 3 m Z 2 ( E ν + E ν ¯ ) .
Even in the case without lepton asymmetry, N ν N ν ¯ 0 due to ρ p ( t ) = ρ ¯ p ( t ) T ρ ¯ p in general, where N ν N ν ¯ is
N ν N ν ¯ = d 3 p ( 2 π ) 3 ρ p ( t ) ρ ¯ p ( t ) .
Though ρ p ( t ) ρ ¯ p ( t ) might be small, forward scattering via neutrino self-interactions could be more dominant than neutrino oscillation in vacuum, with a typical dimensional analysis,
2 G F p 3 10 11 MeV G F 10 5 GeV 2 p 1 MeV 3 M 2 2 p 10 14 MeV M 0.1 eV 2 1 MeV p .
In ref. [28], the authors suggest forward scattering of neutrinos via their self-interactions contributes δ N eff + ( 1 5 ) × 10 4 to N eff by solving a simplified kinetic equations for neutrinos. In the future, relic neutrino spectra and N eff should be estimated including the above forward scattering of neutrinos more precisely.
Though recent estimations might contain uncertainties of | δ N eff | ( 10 3 10 4 ) in N eff , N eff = 3.044 would still be one of very good reference values in N eff .

4. Relic Cosmic Neutrino Spectra in the Current Homogeneous and Isotropic Universe

In the current universe, two neutrino species at least are non-relativistic. Then relic neutrino spectra in the mass basis will be important observable to detect the C ν B in a direct way as discussed in Section 2.1. In this section, we present the spectrum (as a function of comoving momenta), number density and energy density of the C ν B in the current homogeneous and isotropic universe, including non-thermal distortions due to e ± -annihilation during neutrino decoupling.

4.1. Relic Neutrino Spectra in the Mass Basis

We present relic neutrino spectra in the mass basis by solving a set of Equations (58) and (59) in the mass basis directly. We can also obtain the same result by transforming relic neutrino spectra in the flavor basis through Equation (8).
In the mass basis, the neutral and charged currents including left-handed neutrino fields in Equation (24) are given by, using ν α = i = 1 , 2 , 3 U α i ν i as in Equation (6),
J ν ν = 1 4 cos θ W α = e , μ , τ ν ¯ α γ μ ( 1 γ 5 ) ν α = 1 4 cos θ W i = 1 , 2 , 3 ν ¯ i γ μ ( 1 γ 5 ) ν i , J e ν e μ = 1 2 2 ν ¯ e γ μ ( 1 γ 5 ) e = 1 2 2 i = 1 3 U e i ν ¯ i γ μ ( 1 γ 5 ) e .
Then, using the relations of Equation (70) and (34), we obtain the 4-point interaction Hamiltonian (32) in the mass basis
H int e ν | mass = G F 2 d x 3 ν ¯ γ μ ( 1 γ 5 ) Z L ν e ¯ γ μ ( 1 γ 5 ) e + ν ¯ γ μ ( 1 γ 5 ) Z R ν e ¯ γ μ ( 1 + γ 5 ) e , H int ν | mass = G F 4 2 d x 3 ν ¯ γ μ ( 1 γ 5 ) ν ν ¯ γ μ ( 1 γ 5 ) ν ,
with
ν = ν 1 ν 2 ν 3 , Z L = g ˜ L + U e 1 U e 1 U e 1 U e 2 U e 1 U e 3 U e 2 U e 1 g ˜ L + U e 2 U e 2 U e 2 U e 3 U e 3 U e 1 U e 3 U e 2 g ˜ L + U e 3 U e 3 , Z R = Y R = sin 2 θ W × 1 .
Then we obtain the Boltzmann equation for the neutrino density matrix in the mass basis after replacements of Y L , R Z L , R and p U PMNS p U PMNS analogous to M diag 2 = U PMNS M 2 U PMNS in Equation (58) for the flavor basis.
In the left panel of Figure 4, we show the evolution of the neutrino spectra, f ν i / f eq , for a comoving momentum ( y = 5 ) as a function of the normalized scale factor x. In the right panel of Figure 4, we show the asymptotic values of the neutrino spectra 4  f ν i / f eq as a function of y. The differences of distortions for each neutrino species arise from the charged current interactions between neutrinos and electrons weighted by the PMNS matrix with mass species i, U e i , as in Equation (71). Note that neutral currents between neutrinos in the mass basis are the same as that in the flavor basis except for the subscript, J ν ν μ = α = e , μ , τ ν ¯ α γ μ ( 1 γ 5 ) ν α = α = 1 , 2 , 3 ν ¯ i γ μ ( 1 γ 5 ) ν i . Then the scattering and annihilation among neutrinos and electrons and their anti-particles induce the spectral distortions in Figure 4.
Finally we comment on N eff . After we directly solve a set of Equations (58) and (59) in the mass basis, including vacuum three-flavor neutrino oscillations, forward scatterings in e ± -background, and QED corrections to m e , ρ and P up to O ( e 3 ) , we find N eff = 3.04388 , which is an excellent agreement with our calculation in the flavor basis. The tiny difference from N eff in the flavor basis may come from ignoring the off-diagonal parts for self-interaction processes in the Boltzmann equations and/or numerical errors.

4.2. Neutrino Number Density and Energy Density in the Current Homogeneous and Isotropic Universe

In Table 4, we show the final values of the dimensionless photon temperature z fin , the relativistic energy densities ρ ν i / ρ ν eq and number densities n ν i / n ν eq of neutrinos in the mass basis after neutrino decoupling. Note that the expression of energy density for a relativistic particle is not applicable to the first and second heaviest neutrinos today because they are non-relativistic in the current universe.
After neutrino decoupling, the neutrino momentum distribution in the homogeneous and isotropic universe can be parametrized as
f ν i ( p , t ) = 1 e | p | / T ˜ ν ( t ) + 1 1 + δ f ν i ( p , t ) .
T ν ( t ) is the effective neutrino temperature, which is a ( t ) 1 and normalized as T ν T γ in high temperature limit. Under this definition of T ν ( t ) , neutrino spectral distortions, δ f ν i ( p , t ) , can be rewritten as δ f ν i ( y ) given in the right panel of Figure 4. At t 0 = 4.35 × 10 17 s in the current universe, T ν ( t 0 ) satisfies
T γ ( t 0 ) T ν ( t 0 ) = z fin = 1.39797 ,
where T γ ( t 0 ) 2.7255 K is the effective photon temperature in the current universe [56]. Then the effective neutrino temperature in the current universe is
T ν ( t 0 ) = 1.9496 K .
Neutrino number density and energy density per one degree of freedom in the current universe are also parametrized as
n ν i ( t 0 ) = d 3 p ( 2 π ) 3 f ν i ( p , t ) , = n 0 ( 1 + δ n ¯ ν i ) , ρ ν i ( t 0 ) = d 3 p ( 2 π ) 3 E ν i f ν i ( p , t ) , = m i n ν i for non - relativistic ν i ρ 0 ( 1 + δ ρ ¯ ν i ) for relativistic ν i ,
where n 0 and ρ 0 are given by
n 0 = d 3 p ( 2 π ) 3 1 e | p | / T ν ( t 0 ) + 1 = 3 ζ ( 3 ) 4 π 2 T ν ( t 0 ) 3 = 56.376 cm 3 , ρ 0 = d 3 p ( 2 π ) 3 | p | e | p | / T ν ( t 0 ) + 1 = 7 π 2 240 T ν ( t 0 ) 4 = 29.848 meV cm 3 .
Then δ n ¯ ν i and δ ρ ¯ ν i are given in Table 4. The values of neutrino number density in the current universe are listed in Table 5.
In the current universe, two species of cosmic relic neutrinos at least are non-relativistic because of T ν ( t 0 ) Δ m 21 2 8.6 meV , | Δ m 31 2 | 50 meV . On the other hand, the lightest neutrinos might be relativistic in the current universe because the lightest neutrino mass is not yet determined. In Table 6, we show energy density for the lightest neutrinos in the case of m lightest p 0 3.15 T ν ( t 0 ) . Here we consider both the normal mass ordering, m ν 3 > m ν 2 > m ν 1 , and the inverted mass ordering, m ν 2 > m ν 1 > m ν 3 .
To estimate the effects of e ± -annihilation into neutrinos during neutrino decoupling on neutrino number density and energy density, it is useful to compare the neutrino number density and relativistic energy density per one degree of freedom in the case when all e ± -pairs annihilate into photons, n 0 and ρ 0 , respectively,
n 0 = 3 ζ ( 3 ) 4 π 2 T ν ( t 0 ) 3 = 56.01 cm 3 ,
ρ 0 = 7 π 2 240 T ν ( t 0 ) 4 = 29.65 cm 3 ,
where T γ ( t 0 ) / T ν ( t 0 ) = ( 11 / 4 ) 1 / 3 . We show the deviation of neutrino number density from the case when all e ± -pairs annihilate into photons, δ n ν i d n ν i / n 0 1 , in Table 7. The number densities for all neutrino species are enhanced by about 1 % due to e ± -annihilations to neutrinos during neutrino decoupling and the number density for ν 1 is most efficiently enhanced.

4.3. Helicity of Relic Majorana Neutrinos vs. Dirac Neutrinos

The weak interaction is chiral, which is manifest in the Lagrangian. Due to its chirality, the left-chiral states for SM fermions interact with the weak bosons while the right-chiral states do not. In the early universe, only left-chiral neutrinos and right-chiral anti-neutrinos, i.e., left-handed neutrinos and right-handed anti-neutrinos are produced via the weak interaction. Note that chirality is different from helicity in general, which is defined as the projection of the spin vector onto the momentum vector.
During free streaming of relic neutrinos after their decoupling, the chirality for non-relativistic neutrinos is not conserved since the chiral symmetry in the free neutrino Lagrangian is broken due to their masses. On the other hand, the helicity for relic neutrinos is conserved in the homogeneous and isotropic universe. Thus, we should estimate the spectrum for each helicity state of relic cosmic neutrinos in the current universe.
In the early universe, both chirality and helicity for relic neutrinos are conserved and then neutrino helicity and chirality have one-to-one correspondence since neutrinos are approximately massless in the early universe. We define left (right) helical neutrinos with helicity s ν = 1 / 2 ( + 1 / 2 ) such that they correspond to left (right) handed neutrinos in the early universe. Then the spectra for the left-handed neutrinos (right-handed anti-neutrinos) produced in the early universe are translated into the left-helical neutrinos (right-helical anti-neutrinos) [34],
f ν i ( p ν , s ν = 1 / 2 ) = f ν i ( p ν , t ) , f ν i ( p ν , s ν = + 1 / 2 ) 0 , f ν ¯ i ( p ν , s ν = 1 / 2 ) 0 , f ν ¯ i ( p ν , s ν = + 1 / 2 ) = f ν ¯ i ( p ν , t ) f ν i ( p ν , t ) ,
where f ν i ( p ν , t ) is given by Equation (73) and f ν ¯ i ( p ν , t ) f ν i ( p ν , t ) if we neglect lepton asymmetry. Here right-helical neutrinos, ν i with s ν = + 1 / 2 , (left-helical anti-neutrinos, ν ¯ i with s ν = + 1 / 2 ,) corresponds to right-handed neutrinos (left-handed anti-neutrinos), which are sterile states. We assume sterile neutrinos are not produced in the early universe due to very weak interactions with the SM particles or have already decayed if sterile neutrinos are right-handed heavy Majorana particles as required for the see-saw mechanism.
For Majorana neutrinos, right-handed active anti-neutrinos are regarded as right-handed active neutrinos due to the lepton number violation. Then f ν i ( p ν , s ν ) for ν i are given by
f ν i ( p ν , s ν = 1 / 2 ) = f ν i ( p ν , t ) , f ν i s ( p ν , s ν = + 1 / 2 ) 0 , f ν i s ( p ν , s ν = 1 / 2 ) 0 , f ν i ( p ν , s ν = + 1 / 2 ) = f ν ¯ i ( p ν , t ) f ν i ( p ν , t ) ,
where ν i s denotes a sterile state of neutrino. Note that even in the case of Majorana neutrinos lepton asymmetry can be interpreted as chiral asymmetry between left-handed and right-handed neutrinos. Then f ν ¯ i ( p ν , t ) and f ν i ( p ν , t ) are different strictly speaking but almost the same approximately.
For Dirac neutrinos, since right-handed neutrinos and left-handed anti-neutrinos are sterile, f ν i ( p ν , s ν ) for ν i are given by
f ν i ( p ν , s ν = 1 / 2 ) = f ν i ( p ν , t ) , f ν i s ( p ν , s ν = + 1 / 2 ) 0 , f ν ¯ i s ( p ν , s ν = 1 / 2 ) 0 , f ν ¯ i ( p ν , s ν = + 1 / 2 ) = f ν ¯ i ( p ν , t ) f ν i ( p ν , t ) ,
where ν ¯ i s denotes a sterile state of anti-neutrino.
From Equations (81) and (82), the magnitude of relic neutrino spectra summed over helicity for Majorana and Dirac neutrinos differ by a factor of two, which is first pointed out in ref. [34],
s ν = ± 1 / 2 f ν i ( p ν , s ν ) 2 f ν i ( p ν , t ) for Majorana ν i f ν i ( p ν , t ) for Dirac ν i .
Then number density and energy density summed over helicity for Majorana and Dirac neutrinos also differ by a factor of two,
s ν = ± 1 / 2 n ν i ( s ν ) 2 n ν i for Majorana ν i n ν i for Dirac ν i , s ν = ± 1 / 2 ρ ν i ( s ν ) 2 ρ ν i for Majorana ν i ρ ν i for Dirac ν i .

5. Implications for the Capture Rates on Cosmic Neutrino Capture on Tritium

Finally we discuss how neutrino spectral distortions from e ± -annihilations during neutrino decoupling affect direct detection of the C ν B on tritium target, with emphasis on the PTOLEMY-type experiment [30,31], where cosmic neutrinos can be captured on tritium by the inverse beta decay process without threshold energy for neutrinos, ν i + 3 H e + 3 He . Tritium is one of appropriate candidates for the target because of its availability, high capture rate for neutrinos, low Q-value and long half lifetime of t 1 / 2 = 12.32 years. Here we take 100 g of tritium as the target. We take into account gravitational clustering for cosmic neutrinos in our Galaxy and nearby galaxies because we would observe the C ν B directly inside our Galaxy. We also comment on gravitational helicity flipping and annual modulation for the C ν B. Then we discuss the potential of direct measurements of such cosmological effects although it would be still extremely difficult to observe such effects directly. In particular, we compute the capture rates of cosmic relic neutrinos on tritium, including such cosmological effects.

5.1. Gravitational Effects for the C ν B

5.1.1. Clustering for the C ν B by Our Galaxy and Nearby Galaxies

Near the Earth, non-relativistic relic neutrinos cluster locally in the gravitational potential of our Galaxy and nearby galaxies. Then the local distribution function is distorted and the local number density is enhanced compared with the global distribution function and number density. The local number density for relic neutrinos in the current universe is described as
n ν i loc = n ν i ( 1 + δ n ν i c ) ,
where δ n ν i c is an enhancement factor by the gravitational attraction by galaxies, which is estimated in refs. [36,57,58,59,60,61]. For reference, we display some of these values, estimated in a recent numerical study [36], in Table 8, where the authors consider the gravitational potential in the Milky Way, Virgo cluster, and Andromeda galaxy. Note that so far, when evaluating values of δ n ν i c , effects of e ± -annihilations into ν , ν ¯ during neutrino decoupling have not been taken into account simultaneously. For m ν i < 0.15 eV , spectral distortions to the momentum distributions for relic cosmic neutrinos by the gravitational clustering have not also been explicitly estimated (see ref. [58] for spectral distortions by gravitational clustering for relic neutrinos with m ν i 0.15 eV ).
In the following, we discuss only the case where δ n ν i c < 1 and the lightest neutrino mass is quite small because the Planck satellite suggests m ν < 0.12 eV . Then the local number density for relic neutrino can be parametrized as, using linear approximation,
n ν i loc n 0 ( 1 + δ n ν i c + δ n ν i d ) ,
where δ n ν i d is the enhancement factor by e ± -annihilations into ν and ν ¯ during neutrino decoupling given in Table 7.

5.1.2. Helicity Flipping and Annual Modulation for the C ν B

We shortly comment on gravitational helicity flipping and annual modulations for relic neutrinos. Gravitational clustering for massive neutrinos may induce mixing of relic neutrino helicity [34,35,62] since the direction of neutrino momentum would change in the gravitational potential for our Galaxy whereas its spin does not. Although the quantitative calculations have not yet been achieved, the capture rates on tritium would not change since their capture rates depend on neutrino number density summed over helicities at leading order as we will see in the next section. In addition, an annual modulation for relic neutrinos might occur in a direct detection experiment for the C ν B since their velocity relative to the Earth could be anisotropic due to neutrino clustering and the gravitational focusing for the C ν B by the Sun could also occur. The former effect is negligible since the capture rates on tritium target are independent of neutrino velocity as we will see in the next section. The latter effect is expected to change the capture rates by much less than 1% for m ν < 0.15 meV [63]. In the following, we neglect helicity flipping and annual modulation for relic neutrinos.

5.2. Precise Capture Rates on Tritium including Sub-Dominant Cosmological Effects

In Table 6, non-thermal distortions during neutrino decoupling enhance the number density of the C ν B by about 1 % . To properly incorporate such effects into the capture rates of the C ν B on tritium, we discuss the formula of their capture rate with 1 % precision.
Cosmic relic neutrinos can be captured on tritium by the following inverse beta decay process,
ν i + 3 H 3 He + e .
The total capture rate for the C ν B in this process, Γ C ν B , can be written
Γ C ν B = i = 1 N ν Γ i ,
where N ν is the number of (mass) species of neutrinos. Γ i is the capture rate for a given mass-eigenstate of neutrino ν i , given by
Γ i = N T s ν = ± 1 / 2 d 3 p ν ( 2 π ) 3 σ ν i ( p ν , s ν ) v ν i f ν i loc ( p ν , s ν ) ,
where N T = M T / M 3 H is the number of tritium, M T is the total tritium mass in the experimental setup, and M 3 H 2809.432 MeV is the atomic mass of tritium. s ν , v ν i = | p ν | / E ν i and σ ν i are helicity, velocity and the total cross section in the inverse beta decay on tritium, respectively. f ν i loc ( p ν , s ν ) is the local momentum distribution for relic cosmic neutrinos around the Earth, which satisfies n ν i loc ( s ν ) = d p ν 3 ( 2 π ) 3 f ν i loc ( p ν , s ν ) .
In cosmic neutrino capture on tritium, the spins of the outgoing electron and nucleus would not be measured. In addition, the spin of the initial nucleus would not be identified either. On the other hand, the helicity state for cosmic neutrinos in the Dirac case is polarized as in Section 4.3. Then we compute the spin-polarized cross section for ν i . After averaging over the spin of 3 H and summing over the spin of outgoing e and 3 He , the formulae of σ ν i ( p ν , s ν ) with 1 % precision reduces to (see Appendix D for detail calculations)
σ ν i ( p ν , s ν ) G F 2 2 π | V u d | 2 | U e i | 2 m 3 He m 3 H v ν i f F 2 + g A 2 g V 2 g G T 2 × F ( 2 , E e ) E e | p e | ( 1 2 s ν v ν i ) ,
where V u d 0.9740 is a component of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, m 3 H 2808.921 MeV and m 3 He 2808.391 MeV are the nuclear masses of 3 H and 3 He , g A 1.2723 and g V 1 are the axial and vector coupling constant, and f F 0.9998 and g GT 3 × ( 0.9511 ± 0.0013 ) are the reduced matrix elements of the Fermi and Gamow-Teller (GT) operators, respectively. The Fermi function F ( Z , E e ) is an enhancement factor by the Coulombic attraction of the outgoing electron and proton, which is approximately given by [64]
F ( Z , E e ) = 2 π α Z E e / | p e | 1 e 2 π α Z E e / | p e | ,
where α 137.036 is the fine structure constant. Z is the atomic number of the daughter nucleus and Z = 2 for 3 He . The energy and momentum for an emitted electron E e and p e depend on the neutrino masses and momenta strictly because of momentum conservation in the inverse β -decay process. However, since the contributions of the neutrino masses and momenta to E e and p e are very small, E e and | p e | are approximately given by (see Appendix C for details)
E e K end 0 + m e + E ν i K end 0 + m e , | p e | = E e 2 m e 2 ,
where K end 0 is the beta decay endpoint kinetic energy for massless neutrinos given by
K end 0 = ( m 3 H m e ) 2 m 3 He 2 2 m 3 H 18.6 keV .
E ν i is so small compared to K end 0 and m e that we can safely neglect E ν i in Equation (92).
Then we obtain Γ i with 1 % precision substituting Equation (90) into Equation (89),
Γ i N T G F 2 2 π | V u d | 2 | U e i | 2 m 3 He m 3 H f F 2 + g A 2 g V 2 g G T 2 × F ( 2 , E e ) E e | p e | s ν = ± 1 / 2 n ν i ( s ν ) 2 s ν v ν i ,
where v ν i is the (unnormalized) average magnitude of velocity for ν i given by
v ν i = d 3 p ν ( 2 π ) 3 f ν i ( p ν , s ν ) v ν i .
Typically, v ν i contributes more than 1 % to Γ ν i . If m ν i 100 meV , due to v ν i p 0 / m ν i 0.01 , we can drop v ν i in the formula of Equation (94) with 1 % precision. Here p 0 3.15 T ν ( t 0 ) 0.53 meV is the average momentum of the C ν B in the current universe. We also comment on whether we can use further approximations with 1 % precision to write Equation (94) into a simpler form. For massless neutrinos, due to v ν i = | p ν i | / E ν i = 1 , the (unnormalized) velocity is written as v ν i = n ν i . For non-relativistic neutrinos ( m ν 10 meV ) , due to v ν i 1 , v ν i is approximately written as v ν i d 3 p / ( 2 π ν 3 ) f ν 0 ( p , t 0 ) | p ν | / E ν i , where f ν 0 ( p ν , t 0 ) = [ exp ( p ν / T ν ( t 0 ) ) + 1 ] 1 and T ν ( t 0 ) / T γ ( t 0 ) = ( 4 / 11 ) 1 / 3 . We note that gravitational helicity flipping for massive neutrinos by neutrino clustering would be negligible since the helicity-dependent part in Γ i is already suppressed by v ν i .

5.2.1. Majorana vs. Dirac Neutrinos

For non-relativistic neutrinos, i.e., v i 1 , if we set v ν i = 0 in Equation (94), Γ i is porportional to s ν n ν i and left-helical and right-helical components for relic neutrinos interact with tritium with the same magnitude via the weak interaction. Then the capture rate on tritium for Majorana neutrinos Γ i M is twice that for Dirac neutrinos [34],
Γ i M | v ν i 1 2 Γ i D | v ν i 1 .
On the other hand, for relativistic neutrinos, i.e., v i 1 , only the left-helical neutrinos interact with tritium via the weak interaction since helicity coincides with chirality in the relativistic limit. Then in both Majorana and Dirac cases, the capture rates are the same [35],
Γ i M | v ν i 1 Γ i D | v ν i 1 .
Note again that the approximations in Equations (96) and (97) might not be valid for the capture rates with 1 % precision. To estimate the capture rates with 1 % precision, the term that depends on v ν i in Equation (94) should be included precisely.

5.2.2. Values of the Capture Rates on Tritium with m lightest = 0

For references, we show values of the capture rates including cosmological effects discussed in Section 4.2 and Section 5.1 in the case of m lightest = 0 . We choose other neutrino masses and their ordering to satisfy the observed values of neutrino squared-mass differences from neutrino oscillation experiments [51,52],
Normal Ordering ( NO ) : Δ m 21 2 ( 8.6 meV ) 2 Δ m 31 2 ( 50 meV ) 2 Inverted Ordering ( IO ) : Δ m 21 2 ( 8.6 meV ) 2 Δ m 32 2 ( 50 meV ) 2
In both neutrino mass ordering we take the following values of the PMNS matrix,
| U e 1 | 2 0.681 , | U e 2 | 2 0.297 , | U e 3 | 2 0.0222 .
Note that neutrino squared-mass differences and neutrino mixing parameters currently include a few percent (about 10 % ) uncertainties even at 1 σ ( 3 σ ) confidence level.
In Table 9, we show values of the capture rates on 100 g of tritium in both the cases of NO and IO for Majorana and Dirac neutrinos with m lightest = 0 . δ Γ i d denotes the differences between the cases with and without effects of e ± -annihilation during neutrino decoupling and δ Γ i c denotes the differences with and without gravitational clustering for relic neutrinos in nearby galaxies.
For Majorana neutrinos, the capture rates for the first and second heaviest neutrinos are slightly less than twice those for Dirac neutrinos because of v ν i 0 . On the other hand, the capture rates for massless (or almost massless) neutrinos in the cases of Majorana and Dirac neutrinos are the same because of v ν i 1 .

5.2.3. Discussions on Exposure and Uncertainties in the Capture Rates

In this section, we discuss the required amount of tritium to observe the sub-leading cosmological effects themselves, δ Γ i c , d , and the estimated error of the capture rates for relic neutrinos on tritium in more detail.
To observe δ Γ i c , d , we need a large number of events to satisfy typically
i δ Γ i c , d T Γ C ν B T + Γ background T 1 ,
where T is the exposure time and Γ background is a background rate. Even if the background is successfully removed, we need 10 2 –10 4 events of the C ν B signal ( Γ C ν B T 10 2 –10 4 ) because of δ Γ i c , d ( 0.1 0.01 ) × Γ i for i m ν i < 0.12 eV . This requirement corresponds to the need for 10–10 3 kg yr of exposure of tritium. Currently, it is extremely difficult to obtain such amount of the exposure. In the next Section 5.3, we comment on β -decay background, which is one of main background in cosmic neutrino capture on tritium.
The estimated error of the neutrino capture rates mainly comes from the uncertainties of the neutrino mixing parameter, | U e i | 2 , and the undetermined value of the lightest neutrino mass, m lightest . The current errors of PMNS matrix are about a few percent (about 10 % ) at 1 σ ( 3 σ ) confidence level [51,52]. The current upper bound of m lightest is ≲0.8 eV [65]. Thus, unfortunately, it is still difficult to incorporate cosmological sub-dominant contributions into the value of Γ ν i precisely. However, δ Γ i c , d for m lightest = 0 is correctly estimated since uncertainties of | U e i | are canceled out in δ Γ i c , d . Future neutrino oscillation experiments will reduce uncertainties of PMNS matrix (see, e.g., [66,67,68]). In addition, measurement of large β -decay background in the PTOLEMY-type experiment might determine the value of m lightest very precisely [31].
We also note that the theoretical calculation of g GT still includes the uncertainty of a few %, although the estimation of g G T through the observation of the tritium half-life and the value of the Fermi operator, f F , only involves uncertainty of 0.1 % [69].
For a large value of m lightest , gravitational clustering effects of relic neutrinos are typically more dominant than effects of e ± -annihilation during neutrino decoupling. Although the C ν B itself with a large value of m lightest would be easier to observe due to a large gravitational clustering, it is also a very difficult task to distinguish the effects of e ± -annihilation during neutrino decoupling from gravitational clustering effect of relic neutrinos.
Based on the evaluation in this section, it is still extremely difficult to observe e ± -annihilation during neutrino decoupling in the PTOLEMY-type experiment. However, the precise capture rates including cosmological sub-dominant contributions might be useful to distinguish the SM from physics beyond the SM properly in the future.

5.3. β -Decay Background and the Energy Resolution of the Detector to Distinguish the C ν B Signal from It

Finally we comment on β -background and the required energy resolution of the detector to distinguish the C ν B signal from this background, which is one of main difficulties to observe the C ν B directly in the inverse β -decay process.
The main background comes from tritium β -decay process,
3 H 3 He + e + ν ¯ i .
The β -decay spectrum and the capture rate for the β -decay process are given by [70] (see also Appendix D)
d Γ β d E e = N T G F 2 2 π 3 | V u d | 2 | U e i | 2 m 3 He m 3 H f F 2 + g A 2 g V 2 g G T 2 × F ( 2 , E e ) E e | p e | i = 1 3 | U e i | 2 H ( E e , m ν i ) ,
where
H ( E e , m ν i ) = 1 m e 2 / ( E e m 3 H ) ( 1 2 E e / m 3 H + m e 2 / m 3 H 2 ) 2 ( E e max , i E e ) E e max , i E e + 2 m ν i m 3 He m 3 H × E e max , i E e + m ν i m 3 H ( m 3 He + m ν i ) ,
E e max , i is the maximal energy of the emitted electron for 3 H 3 He + e + ν ¯ i , where the electron is emitted in opposite direction to both 3 He and ν e ¯ (see also Appendix C),
E e max , i K end 0 + m e m ν i .
Then the maximal energy for the emitted electron in the β -decay process called the energy at β -decay endpoint is
E e end K end 0 + m e m ν lightest ,
where m lightest is the lightest neutrino mass. We can see that the β -decay spectrum d Γ β / d E e vanishes for E e = E e end . Then the total tritium β -decay rate is obtained as
Γ β = m e E e end d E e d Γ β d E e 10 24 M T 100 g yr 1 .
Since the event number of β -decay background is extremely larger than that of the C ν B signal, we must distinguish the two signals clearly.
To distinguish the C ν B signal and β -decay background, we need a tiny energy resolution of the detector Δ . The energy resolution of a detector characterizes the smallest separation where two signals can be distinguished. The β -decay background closest to the C ν B signal is the electron signal with the maximal energy E e max . To distinguish the C ν B signal for a mass species ν i from β -decay background near the endpoint, the required energy resolution Δ i is expected to be (see Appendix C for details)
Δ i E e C ν B , i E e end m lightest + E ν i ,
where E e C ν B , i is the emitted electron energy from the C ν B signal, ν i + 3 H e + 3 He given by Equation (92).
To take into account the energy resolution of the detector Δ in the spectrum and the number of events for the C ν B signal and the β -decay background, we model the would-be observed spectrum of the emitted electron as a Gaussian-smeared version of the actual spectrum. This is achieved by convolving both the C ν B signal and the β -decay background with a Gaussian of full width at half maximum (FWHM) equal to Δ = 8 ln 2 σ , where σ is the Gaussian standard deviation,
d Γ i d E e = 1 2 π σ d E e Γ i ( E e ) δ [ E e ( E end + E ν i + m lightest ) ] exp ( E e E e ) 2 2 σ 2 ,
d Γ β d E e = 1 2 π σ d E e d Γ β d E e ( E e ) exp ( E e E e ) 2 2 σ 2 ,
Substituting Equation (88) into Equation (108), the smeared spectrum of the emitted electron from the C ν B signal can be written as
d Γ i d E e = N T 2 π σ s ν = ± 1 2 d 3 p ν ( 2 π ) 3 σ ν i ( p ν , s ν ) v ν i f ν i ( p , s ν ) × exp [ E e ( E end + m lightest + E ν i ) ] 2 2 σ 2 ,
where
σ ν i ( p ν , s ν ) = σ ν i ( p ν , s ν , E e ) = σ ν i ( p ν , s ν , E end + m lightest + E ν i ) .
Equation (110) is a Fredholm integral equation of the first kind and d Γ i d E e is a would-be observed quantity. After solving Equation (110) inversely, the spectrum of the C ν B, f ν i ( p , s ν ) , can be in principle reconstructed though we might need a significantly large number of observations for the C ν B events. We leave the detailed study for the reconstruction of the C ν B spectrum f ν i ( p , s ν ) on tritium as future work.
In Figure 5, we show the expected spectra for the emitted electrons from the C ν B signals (solid lines) and the β -decay background (dashed lines) with m lightest = 0 meV and 100 g of tritium, the energy resolution Δ = 20 meV (left panel) and Δ = 0.4 meV (right panel) considering the case of Dirac neutrinos and both the normal (fine red) and inverted (bold blue) mass hierarchies. In these figures, we neglect spectral distortions for the C ν B from e ± -annihilation during their neutrino decoupling and the gravitational clustering for simplicity. We can see that the C ν B signal is distinguished from the β -decay background if Δ E ν i . It is easier to distinguish the C ν B signal from the β -decay background in the inverted mass ordering than the normal ordering. This is because we can obtain a larger number of events for the heaviest neutrinos in the inverted case due to the large value of | U e 1 | . In addition, β -decay spectrum near the endpoint is smaller in the inverted case because in the inverted case the β -decay spectrum near the endpoint is composed of ν 3 with small | U e 3 | while in the normal ordering that is composed of ν 1 with large | U e 1 | .

Comments on Statistical Analysis

To estimate the required energy resolution of the detector Δ and exposure of tritium to discover the C ν B in a qualitative way, we need statistical analysis. In ref. [31], the authors estimated statistical significance for the detection of the C ν B on tritium as a function of the lightest neutrino mass and the energy resolution in an exposure of 100 g yr of tritium using a χ 2 -analysis (see Figure 5 in ref. [31]). Here a fiducial value of constant number events of background of N b = Γ b T , where Γ b = 10 5 Hz in the 15 eV region around the β -decay endpoint energy, is introduced in addition to the β -decay background. If we would obtain a larger exposure of tritium, the result of Figure 5 in ref. [31] will be improved. The reduction of the constant background N b might improve the result. A more quantitative discussion will be possible when the more concrete setup of the PTOLEMY-type experiment is decided, and the neutrino mass ordering and the lightest neutrino mass are constrained more severely from complementary future neutrino experiments.
We leave as future work the statistical analysis to estimate the required energy resolution Δ and exposures to observe the C ν B spectral distortions due to e ± -annihilation in neutrino decoupling and gravitational clustering by nearby galaxies. However, the required energy resolution would not change drastically compared to observing the C ν B itself since their spectral distortions are sub-leading contributions. As discussed in Section 5.2.3, to observe 1–10% modifications in Γ i due to their spectral distortions, one will need 10 2 –10 4 events of the C ν B. The required exposures correspond to 10–10 3 kg yr of the exposure of tritium. It is extremely difficult to achieve this exposure at present. Note that here we consider neutrino masses small enough to satisfy m ν < 0.12 eV . If neutrino masses are enough large, the required exposure will be smaller due to large neutrino clustering. However, it would be difficult to distinguish the C ν B spectral distortions due to e ± -annihilation in neutrino decoupling from such large neutrino clustering experimentally. We also leave as future work how to distinguish the two contributions to the C ν B spectral distortions by numerical simulations and actual experiments.

6. Conclusions

In the near future, CMB-S4 will determine N eff with a very good precision of ∼0.03 at 68 % C.L., and consequently confirm neutrino decoupling process in the SM and/or impose severe constraints on many scenarios in physics beyond the SM. In addition, in the future, a direct observation of the C ν B might bring us more information about the early universe and neutrino physics. In both observations, the C ν B spectrum is one of crucial ingredients to estimate N eff and a direct detection rate.
In this article, we review the formula of kinetic equations for neutrinos in the early universe, which are the quantum Boltzmann equations for neutrinos and the continuity equation and the possible spectral distortions due to e ± -annihilation in neutrino decoupling. We also discuss the impact of the distortion of the C ν B spectrum in neutrino decoupling on direct observation of the C ν B on tritium, with emphasis on the PTOLEMY-type experiment.
We find N eff = 3.044 [25,26,27] by solving the kinetic equations for neutrino density matrix in the early universe, including vacuum three-flavor oscillations, oscillations in e ± -background, finite temperature corrections to m e , ρ and P up to the next-to-leading order O ( e 3 ) (see also ref. [24] for the first suggestion on the importance of this contribution), and the collision term where we consider full diagonal parts and off-diagonal parts derived from charged current interactions but neglect off-diagonal parts derived from neutral current interactions. Later, the authors in refs. [26,27] also find N eff = 3.0440 and 3.0440 ± 0.0002 , respectively, including off-diagonal parts in the collision term derived from neutrino neutral current interactions. Effects of their off-diagonal parts, and the choice of neutrino mass and mixing parameters on N eff are quite small, δ N eff ± ( 1 2 ) × 10 4 [27]. In refs. [25,26,27], the Dirac CP-violating phase in neutrino mixing parameters is neglected. This contribution to N eff is expected to be also quite small since increases and decreases for the energy densities of neutrinos and anti-neutrinos due to the Dirac CP-violating phase would be canceled out (see also ref. [54]). However, QED corrections to weak interaction rates at order O ( e 2 G F 2 ) and forward scattering of neutrinos via their self-interactions have not been precisely taken into account. Recent studies [23,28] suggest that these neglects might still induce uncertainties of ± ( 10 3 10 4 ) in N eff . Although we should consider their contributions to N eff in the future, N eff = 3.044 is still a very good reference value.
We have revealed the spectrum, number and energy density of the C ν B in the current homogeneous and isotropic universe, including the spectral distortions in neutrino decoupling, as in the right panel of Figure 4 and Table 4 and Table 5. Then we have discussed the capture rates of the C ν B on tritium with 1 % precision to observe effects of 1 % enhancement of the number density of the C ν B by the spectral distortions due to e ± -annihilation during neutrino decoupling. Unfortunately, it is extremely difficult to observe such sub-dominant effects since we will need more than 10 kg of tritium. The precise capture rates of the C ν B on tritium will be also useful to distinguish the SM from physics beyond the SM properly.
If observations and theoretical estimations of the C ν B spectrum are improved significantly, we will obtain much richer information about neutrino physics and the early universe. Through direct observations of the C ν B, one can impose significant constraints on neutrino decays and lifetimes in the region of the age of the universe, t 0 = 4.35 × 10 17 s [34,71]. The C ν B spectrum would also have fluctuations imprinted by inflationary perturbations. Towards a precise estimation of anisotropy of the C ν B as the CMB, one would need to solve kinetic equations for neutrinos in an anisotropic background, develop a detection method of the anisotropy, and reduce uncertainties of physical constants such as neutrino mass and mixing parameters, and Newton constant.

Author Contributions

K.A. and M.Y. contributed equally to this work at all stages. All authors have read and agreed to the published version of the manuscript.

Funding

K.A. acknowledges financial support from IBS under the project code, IBS-R018-D1. M.Y. acknowledges financial support from JSPS Grant-in-Aid for Scientific Research No. JP18K18764, JP21H01080, JP21H00069.

Data Availability Statement

Not applicable.

Acknowledgments

We are grateful to Saul Hurwitz for the collaboration in the work [29].

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Kinetic Equations for Neutrinos in Comoving Variables

In this appendix, we write the Boltzmann equations for the neutrino density matrix (58) and the continuity Equation (59) in terms of the comoving variables, x = m e a , y = p a , z = T γ a . In terms of these variables, we can write the Boltzmann Equations (58) as in ref. [20],
d ρ y ( x ) d x = m Pl 3 8 π ρ ¯ i x 2 m e 3 ¯ y ( x ) , ρ y ( x ) + m e 3 x 4 C ¯ [ ρ y ( x ) ] .
where ρ ¯ , ¯ y ( x ) , and C ¯ [ ρ y ( x ) ] are quantities written in the comoving variables, x , y , z . Here we have used the following relations for the Hubble parameter,
H = 1 m Pl 8 π ρ 3 , ρ = ρ ¯ m e x 4 .
The effective Hamiltonian for neutrino oscillations in vacuum and the forward scattering of neutrinos with the e ± , ν , ν ¯ -background (multiplied by m e / x ), ¯ y ( x ) , is given by
¯ y ( x ) = M 2 2 y + 2 G F m e x 4 ( N ¯ e N ¯ e + ) + 2 G F m e x 4 ( N ¯ ν N ¯ ν ¯ ) 2 2 G F y m W 2 m e x 6 ( E ¯ e + P ¯ e + E ¯ e + + P ¯ e + ) 8 2 G F y 3 m Z 2 m e x 6 ( E ¯ ν + E ¯ ν ¯ ) ,
where N ¯ e ± , N ¯ ν , ν ¯ , E ¯ e ± , P ¯ e ± , E ¯ ν , ν ¯ are written in the flavor basis around the temperature of MeV scale as
N ¯ e N ¯ e + = diag ( n ¯ e n ¯ e + , 0 , 0 ) , n e ± = 2 d 3 y ( 2 π ) 3 f e ± ( y ) , N ¯ ν N ¯ ν ¯ = d 3 y ( 2 π ) 3 ρ y ( x ) ρ ¯ y ( x ) , E ¯ e ± + P ¯ e ± = diag ( ρ ¯ e ± + P ¯ e ± , 0 , 0 ) , ρ ¯ e ± + P ¯ e ± = d 3 y ( 2 π ) 3 E ¯ e + y 2 3 E ¯ e f e ± ( y ) , E ¯ ν + E ¯ ν ¯ = d 3 y ( 2 π ) 3 y ρ y ( x ) + ρ ¯ y ( x ) ,
where, neglecting the chemical potential for e ± ,
f e ± ( y ) = 1 e E ¯ e / z + 1 , E ¯ e = y 2 + x 2 + δ m ¯ e 2 .
δ m ¯ e 2 is the finite temperature correction to the electron mass up to O ( e 2 ) in the comoving variables, ignoring the logarithmic term in Equation (57) and the chemical potential for e ± ,
δ m ¯ e 2 = e 2 z 6 + e 2 π 2 d y y 2 y 2 + x 2 1 exp ( y 2 + x 2 / z ) + 1 .
The collision term in the comoving variables can be also decomposed as in Equation (48)
C ¯ [ ρ y ( x ) ] = C ¯ ν ν ¯ e e + + C ¯ ν e ν e + C ¯ ν e + ν e + + C ¯ ν ν ν ν + C ¯ ν ν ¯ ν ν ¯ .
The collision terms from the annihilation and scattering processes including both ν and e ± are, neglecting the chemical potential for e ± and reducing nine-dimensional collision integrals in Equation (37) to two integrals as in Appendix B,
C ¯ ν ν ¯ e e + [ ρ y 1 ( x ) ] = G F 2 2 π 3 y 1 d y 2 d y 3 y 2 y 3 E ¯ 4 × [ Π ann 1 F ann L L ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) + Π ann 2 F ann R R ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) + Π ann 3 ( F ann R L ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) + F ann L R ν ( 1 ) , ν ¯ ( 2 ) , e ( 3 ) , e + ( 4 ) ) ] , C ¯ ν e ν e [ ρ y 1 ( x ) ] + C ¯ ν e + ν e + [ ρ y 1 ( x ) ] = G F 2 2 π 3 y 1 d y 2 d y 3 y 2 y 3 E ¯ 4 × [ Π sc 1 ( F sc L L ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) + F sc R R ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) )
Π sc 2 ( F sc L R ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) + F sc R L ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) ) ] ,
where E ¯ i = y i 2 + x 2 + δ m ¯ e 2 and F sc a b ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) = F sc a b ν ( 1 ) , e + ( 2 ) , ν ( 3 ) , e + ( 4 ) = F sc a b ν ( 1 ) , e ( 2 ) , ν ( 3 ) , e ( 4 ) due to no lepton asymmetry. F ann a b and F sc a b are given by Equations (39) and (42). Similarly, the collision terms from the self-interaction processes in the comoving variables are
C ¯ ν ν ν ν [ ρ y 1 ( x ) ] + C ¯ ν ν ¯ ν ν ¯ [ ρ y 1 ( x ) ] = G F 2 2 π 3 y 1 d y 2 d y 3 y 2 y 3 y 4 × [ Π self 1 F sc ( ν ( 1 ) , ν ( 2 ) , ν ( 3 ) , ν ( 4 ) ) + Π self 2 F sc ( ν ( 1 ) , ν ( 2 ) , ν ( 3 ) , ν ( 4 ) ) + F ann ( ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ ( 4 ) ) ] .
F sc ν ( 1 ) , ν ( 2 ) , ν ( 3 ) , ν ( 4 ) , F sc ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ + ( 4 ) and F ann ν ( 1 ) , ν ¯ ( 2 ) , ν ( 3 ) , ν ¯ ( 4 ) are given by Equations (45)–(47). The functions Π self , ann , sc 1 , 2 , 3 in Equations (A8)–(A9) take the following forms,
Π ann 1 = 2 D 1 2 D 2 ( y 2 , y 3 ) y 2 E ¯ 3 2 D 2 ( y 1 , y 4 ) y 1 E ¯ 4 + 2 D 3 y 1 y 2 E ¯ 3 E ¯ 4 , Π ann 2 = 2 D 1 2 D 2 ( y 2 , y 4 ) y 2 E ¯ 4 2 D 2 ( y 1 , y 3 ) y 1 E ¯ 3 + D 3 y 1 y 2 E ¯ 3 E ¯ 4 , Π ann 3 = ( x 2 + δ m ¯ e 2 ) D 1 + D 2 ( y 1 , y 2 ) y 1 y 2 1 E ¯ 3 E ¯ 4 , Π sc 1 = 4 D 1 2 D 2 ( y 2 , y 3 ) E ¯ 2 y 3 2 D 2 ( y 1 , y 4 ) y 1 E ¯ 4 + 2 D 2 ( y 3 , y 4 ) y 3 E ¯ 4 + 2 D 2 ( y 1 , y 2 ) y 1 E ¯ 2 + 4 D 3 y 1 E ¯ 2 y 3 E ¯ 4 , Π sc 2 = 2 ( x 2 + δ m ¯ e 2 ) D 1 D 2 ( y 1 , y 3 ) y 1 y 3 1 E ¯ 2 E ¯ 4 , Π self 1 = D 1 + D 2 ( y 1 , y 2 ) y 1 y 2 + D 2 ( y 3 , y 4 ) y 3 y 4 + D 3 y 1 y 2 y 3 y 4 , Π self 2 = D 1 D 2 ( y 2 , y 3 ) y 2 y 3 D 2 ( y 1 , y 4 ) y 1 y 4 + D 3 y 1 y 2 y 3 y 4 .
The functions of D 1 , 2 , 3 are written as,
D 1 = 4 π 0 d λ λ 2 sin ( λ y 1 ) sin ( λ y 2 ) sin ( λ y 3 ) sin ( λ y 4 ) , D 2 ( y 3 , y 4 ) = 4 y 3 y 4 π 0 d λ λ 2 sin ( λ y 1 ) sin ( λ y 2 ) cos ( λ y 3 ) sin ( λ y 3 ) λ y 3 cos ( λ y 4 ) sin ( λ y 4 ) λ y 4 , D 3 = 4 y 1 y 2 y 3 y 4 π 0 d λ λ 2 cos ( λ y 1 ) sin ( λ y 1 ) λ y 1 cos ( λ y 2 ) sin ( λ y 2 ) λ y 2 × cos ( λ y 3 ) sin ( λ y 3 ) λ y 3 cos ( λ y 4 ) sin ( λ y 4 ) λ y 4 ,
which can be integrated out analytically as in Appendix B.
If we neglect the off-diagonal components of ρ y ( x ) in the collision terms from neutrino self-interactions, which could have a negligible effect on N eff with 10 3 precision, their collision terms are reduced to
C ¯ ν ν ν ν [ ρ y 1 ( x ) ] + C ¯ ν ν ¯ ν ν ¯ [ ρ y 1 ( x ) ] | diag = G F 2 2 π 3 y 1 d y 2 d y 3 y 2 y 3 y 4 [ 2 Π self 1 + 4 Π self 2 ( ν α ( 1 ) , ν α ( 2 ) , ν α ( 3 ) , ν α ( 4 ) ) + Π self 1 + Π self 2 F ( ν α ( 1 ) , ν β ( 2 ) , ν α ( 3 ) , ν β ( 4 ) ) + Π self 2 F ( ν α ( 1 ) , ν α ( 2 ) , ν β ( 3 ) , ν β ( 4 ) ) + Π self 1 + Π self 2 F ( ν α ( 1 ) , ν γ ( 2 ) , ν α ( 3 ) , ν γ ( 4 ) ) + Π self 2 F ( ν α ( 1 ) , ν α ( 2 ) , ν γ ( 3 ) , ν γ ( 4 ) ) ] .
where
F ( ν α ( 1 ) , ν β ( 2 ) , ν γ ( 3 ) , ν δ ( 4 ) ) = f ν γ ( y 3 ) f ν δ ( y 4 ) 1 f ν α ( y 1 ) 1 f ν β ( y 2 ) f ν α ( y 1 ) f ν β ( y 2 ) 1 f ν γ ( y 3 ) 1 f ν δ ( y 4 ) .
Finally, the continuity Equation (59) is translated into the evolution equation for z, including finite temperature corrections from QED up to O ( e 3 ) but neglecting the logarithmic O ( e 2 ) corrections [14,24],
d z d x = x z J ( x / z ) 1 2 π 2 z 3 0 d y y 3 d f ν e d x + d f ν μ d x + d f ν τ d x + G 1 ( 2 ) ( x / z ) + G 1 ( 3 ) ( x / z ) x 2 z 2 J ( x / z ) + Y ( x / z ) + 2 π 2 15 + G 2 ( 2 ) ( x / z ) + G 2 ( 3 ) ( x / z ) ,
where
G 1 ( 2 ) ( ω ) = 2 π α 1 ω ( K ( ω ) 3 + 2 K ( ω ) 2 J ( ω ) 6 K ( ω ) J ( ω ) + K ( ω ) 6 K ( ω ) K ( ω ) + J ( ω ) 6 + J ( ω ) K ( ω ) + J ( ω ) K ( ω ) ] , G 2 ( 2 ) ( ω ) = 8 π α K ( ω ) 6 + J ( ω ) 6 1 2 K ( ω ) 2 + K ( ω ) J ( ω ) + 2 π α ω K ( ω ) 6 K ( ω ) K ( ω ) + J ( ω ) 6 + J ( ω ) K ( ω ) + J ( ω ) K ( ω ) , G 1 ( 3 ) ( ω ) = e 3 4 π K ( ω ) + ω 2 2 k ( ω ) 1 / 2 [ 1 ω 2 J ( ω ) 4 K ( ω ) 2 J ( ω ) ω 2 j ( ω ) ω 2 k ( ω ) + j ( ω ) 2 J ( ω ) + ω 2 j ( ω ) ω k ( ω ) j ( ω ) + K ( ω ) 2 2 K + ω 2 k ( ω ) ] , G 2 ( 3 ) ( ω ) = e 3 4 π K ( ω ) + ω 2 2 k ( ω ) 1 / 2 ( 2 J ( ω ) + ω 2 j ( ω ) ) 2 2 ( 2 K ( ω ) + ω 2 k ( ω ) ) 2 ω Y ( ω ) ω 3 J ( ω ) + ω 2 j ( ω ) ,
with
K ( ω ) = 1 π 2 0 d u u 2 u 2 + ω 2 1 exp u 2 + ω 2 + 1 , J ( ω ) = 1 π 2 0 d u u 2 exp u 2 + ω 2 exp u 2 + ω 2 + 1 2 , Y ( ω ) = 1 π 2 0 d u u 4 exp u 2 + ω 2 exp u 2 + ω 2 + 1 2 , k ( ω ) = 1 π 2 0 d u 1 u 2 + ω 2 1 exp u 2 + ω 2 + 1 , j ( ω ) = 1 π 2 0 d u exp u 2 + ω 2 exp u 2 + ω 2 + 1 2 .
The prime represents the derivative with respect to ω . G ( 2 ) ( ω ) and G ( 3 ) ( ω ) denote QED finite temperature corrections at O ( e 2 ) and O ( e 3 ) , respectively.

Appendix B. Reduction of the Collision Integrals

In this appendix, we analytically perform seven out of nine integrations in the collision terms for four-Fermi interaction processes at order of O ( G F 2 ) in the homogeneous and isotropic universe, following refs. [9,39]. We consider the general form of the collision term in this case,
C coll = 1 2 E 1 ( 2 π ) 4 δ 4 ( i p i ) | M | 2 part F ρ p i = 2 4 d 3 p i ( 2 π ) 3 2 E i ,
where E i is the energy of i-th particle. The matrix F ρ p is a function of neutrino density matrix and | M | 2 part is a part of the possible squared matrix elements summed over spin degrees of freedom of all particles except for the first particle | M | 2 . We change the delta function for 3-momentum into the exponential representation:
δ ( 3 ) ( i p i ) = e λ · ( p 1 + p 2 p 3 p 4 ) d 3 λ ( 2 π ) 3 ,
and decompose momentum integrations into the radial and angle components,
d 3 p i = p i 2 d p i sin θ i d θ i d ϕ i p i 2 d p i d Ω i .
Using Equations (A19) and (A20), we rewrite the general collision term (A18) to
C coll = 1 64 π 3 E 1 p 1 δ ( E 1 + E 2 E 3 E 4 ) F ( ρ p ( t ) ) D ( p 1 , p 2 , p 3 , p 4 ) p 2 d p 2 E 2 p 3 d p 3 E 3 p 4 d p 4 E 4 ,
where
D ( p 1 , p 2 , p 3 , p 4 ) = p 1 p 2 p 3 p 4 64 π 5 0 λ 2 d λ e i λ · p 1 d Ω λ e i λ · p 2 d Ω p 2 × e i λ · p 3 d Ω p 3 e i λ · p 4 d Ω p 4 | M | 2 .
For four-Fermi interaction processes at order of O ( G F 2 ) , all of | M | 2 have two kinds of forms,
K 1 ( q 1 μ q 2 μ ) ( q 3 ν q 4 ν ) = K 1 ( E 1 E 2 q 1 · q 2 ) ( E 3 E 4 q 3 · q 4 ) ,
K 2 m 2 ( q 3 μ q 4 μ ) = K 2 m 2 ( E 3 E 4 q 3 · q 4 ) ,
where q i corresponds to one of p j and the angle between q i and q j is written in terms of the integration variables of angle,
cos ψ i j = sin θ i sin θ j cos ( ϕ i ϕ j ) + cos θ i cos θ j .
In both cases of Equations (A23) and (A24), we can perform all integrals for angle components in Equation (A22) so that D ( q 1 , q 2 , q 3 , q 4 ) in the case of Equation (A23) reduces to
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 ] ,
while in the case of Equation (A24), D ( q 1 , q 2 , q 3 , q 4 ) is given by
D = K 2 E 1 E 2 [ E 3 E 4 D 1 + D 2 ( q 3 , q 4 ) ] ,
where D 1 , 2 , 3 are defined in Equation (A12).
In the following we only consider D 1 , D 2 ( q 3 , q 4 ) , D 3 . For simplicity we assume that q 1 > q 2 and q 3 > q 4 without loss of generality though we can perform the integrals in D 1 , 2 , 3 without this assumption and obtain the exact expressions given in ref. [39]. Then we obtain the simplified expressions of D 1 , 2 , 3 in four cases:
(1)
q 1 + q 2 > q 3 + q 4 , q 1 + q 4 > q 2 + q 3 and q 1 q 2 + q 3 + q 4
D 1 = 1 2 ( q 2 + q 3 + q 4 q 1 ) , D 2 ( q 3 , q 4 ) = 1 12 ( q 1 q 2 ) 3 + 2 ( q 3 3 + q 4 3 ) 3 ( q 1 q 2 ) ( q 3 2 + q 4 2 ) , D 3 = 1 60 ( q 1 5 5 q 1 3 q 2 2 + 5 q 1 2 q 2 3 q 2 5 5 q 1 3 q 3 2 + 5 q 2 3 q 3 2 + 5 q 1 2 q 3 3 + 5 q 2 2 q 3 3 q 3 5 5 q 1 3 q 4 2 + 5 q 2 3 q 4 2 + 5 q 3 3 q 4 2 + 5 q 1 2 q 4 3 + 5 q 2 2 q 4 3 + 5 q 3 2 q 4 3 q 4 5 ) .
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
D 1 = q 4 , D 2 ( q 3 , q 4 ) = 1 3 q 4 3 , D 3 = 1 30 q 4 3 5 q 1 2 + 5 q 2 2 + 5 q 3 2 q 4 2 .
(3)
q 1 + q 2 < q 3 + q 4 , q 1 + q 4 < q 2 + q 3 and q 3 q 1 + q 2 + q 4
D 1 = 1 2 ( q 1 + q 2 + q 4 q 3 ) , D 2 ( q 3 , q 4 ) = 1 12 ( q 1 + q 2 ) 3 2 q 3 3 + 2 q 4 3 + 3 ( q 1 + q 2 ) ( q 3 3 + q 4 3 ) .
D 3 is equal to that in Equation (A28) with the replacement of variables q 1 q 3 and q 2 q 4 and the case of 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
D 1 = q 2 , D 2 ( q 3 , q 4 ) = 1 6 q 2 3 q 3 2 + 3 q 4 2 3 q 1 2 q 2 2 , D 3 = 1 30 q 2 3 5 q 1 2 + 5 q 3 2 + 5 q 4 2 q 2 2 .
After we have integrated the δ -function, we obtain the simplified expression of the collision term, leaving two integrals,
C coll = 1 64 π 3 E 1 p 1 F ρ p ( t ) D ( p 1 , p 2 , p 3 , p 4 ) p 2 d p 2 E 2 p 3 d p 3 E 3 ,
where E 4 = E 1 + E 2 E 3 and p 4 = E 4 2 m 4 2 .

Appendix C. Kinematics for νi + 3H→e + 3He and 3H→e + 3He + ν ¯ i

In this appendix, we estimate the kinematics of inverse tritium β -decay for the C ν B, ν i + 3 H e + 3 He , and tritium β -decay 3 H e + 3 He + ν ¯ i . We also discuss the kinematic relations between the two processes. In particular, we investigate the maximal energy of the electron emitted from β -decay, called the β -decay endpoint energy, and the energy of the electron emitted from the inverse β -decay process for the C ν B. Here we consider the nuclear process and use the nuclear masses of 3 H and 3 He , m 3 H and m 3 He .
We first consider the kinematics of tritium beta decay, 3 H 3 He + e + ν ¯ i , in the rest frame of 3 H . From 4-momentum conservation, the energy of the electron is
E e = m 3 H 2 + m e 2 m ν i 2 m 3 He 2 2 E ν i E 3 He + 2 | p ν | | p 3 He | cos θ ν 3 He 2 m 3 H .
The maximal energy, E end , is achieved when the emitted anti-neutrino is the lightest and cos θ ν 3 He = 1 ( θ ν 3 He = 0 ) . When the neutrino and the helium-3 nucleus are emitted in parallel, the electron is produced in opposite direction. In addition, the maximization condition of the electron energy corresponds to the minimization condition of ( E ν + E 3 He ) , which yields
E ν i E 3 He = | p ν | | p 3 He | = m ν i m 3 He .
From these conditions, the maximal energy of the electron for 3 H e + 3 He + ν ¯ i is given by
E e max , i = m 3 H 2 + m e 2 ( m ν i + m 3 He ) 2 2 m 3 H .
The endpoint energy of the electron for the tritium β -decay is also given by
E e end = m 3 H 2 + m e 2 ( m lightest + m 3 He ) 2 2 m 3 H .
If the lightest neutrino is massless, the endpoint energy is identified as
E e end , 0 = m 3 H 2 + m e 2 m 3 He 2 2 m 3 H .
Due to m 3 H m 3 He , the difference between the endpoint energy for the massive and massless lightest neutrinos is
E e end E e end , 0 m lightest .
Next we investigate the kinematics of inverse tritium beta decay for relic cosmic neutrinos, ν i + 3 H 3 He + e . In the rest-frame of 3 H , we similarly obtain the energy of the electron as
E e C ν B , i = ( E ν i + m 3 H ) 2 + m e 2 | p ν | 2 + 2 | p ν | | p e | cos θ e ν m 3 He 2 2 ( E ν i + m 3 H ) ( E ν i + m 3 H ) 2 + m e 2 m 3 He 2 2 ( E ν i + m 3 H ) .
where we neglect the terms proportional to | p ν | 2 and | p ν | | p e | and leave the term proportional to E ν i m 3 H because of m 3 H | p e | | p ν | . For m 3 H m e , the difference between E e C ν B , i and E end is
E e C ν B , i E e end E ν i + m lightest .
Since E e C ν B , i E end is (approximately) not function of any nuclear masses, it is insensitive to the uncertainties in the nuclear masses which are calculated from the measured values of atomic masses.

Appendix D. Cross Section for νi + 3H→e + 3He and Decay Rate for 3H→e + 3He + ν ¯ i

In this section, we derive the cross section with 1 % precision for ν i + 3 H e + 3 He , σ ν i , following ref. [34] and the decay rate for 3 H e + 3 He + ν ¯ i , Γ β . We also discuss the spectrum for the tritium β -decay, d Γ β / d E e .

Appendix D.1. Cross Section for νi + 3H→e + 3He

In this section, we follow ref. [34]. The differential cross section for ν i + 3 H e + 3 He takes the following Lorentz invariant form:
d σ ν i d t = 1 16 π | M i | 2 [ s ( m ν i + m 3 H ) 2 ] [ s ( m ν i m 3 H ) 2 ] ,
where s = ( p ν i + p 3 H ) 2 and t = ( p ν i p e ) 2 are the Mandelstam variables, and | M i | 2 is the squared matrix element for the inverse β -decay. In the rest frame of 3 H , s and t are expressed as
s = ( m 3 H + E ν i ) 2 | p ν | 2 = m 3 H 2 + 2 m 3 H E ν i + m ν i 2 , t = ( E e E ν i ) 2 | p e p ν | 2 ( m e m ν i ) 2 + 2 | p e | | p ν | cos θ .
Using also d t / d cos θ = 2 | p e | | p ν | , we obtain
d σ ν i d cos θ = 1 32 π 1 m 3 H 2 | p e | | p ν | | M i | 2 .
The matrix element for ν i + 3 H e + 3 He is effectively given by
i M i = i G F 2 V u d U e i u ¯ e γ μ ( 1 γ 5 ) u ν i u ¯ 3 He γ μ F G γ 5 u 3 H ,
where
F = f F , G = g A 3 g V g G T .
u α denotes the Dirac spinor for species α , g A 1.2723 and g V 1 are the axial and vector coupling constants, respectively, and f F 0.9998 and g G T 3 × ( 0.9511 ± 0.0013 ) denote the reduced matrix elements of the Fermi and Gamow-Teller (GT) operators, respectively, [69].
After averaging over the spins of 3 H and summing over the spins of the outgoing e and 3 He , the squared matrix element is given by
1 2 s e , s 3 H , s 3 He = ± 1 2 | M i | 2 = G F 2 4 | V u d | 2 | U e i | 2 T 1 α β T 2 α β ,
where
T 1 α β = s e = ± 1 / 2 tr γ α ( 1 γ 5 ) u ν i u ¯ ν i γ β ( 1 γ 5 ) u e u ¯ e , T 2 γ δ = s 3 H , s 3 He = ± 1 / 2 tr [ γ γ F G γ 5 u 3 H u ¯ 3 H γ δ F G γ 5 u 3 He u ¯ 3 He ] .
Using the completeness relations, we obtain the relation of Dirac spinors for 3 H , 3 He , and e ,
s j = ± 1 / 2 u j u ¯ j = ( p j + m j ) ,
and for neutrinos with their helicity s ν ,
u ν i u ¯ ν i = 1 2 ( p ν i + m ν i ) ( 1 + 2 s ν γ 5 S ν i ) ,
where S ν i is the spin vector for neutrinos given by
( S ν i ) α = | p ν | m ν i , E ν m ν i p ν | p ν | .
In the massless limit, the previous relation of the Dirac spinor for neutrinos becomes
u ν i u ¯ ν i = 1 2 p ν i 1 2 s γ 5 ,
where we used m S μ = p μ and p μ S μ = 0 . Using the above relations, we rewrite Equation (A47) as
T 1 α β = 1 2 tr γ α ( 1 γ 5 ) ( p ν i + m ν i ) ( 1 + 2 s ν γ 5 S ν i ) γ β ( 1 γ 5 ) ( p e + m e ) ,
T 2 γ δ = tr [ γ γ F G γ 5 p n + m n ) γ δ F G γ 5 ( p p + m p ) ] .
Then we obtain T 1 α β T 2 α β as
T 1 α β T 2 α β = 32 { G + F 2 p e · p 3 He p ν i · p 3 H + G F 2 p e · p 3 H p ν i · p 3 He + G 2 F 2 m 3 H m 3 He p e · p ν i } 64 s ν m ν i { G + F 2 p e · p 3 He S ν i · p 3 H + G F 2 p e · p 3 H S ν i · p 3 He + G 2 F 2 m 3 H m 3 He p e · S ν i } .
In the rest frame of 3 H , neglecting the momentum of 3-helium | p 3 He | / m 3 He ( m 3 H m 3 He ) / m 3 He O ( 10 4 ) , T 1 α β T 2 α β is given by
T 1 α β T 2 α β = 32 m 3 H E 3 He E e E ν i F 2 + 3 G 2 1 2 s ν v ν i + F 2 G 2 v ν i 2 s ν v e cos θ .
We note that θ is the angle between p e and p ν . Finally we obtain the differential cross section for ν i + 3 H e + 3 He , including the enhancement factor due to the Coulombic attraction between e and 3 He , F ( 2 , E e ) , and using also F = f A and G = g A 3 g V g G T ,
d σ ν i d cos θ = G F 2 4 π | V u d | 2 | U e i | 2 F ( 2 , E e ) m 3 He m 3 H v ν i E e | p e | × f A 2 + g A 2 g V 2 g G T 2 ( 1 2 s ν v ν i ) + f A 2 g A 2 3 g V 2 g G T 2 ( v ν i 2 s ν ) v e cos θ .

Appendix D.2. Decay Rate for 3H→e + 3He + ν ¯ i

The decay rate of the β -decay follows the standard formula at the rest frame of tritium,
Γ β = 1 2 9 π 5 m 3 H d 3 p e d 3 p ν i d 3 p 3 He E e E ν i E 3 He | M | 2 δ 4 ( p 3 H p e p ν i p 3 He ) , = 1 2 6 π 4 m 3 H d E e d E ν i | M β | 2 ,
where | M β | 2 is the effective squared matrix element for β -decays summed over spins for the final states and averaged over spins for the initial state,
| M β | 2 = 1 2 i = 1 3 s 3 H , s 3 He , s ν i = ± 1 / 2 | M i | 2 ,
where
i M i = i G F 2 V u d U e i [ u ¯ e γ μ ( 1 γ 5 ) v ν i ] u ¯ 3 H γ μ f F g A 3 g V g G T γ 5 u 3 He .
Then we integrate over E ν i for each E e in Equation (A57). The upper (lower) limit of the integral denotes E ν i max ( E ν i min ) . After some calculations, E ν i max E ν i min and E ν i max + E ν i min are given by
E ν i max E ν i min = 2 m 3 H | p e | M 2 ( E e max , i E e ) 1 / 2 E e max , i E e + 2 m ν i m 3 He m 3 H 1 / 2 , E ν i max + E ν i min = 2 m 3 H M 2 ( m 3 H E e ) E e max , i E e + m ν i m 3 H ( m 3 He + m ν i ) ,
where E e max , i is the maximal energy of the emitted electron for 3 H e + 3 He + ν ¯ i given by Equation (A35) in Appendix C.
M 2 = m 3 H 2 2 m 3 H E e + m e 2 .
Then d Γ β / d E e is given by
d Γ β d E e = 1 2 6 π 3 m 3 H E ν i min E ν i max d E ν i | M β | 2 .
After similar calculations in Appendix D.1 , | M β | 2 for β -decays at rest of tritium is written as
| M β | 2 16 G F 2 | V u d | 2 i = 1 3 | U e i | 2 m 3 H m 3 He E e E ν i × f F 2 + g A 2 g V 2 g G T 2 + f F 2 g A 2 3 g V 2 g G T 2 p ν i · p e E ν i E e ,
where we neglect the momentum of 3 He due to p 3 He m 3 He . In addition, we neglect the second term in Equation (A63) since | p e | m 3 H m 3 He E e . Thus, | M β | 2 approximately becomes
| M β | 2 16 G F 2 | V u d | 2 i = 1 3 | U e i | 2 m 3 H m 3 He E e E ν i f F 2 + g A 2 g V 2 g G T 2 .
Plugging Equation (A64) into Equation (A62), we obtain
d Γ β d E e = G F 2 8 π 3 | V u d | 2 m 3 He E e f F 2 + g A 2 g V 2 g G T 2 × i = 1 3 | U e i | 2 ( E ν i max + E ν i min ) ( E ν i max E ν i min ) .
Finally, substituting Equation (A60) into Equation (A65), we obtain the electron spectrum from the β -decays as
d Γ β d E e = σ ¯ π 2 N T i = 1 3 | U e i | 2 H ( E e , m ν i ) ,
where σ ¯ is the average cross section at the leading order for neutrino capture, including the enhancement due to the Coulombic attraction between e and 3 He , F ( 2 , E e ) ,
σ ¯ = G F 2 2 π | V u d | 2 m 3 He m 3 H f F 2 + g A 2 g V 2 g G T 2 F ( 2 , E e ) E e | p e | .
F ( Z , E e ) is given in Equation (91) and H ( E e , m ν i ) takes the following form,
H ( E e , m ν i ) = 1 E e / m 3 H ( 1 2 E e / m 3 H + m e 2 / m 3 H 2 ) 2 ( E e max , i E e ) E e max , i E e + 2 m ν i m 3 He m 3 H × E e max , i E e + m ν i m 3 H ( m 3 He + m ν i ) .
Then we obtain Γ β ,
Γ β = m e E e end d E e d Γ β d E e ,
where E e end = max { E e max , 1 , E e max , 2 , E e max , 3 } is the endpoint energy of the tritium β -decay given by Equation (A36) in Appendix C.

Notes

1
If we follow the evolution of neutrinos until today, it is also easier to follow the evolution of negative-helicity neutrinos in the mass-diagonal basis since the helicity states of neutrinos are conserved while non-relativistic neutrinos are freely streaming. On the other hand, the chiral states for non-relativistic neutrinos are not conserved.
2
Note that Equation (8) is different from Equation (13) in ref. [19]
3
For forward scattering with background in an anisotropic universe, see ref. [28].
4
The result in the right panel of Figure 4 is quite different from Figure 4 in ref. [20]. Our results are confirmed by Equation (8) and the numerical results in the flavor basis.

References

  1. Aghanim, N. et al. [Planck Collaboration] Planck 2018 results. VI. Cosmological parameters. Astron. Astrophys. 2020, 641, A6, Erratum in Astron. Astrophys. 2021, 652, C4. [Google Scholar] [CrossRef] [Green Version]
  2. Abazajian, K.; Addison, G.; Adshead, P.; Ahmed, Z.; Allen, S.W.; Alonso, D.; Alvarez, M.; Anderson, A.; Arnold, K.S.; Baccigalupi, C.; et al. CMB-S4 Science Case, Reference Design, and Project Plan. arXiv 2019, arXiv:1907.04473. [Google Scholar]
  3. Weinberg, S. Universal Neutrino Degeneracy. Phys. Rev. 1962, 128, 1457–1473. [Google Scholar] [CrossRef]
  4. Cocco, A.G.; Mangano, G.; Messina, M. Probing low energy neutrino backgrounds with neutrino capture on beta decaying nuclei. JCAP 2007, 6, 015. [Google Scholar] [CrossRef] [Green Version]
  5. Dodelson, S.; Turner, M.S. Nonequilibrium neutrino statistical mechanics in the expanding universe. Phys. Rev. D 1992, 46, 3372–3387. [Google Scholar] [CrossRef]
  6. Dolgov, A.D.; Fukugita, M. Nonequilibrium effect of the neutrino distribution on primordial helium synthesis. Phys. Rev. D 1992, 46, 5378–5382. [Google Scholar] [CrossRef]
  7. Fields, B.D.; Dodelson, S.; Turner, M.S. Effect of neutrino heating on primordial nucleosynthesis. Phys. Rev. D 1993, 47, 4309–4314. [Google Scholar] [CrossRef] [Green Version]
  8. Hannestad, S.; Madsen, J. Neutrino decoupling in the early universe. Phys. Rev. D 1995, 52, 1764–1769. [Google Scholar] [CrossRef] [Green Version]
  9. Dolgov, A.D.; Hansen, S.H.; Semikoz, D.V. Nonequilibrium corrections to the spectra of massless neutrinos in the early universe. Nucl. Phys. B 1997, 503, 426–444. [Google Scholar] [CrossRef] [Green Version]
  10. Dolgov, A.D.; Hansen, S.H.; Semikoz, D.V. Nonequilibrium corrections to the spectra of massless neutrinos in the early universe: Addendum. Nucl. Phys. B 1999, 543, 269–274. [Google Scholar] [CrossRef] [Green Version]
  11. Esposito, S.; Miele, G.; Pastor, S.; Peloso, M.; Pisanti, O. Nonequilibrium spectra of degenerate relic neutrinos. Nucl. Phys. B 2000, 590, 539–561. [Google Scholar] [CrossRef]
  12. Hannestad, S. Oscillation effects on neutrino decoupling in the early universe. Phys. Rev. D 2002, 65, 083006. [Google Scholar] [CrossRef] [Green Version]
  13. Fornengo, N.; Kim, C.W.; Song, J. Finite temperature effects on the neutrino decoupling in the early universe. Phys. Rev. D 1997, 56, 5123–5134. [Google Scholar] [CrossRef] [Green Version]
  14. Mangano, G.; Miele, G.; Pastor, S.; Peloso, M. A Precision calculation of the effective number of cosmological neutrinos. Phys. Lett. B 2002, 534, 8–16. [Google Scholar] [CrossRef] [Green Version]
  15. Birrell, J.; Yang, C.T.; Rafelski, J. Relic Neutrino Freeze-out: Dependence on Natural Constants. Nucl. Phys. B 2014, 890, 481–517. [Google Scholar] [CrossRef] [Green Version]
  16. Grohs, E.; Fuller, G.M.; Kishimoto, C.T.; Paris, M.W.; Vlasenko, A. Neutrino energy transport in weak decoupling and big bang nucleosynthesis. Phys. Rev. D 2016, 93, 083522. [Google Scholar] [CrossRef] [Green Version]
  17. Grohs, E.; Fuller, G.M. Insights into neutrino decoupling gleaned from considerations of the role of electron mass. Nucl. Phys. B 2017, 923, 222–244. [Google Scholar] [CrossRef]
  18. Froustey, J.; Pitrou, C. Incomplete neutrino decoupling effect on big bang nucleosynthesis. Phys. Rev. D 2020, 101, 043524. [Google Scholar] [CrossRef] [Green Version]
  19. Mangano, G.; Miele, G.; Pastor, S.; Pinto, T.; Pisanti, O.; Serpico, P.D. Relic neutrino decoupling including flavor oscillations. Nucl. Phys. B 2005, 729, 221–234. [Google Scholar] [CrossRef] [Green Version]
  20. de Salas, P.F.; Pastor, S. Relic neutrino decoupling with flavour oscillations revisited. JCAP 2016, 7, 051. [Google Scholar] [CrossRef]
  21. Gariazzo, S.; de Salas, P.F.; Pastor, S. Thermalisation of sterile neutrinos in the early Universe in the 3+1 scheme with full mixing matrix. JCAP 2019, 7, 014. [Google Scholar] [CrossRef]
  22. Escudero, M. Neutrino decoupling beyond the Standard Model: CMB constraints on the Dark Matter mass with a fast and precise Neff evaluation. JCAP 2019, 2, 007. [Google Scholar] [CrossRef] [Green Version]
  23. Escudero Abenza, M. Precision early universe thermodynamics made simple: Neff and neutrino decoupling in the Standard Model and beyond. JCAP 2020, 5, 048. [Google Scholar] [CrossRef]
  24. Bennett, J.J.; Buldgen, G.; Drewes, M.; Wong, Y.Y.Y. Towards a precision calculation of the effective number of neutrinos Neff in the Standard Model I: The QED equation of state. JCAP 2020, 3, 3. [Google Scholar] [CrossRef]
  25. Akita, K.; Yamaguchi, M. A precision calculation of relic neutrino decoupling. JCAP 2020, 8, 012. [Google Scholar] [CrossRef]
  26. Froustey, J.; Pitrou, C.; Volpe, M.C. Neutrino decoupling including flavour oscillations and primordial nucleosynthesis. JCAP 2020, 12, 015. [Google Scholar] [CrossRef]
  27. Bennett, J.J.; Buldgen, G.; De Salas, P.F.; Drewes, M.; Gariazzo, S.; Pastor, S.; Wong, Y.Y.Y. Towards a precision calculation of Neff in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED. JCAP 2021, 4, 073. [Google Scholar] [CrossRef]
  28. Hansen, R.S.L.; Shalgar, S.; Tamborra, I. Neutrino flavor mixing breaks isotropy in the early universe. JCAP 2021, 7, 017. [Google Scholar] [CrossRef]
  29. Akita, K.; Hurwitz, S.; Yamaguchi, M. Precise Capture Rates of Cosmic Neutrinos and Their Implications on Cosmology. Eur. Phys. J. C 2021, 81, 344. [Google Scholar] [CrossRef]
  30. Betti, M.G.; Biasotti, M.; Boscá, A.; Calle, F.; Carabe-Lopez, J.; Cavoto, G.; Chang, C.; Chung, W.; Cocco, A.G.; Colijn, A.P.; et al. A design for an electromagnetic filter for precision energy measurements at the tritium endpoint. Prog. Part. Nucl. Phys. 2019, 106, 120–131. [Google Scholar] [CrossRef] [Green Version]
  31. Betti, M.G.; Biasotti, M.; Boscá, A.; Calle, F.; Canci, N.; Cavoto, G.; Chang, C.; Cocco, A.G.; Colijn, A.P.; Conrad, J.; et al. Neutrino physics with the PTOLEMY project: Active neutrino properties and the light sterile case. JCAP 2019, 7, 047. [Google Scholar] [CrossRef]
  32. Lazauskas, R.; Vogel, P.; Volpe, C. Charged current cross section for massive cosmological neutrinos impinging on radioactive nuclei. J. Phys. G 2008, 35, 025001. [Google Scholar] [CrossRef]
  33. Blennow, M. Prospects for cosmic neutrino detection um experiments in the case of hierarchical neutrino masses. Phys. Rev. D 2008, 77, 113014. [Google Scholar] [CrossRef] [Green Version]
  34. Long, A.J.; Lunardini, C.; Sabancilar, E. Detecting non-relativistic cosmic neutrinos by capture on tritium: Phenomenology and physics potential. JCAP 2014, 8, 038. [Google Scholar] [CrossRef] [Green Version]
  35. Roulet, E.; Vissani, F. On the capture rates of big bang neutrinos by nuclei within the Dirac and Majorana hypotheses. JCAP 2018, 10, 049. [Google Scholar] [CrossRef] [Green Version]
  36. Mertsch, P.; Parimbelli, G.; de Salas, P.F.; Gariazzo, S.; Lesgourgues, J.; Pastor, S. Neutrino clustering in the Milky Way and beyond. JCAP 2020, 1, 015. [Google Scholar] [CrossRef] [Green Version]
  37. Sigl, G.; Raffelt, G. General kinetic description of relativistic mixed neutrinos. Nucl. Phys. B 1993, 406, 423–451. [Google Scholar] [CrossRef]
  38. Notzold, D.; Raffelt, G. Neutrino Dispersion at Finite Temperature and Density. Nucl. Phys. B 1988, 307, 924–936. [Google Scholar] [CrossRef]
  39. Blaschke, D.N.; Cirigliano, V. Neutrino Quantum Kinetic Equations: The Collision Term. Phys. Rev. D 2016, 94, 033009. [Google Scholar] [CrossRef] [Green Version]
  40. Froustey, J. The Universe at the MeV Era: Neutrino Evolution and Cosmological Observables. Ph.D. Thesis, Paris Institution of Astrophysis, Paris, France, 2022. [Google Scholar]
  41. Kapusta, J.I.; Gale, C. Finite-Temperature Field Theory: Principles and Applications; Cambridge Monographs on Mathematical Physics; Cambridge University Press: Cambridge, UK, 2011. [Google Scholar] [CrossRef]
  42. Heckler, A.F. Astrophysical applications of quantum corrections to the equation of state of a plasma. Phys. Rev. D 1994, 49, 611–617. [Google Scholar] [CrossRef]
  43. Lopez, R.E.; Turner, M.S. An Accurate Calculation of the Big Bang Prediction for the Abundance of Primordial Helium. Phys. Rev. D 1999, 59, 103502. [Google Scholar] [CrossRef] [Green Version]
  44. Dolgov, A.D.; Hansen, S.H.; Pastor, S.; Petcov, S.T.; Raffelt, G.G.; Semikoz, D.V. Cosmological bounds on neutrino degeneracy improved by flavor oscillations. Nucl. Phys. B 2002, 632, 363–382. [Google Scholar] [CrossRef] [Green Version]
  45. Wong, Y.Y.Y. Analytical treatment of neutrino asymmetry equilibration from flavor oscillations in the early universe. Phys. Rev. D 2002, 66, 025015. [Google Scholar] [CrossRef] [Green Version]
  46. Abazajian, K.N.; Beacom, J.F.; Bell, N.F. Stringent Constraints on Cosmological Neutrino Antineutrino Asymmetries from Synchronized Flavor Transformation. Phys. Rev. D 2002, 66, 013008. [Google Scholar] [CrossRef] [Green Version]
  47. Mangano, G.; Miele, G.; Pastor, S.; Pisanti, O.; Sarikas, S. Updated BBN bounds on the cosmological lepton asymmetry for non-zero θ13. Phys. Lett. B 2012, 708, 1–5. [Google Scholar] [CrossRef] [Green Version]
  48. Castorina, E.; Franca, U.; Lattanzi, M.; Lesgourgues, J.; Mangano, G.; Melchiorri, A.; Pastor, S. Cosmological lepton asymmetry with a nonzero mixing angle θ13. Phys. Rev. D 2012, 86, 023517. [Google Scholar] [CrossRef] [Green Version]
  49. Oldengott, I.M.; Schwarz, D.J. Improved constraints on lepton asymmetry from the cosmic microwave background. EPL 2017, 119, 29001. [Google Scholar] [CrossRef] [Green Version]
  50. Escudero, M.; Ibarra, A.; Maura, V. Primordial Lepton Asymmetries in the Precision Cosmology Era: Current Status and Future Sensitivities from BBN and the CMB. arXiv 2022, arXiv:2208.03201. [Google Scholar]
  51. Esteban, I.; Gonzalez-Garcia, M.C.; Maltoni, M.; Schwetz, T.; Zhou, A. The fate of hints: Updated global analysis of three-flavor neutrino oscillations. JHEP 2020, 9, 178. [Google Scholar] [CrossRef]
  52. de Salas, P.F.; Forero, D.V.; Gariazzo, S.; Martínez-Miravé, P.; Mena, O.; Ternes, C.A.; Tórtola, M.; Valle, J.W.F. 2020 global reassessment of the neutrino oscillation picture. JHEP 2021, 2, 071. [Google Scholar] [CrossRef]
  53. Esteban, I.; Gonzalez-Garcia, M.C.; Hernandez-Cabezudo, A.; Maltoni, M.; Schwetz, T. Global analysis of three-flavour neutrino oscillations: Synergies and tensions in the determination of θ23, δCP, and the mass ordering. JHEP 2019, 1, 106. [Google Scholar] [CrossRef]
  54. Froustey, J.; Pitrou, C. Primordial neutrino asymmetry evolution with full mean-field effects and collisions. JCAP 2022, 3, 065. [Google Scholar] [CrossRef]
  55. Esposito, S.; Mangano, G.; Miele, G.; Picardi, I.; Pisanti, O. Neutrino energy loss rate in a stellar plasma. Nucl. Phys. B 2003, 658, 217–253. [Google Scholar] [CrossRef] [Green Version]
  56. Fixsen, D.J. The Temperature of the Cosmic Microwave Background. Astrophys. J. 2009, 707, 916–920. [Google Scholar] [CrossRef] [Green Version]
  57. Singh, S.; Ma, C.P. Neutrino clustering in cold dark matter halos: Implications for ultrahigh-energy cosmic rays. Phys. Rev. D 2003, 67, 023506. [Google Scholar] [CrossRef] [Green Version]
  58. Ringwald, A.; Wong, Y.Y.Y. Gravitational clustering of relic neutrinos and implications for their detection. JCAP 2004, 12, 005. [Google Scholar] [CrossRef] [Green Version]
  59. de Salas, P.F.; Gariazzo, S.; Lesgourgues, J.; Pastor, S. Calculation of the local density of relic neutrinos. JCAP 2017, 9, 034. [Google Scholar] [CrossRef] [Green Version]
  60. Zhang, J.; Zhang, X. Gravitational clustering of cosmic relic neutrinos in the Milky Way. Nat. Commun. 2018, 9, 1833. [Google Scholar] [CrossRef] [Green Version]
  61. Alvey, J.; Escudero, M.; Sabti, N.; Schwetz, T. Cosmic neutrino background detection in large-neutrino-mass cosmologies. Phys. Rev. D 2022, 105, 063501. [Google Scholar] [CrossRef]
  62. Duda, G.; Gelmini, G.; Nussinov, S. Expected signals in relic neutrino detectors. Phys. Rev. D 2001, 64, 122001. [Google Scholar] [CrossRef] [Green Version]
  63. Safdi, B.R.; Lisanti, M.; Spitz, J.; Formaggio, J.A. Annual Modulation of Cosmic Relic Neutrinos. Phys. Rev. D 2014, 90, 043001. [Google Scholar] [CrossRef]
  64. Primakoff, H.; Rosen, S.P. Double beta decay. Rept. Prog. Phys. 1959, 22, 121–166. [Google Scholar] [CrossRef]
  65. Aker, M.; Beglarian, A.; Behrens, J.; Berlev, A.; Besserer, U.; Bieringer, B.; Block, F.; Bobien, S.; Boettcher, M.; Bornschein, B.; et al. Direct neutrino-mass measurement with sub-electronvolt sensitivity. Nat. Phys. 2022, 18, 160–166. [Google Scholar] [CrossRef]
  66. Abe, K. et al. [T2K Collaboration] Proposal for an Extended Run of T2K to 20×1021 POT. arXiv 2016, arXiv:1609.04111. [Google Scholar]
  67. Abe, K. et al. [Hyper-Kamiokande Collaboration] Hyper-Kamiokande Design Report. arXiv 2018, arXiv:1805.04163. [Google Scholar]
  68. Abi, B. et al. [DUNE Collaboration] Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume II: DUNE Physics. arXiv 2020, arXiv:2002.03005. [Google Scholar]
  69. Baroni, A.; Girlanda, L.; Kievsky, A.; Marcucci, L.E.; Schiavilla, R.; Viviani, M. Tritium β-decay in chiral effective field theory. Phys. Rev. C 2016, 94, 024003, Erratum in Phys. Rev. C 2017, 95, 059902. [Google Scholar] [CrossRef] [Green Version]
  70. Masood, S.S.; Nasri, S.; Schechter, J.; Tortola, M.A.; Valle, J.W.F.; Weinheimer, C. Exact relativistic beta decay endpoint spectrum. Phys. Rev. C 2007, 76, 045501. [Google Scholar] [CrossRef] [Green Version]
  71. Akita, K.; Lambiase, G.; Yamaguchi, M. Unstable cosmic neutrino capture. JHEP 2022, 2, 132. [Google Scholar] [CrossRef]
Figure 1. One-loop thermal contributions to forward scattering of neutrinos in the flavor basis with α , β = e ± , μ ± and τ ± . (Left) Tadpole diagram with all flavors in the one-loop. (Right) Babble diagram with the same flavor in the one-loop.
Figure 1. One-loop thermal contributions to forward scattering of neutrinos in the flavor basis with α , β = e ± , μ ± and τ ± . (Left) Tadpole diagram with all flavors in the one-loop. (Right) Babble diagram with the same flavor in the one-loop.
Universe 08 00552 g001
Figure 2. (Left panel): 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 with QED finite temperature corrections up to O ( e 3 ) . (Right panel): Final distortions of flavor neutrino spectra as a function of the comoving momentum y with QED finite temperature corrections up to O ( e 3 ) . Upper (lower) dotted line is for ν e ( ν μ , τ ) without neutrino oscillations, while inner solid and dashed lines represent those for flavor neutrinos with neutrino oscillations.
Figure 2. (Left panel): 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 with QED finite temperature corrections up to O ( e 3 ) . (Right panel): Final distortions of flavor neutrino spectra as a function of the comoving momentum y with QED finite temperature corrections up to O ( e 3 ) . Upper (lower) dotted line is for ν e ( ν μ , τ ) without neutrino oscillations, while inner solid and dashed lines represent those for flavor neutrinos with neutrino oscillations.
Universe 08 00552 g002
Figure 3. Feynman diagrams that contribute the weak interaction rates up to O ( e 2 G F 2 ) [27,55]. (0): 4-Fermi interactions. QED finite temperature corrections (i): additional photon emissions and absorptions, (ii): corrections to the dispersion relation for e ± , (iii): vertex corrections, (iv): corrections mediated by photon propagator. Matrix elements multiplied by (0) and one of (iiiv), and squared matrix elements for (i) contribute the weak interaction rates at O ( e 2 G F 2 ) .
Figure 3. Feynman diagrams that contribute the weak interaction rates up to O ( e 2 G F 2 ) [27,55]. (0): 4-Fermi interactions. QED finite temperature corrections (i): additional photon emissions and absorptions, (ii): corrections to the dispersion relation for e ± , (iii): vertex corrections, (iv): corrections mediated by photon propagator. Matrix elements multiplied by (0) and one of (iiiv), and squared matrix elements for (i) contribute the weak interaction rates at O ( e 2 G F 2 ) .
Universe 08 00552 g003
Figure 4. (Left panel): Time evolution of the distortions of neutrinos in the mass basis for a fixed momentum ( y = 5 ) with QED finite temperature corrections up to O ( e 3 ) . (Right panel): Final distortions of neutrino spectra in the mass basis as a function of the comoving momentum y with QED finite temperature corrections up to O ( e 3 ) .
Figure 4. (Left panel): Time evolution of the distortions of neutrinos in the mass basis for a fixed momentum ( y = 5 ) with QED finite temperature corrections up to O ( e 3 ) . (Right panel): Final distortions of neutrino spectra in the mass basis as a function of the comoving momentum y with QED finite temperature corrections up to O ( e 3 ) .
Universe 08 00552 g004
Figure 5. The expected spectra as a function of the electron kinetic energy, K e = E e m e , for the emitted electrons from the C ν B signals (solid lines) and the β -decay background (dashed lines) in a tritium experiment, assuming m lightest = 0 meV and 100 g of tritium, with the energy resolution Δ = 20 meV (left panel) and Δ = 0.4 meV (right panel) in the case of Dirac neutrinos. Bold blue lines represent the NH case and fine red lines represent the IH case.
Figure 5. The expected spectra as a function of the electron kinetic energy, K e = E e m e , for the emitted electrons from the C ν B signals (solid lines) and the β -decay background (dashed lines) in a tritium experiment, assuming m lightest = 0 meV and 100 g of tritium, with the energy resolution Δ = 20 meV (left panel) and Δ = 0.4 meV (right panel) in the case of Dirac neutrinos. Bold blue lines represent the NH case and fine red lines represent the IH case.
Universe 08 00552 g005
Table 1. Squared matrix elements with the symmetric factor S | M | 2 for processes ν e ( p 1 ) + b ( p 2 ) c ( p 3 ) + d ( p 4 ) . g L = 1 2 + sin 2 θ W 2 and g R = sin 2 θ W correspond ( Y L ) 11 and ( Y R ) 11 in Equation (33). For processes of ν μ and ν τ , ν μ ( τ ) ( p 1 ) + b ( p 2 ) c ( p 3 ) + d ( p 4 ) , squared matrix elements are obtained by the substitutions of g L g L 1 = 1 2 + sin 2 θ W , which corresponds ( Y L ) 22 ( 33 ) in Equation (33) [9].
Table 1. Squared matrix elements with the symmetric factor S | M | 2 for processes ν e ( p 1 ) + b ( p 2 ) c ( p 3 ) + d ( p 4 ) . g L = 1 2 + sin 2 θ W 2 and g R = sin 2 θ W correspond ( Y L ) 11 and ( Y R ) 11 in Equation (33). For processes of ν μ and ν τ , ν μ ( τ ) ( p 1 ) + b ( p 2 ) c ( p 3 ) + d ( p 4 ) , squared matrix elements are obtained by the substitutions of g L g L 1 = 1 2 + sin 2 θ W , which corresponds ( Y L ) 22 ( 33 ) in Equation (33) [9].
Process 2 5 G F 2 S | M | 2
ν e + ν ¯ e ν e + ν e ¯ 4 ( p 1 · p 4 ) ( p 2 · p 3 )
ν e + ν e ν e + ν e 2 ( p 1 · p 2 ) ( p 3 · p 4 )
ν e + ν ¯ e ν μ ( τ ) + ν ¯ μ ( τ ) ( p 1 · p 4 ) ( p 2 · p 3 )
ν e + ν ¯ μ ( τ ) ν e + ν ¯ μ ( τ ) ( p 1 · p 4 ) ( p 2 · p 3 )
ν e + ν μ ( τ ) ν e + ν μ ( τ ) ( p 1 · p 2 ) ( p 3 · p 4 )
ν e + ν ¯ e e + e + 4 [ g L 2 ( p 1 · p 4 ) ( p 2 · p 3 ) + g R 2 ( p 1 · p 3 ) ( p 2 · p 4 ) + g L g R m e 2 ( p 1 · p 2 ) ]
ν e + e ν e + e 4 [ g L 2 ( p 1 · p 2 ) ( p 3 · p 4 ) + g R 2 ( p 1 · p 4 ) ( p 2 · p 3 ) g L g R m e 2 ( p 1 · p 3 ) ]
ν e + e + ν e + e + 4 [ g R 2 ( p 1 · p 2 ) ( p 3 · p 4 ) + g L 2 ( p 1 · p 4 ) ( p 2 · p 3 ) g L g R m e 2 ( p 1 · p 3 ) ]
Table 2. Final values of comoving photon temperature and the effective number of neutrinos for flavor neutrinos in several cases.
Table 2. Final values of comoving photon temperature and the effective number of neutrinos for flavor neutrinos in several cases.
Case z fin N eff
Instantaneous decoupling1.401023.00000
No mixing + No QED1.399103.03404
No mixing + QED up to O ( e 2 ) 1.397893.04430
No mixing + QED up to O ( e 3 ) 1.398003.04335
mixing + QED up to O ( e 2 ) 1.397863.04486
mixing + QED up to O ( e 3 ) 1.397973.04391
Table 3. Final values of the distortions of energy densities δ ρ ¯ ν α ( ρ ν α ρ ν eq ) / ρ ν eq and number densities δ n ¯ ν α ( n ν α n ν eq ) / n ν eq for flavor neutrinos in several cases.
Table 3. Final values of the distortions of energy densities δ ρ ¯ ν α ( ρ ν α ρ ν eq ) / ρ ν eq and number densities δ n ¯ ν α ( n ν α n ν eq ) / n ν eq for flavor neutrinos in several cases.
Case δ ρ ¯ ν e ( % ) δ ρ ¯ ν μ ( % ) δ ρ ¯ ν τ ( % ) δ n ¯ ν e ( % ) δ n ¯ ν μ ( % ) δ n ¯ ν τ ( % )
Instantaneous decoupling000000
No mixing + No QED0.9490.3970.3970.5830.2400.240
No mixing + QED up to O ( e 2 ) 0.9370.3910.3910.5750.2360.236
No mixing + QED up to O ( e 3 ) 0.9370.3910.3910.5750.2360.236
mixing + QED up to O ( e 2 ) 0.7120.5110.5230.4350.3110.319
mixing + QED up to O ( e 3 ) 0.7120.5110.5230.4360.3120.319
Table 4. Final values of the distortions of “relativistic” energy densities δ ρ ¯ ν i δ ρ ν i / ρ ν 0 and number densities δ n ¯ ν i ( n ν i n ν 0 ) / n ν 0 for neutrinos in the mass basis after neutrino decoupling.
Table 4. Final values of the distortions of “relativistic” energy densities δ ρ ¯ ν i δ ρ ν i / ρ ν 0 and number densities δ n ¯ ν i ( n ν i n ν 0 ) / n ν 0 for neutrinos in the mass basis after neutrino decoupling.
z fin δ ρ ¯ ν 1 ( % ) δ ρ ¯ ν 2 ( % ) δ ρ ¯ ν 3 ( % ) δ n ¯ ν 1 ( % ) δ n ¯ ν 2 ( % ) δ n ¯ ν 3 ( % )
1.397970.7640.5740.4090.4680.3500.248
Table 5. Neutrino number density per one degree of freedom in the current homogeneous and isotropic universe including non-thermal distortions due to e ± -annihilation during neutrino decoupling.
Table 5. Neutrino number density per one degree of freedom in the current homogeneous and isotropic universe including non-thermal distortions due to e ± -annihilation during neutrino decoupling.
n ν 1 ( cm 3 ) n ν 2 ( cm 3 ) n ν 3 ( cm 3 )
56.6456.5756.52
Table 6. Energy density per one degree of freedom for the lightest neutrinos with m ν lightest = 0 in the current homogeneous and isotropic universe including non-thermal distortions due to e ± -annihilation during neutrino decoupling.
Table 6. Energy density per one degree of freedom for the lightest neutrinos with m ν lightest = 0 in the current homogeneous and isotropic universe including non-thermal distortions due to e ± -annihilation during neutrino decoupling.
Case ρ ν lightest ( meV cm 3 )
Normal Ordering ( ν lightest = ν 1 , m ν 1 = 0 )30.08
Inverted Ordering ( ν lightest = ν 3 , m ν 3 = 0 )29.97
Table 7. Deviation of relic neutrino number density including non-thermal distortions during neutrino decoupling from the case when neutrinos decoupled instantaneously and all e ± -pairs annihilated into photons.
Table 7. Deviation of relic neutrino number density including non-thermal distortions during neutrino decoupling from the case when neutrinos decoupled instantaneously and all e ± -pairs annihilated into photons.
δ n ν 1 d (%) δ n ν 2 d (%) δ n ν 3 d (%)
1.131.010.91
Table 8. The enhancement factor, δ n ν i c , due to neutrino clustering by our Galaxy and nearby galaxies for given values of neutrino masses [36].
Table 8. The enhancement factor, δ n ν i c , due to neutrino clustering by our Galaxy and nearby galaxies for given values of neutrino masses [36].
m ν i (meV) δ n ν i c ( % )
100.53
5012
10050
200300
Table 9. Capture rates of relic cosmic neutrinos on 100 g of tritium in unit of year 1 with m lightest = 0 . δ Γ i d is the differences between the cases with and without effects of e ± -annihilation during neutrino decoupling and δ Γ i c is the differences with and without gravitational clustering for relic neutrinos in nearby galaxies.
Table 9. Capture rates of relic cosmic neutrinos on 100 g of tritium in unit of year 1 with m lightest = 0 . δ Γ i d is the differences between the cases with and without effects of e ± -annihilation during neutrino decoupling and δ Γ i c is the differences with and without gravitational clustering for relic neutrinos in nearby galaxies.
OrderingCase Γ 1 δ Γ 1 d δ Γ 1 c Γ 2 δ Γ 2 d δ Γ 2 c Γ 3 δ Γ 3 d δ Γ 3 c
NOMajorana5.480.06102.400.0240.0130.200 1.6 × 10 3 0.021
Dirac5.480.06101.270.012 6.3 × 10 3 0.101 8.0 × 10 4 0.011
IOMajorana6.130.0610.652.670.0240.280.178 1.6 × 10 3 0
Dirac3.100.0310.331.350.0120.140.178 1.6 × 10 3 0
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Akita, K.; Yamaguchi, M. A Review of Neutrino Decoupling from the Early Universe to the Current Universe. Universe 2022, 8, 552. https://doi.org/10.3390/universe8110552

AMA Style

Akita K, Yamaguchi M. A Review of Neutrino Decoupling from the Early Universe to the Current Universe. Universe. 2022; 8(11):552. https://doi.org/10.3390/universe8110552

Chicago/Turabian Style

Akita, Kensuke, and Masahide Yamaguchi. 2022. "A Review of Neutrino Decoupling from the Early Universe to the Current Universe" Universe 8, no. 11: 552. https://doi.org/10.3390/universe8110552

APA Style

Akita, K., & Yamaguchi, M. (2022). A Review of Neutrino Decoupling from the Early Universe to the Current Universe. Universe, 8(11), 552. https://doi.org/10.3390/universe8110552

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop