A review of neutrino decoupling from the early universe to the current universe

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 $\mathcal{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$\nu$B) in the current universe, and discuss the implications of these effects on the capture rates in direct detection of the C$\nu$B on tritium, with emphasis on the PTOLEMY-type experiment. In addition, we find a precise value of the effective number of neutrinos, $N_{\rm eff}=3.044$. However, QED corrections to weak interaction rates at order $\mathcal{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 $\pm(10^{-3} - 10^{-4})$ in $N_{\rm eff}$.


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 ≡ 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 of 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 equa-tions 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 2 F ) 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 2 21 , |∆m 2 31 |, two massive neutrinos at least are nonrelativistic. 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.

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.

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 where a α (p, t) = e iHt a α (p)e −iHt and b α (p, t) = e −iHt b α (p)e iHt 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 = 2p 0 . The annihilation and creation operators satisfy the anti-commutation relations, For freely evolving massless neutrinos without any interactions, a 0 α (p, t) = a α (p)e −ip 0 t and b 0 α (p, t) = b α (p)e −ip 0 t and the Dirac spinors satisfy free Dirac equations, / pu 0 p = 0, / pv 0 p = 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, 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, p ) denotes the Dirac spinor for negative-helicity particles (positive-helicity anti-particles), which is also normalized to be u 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, where U αi represents a component of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix U PMNS . Due to eq. (2.5), in the ultra-relativistic limit, the relation of the density matrices in the flavor and the mass bases is described as In addition, after neutrino decoupling, the off-diagonal parts of the density matrix in the mass basis are zero, (ρ p ) ij 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 Note that eq. (2.7) is only valid when neutrinos are relativistic and decoupled with thermal plasma. Our numerical calculations also confirm eq. (2.7).

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.

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, where [·, ·] represents the commutator of matrices with a flavor (or mass) index and N αβ is the neutrino density operator, H is the full Hamiltonian in a system, which can be separated into 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 , where t 0 is the initial time and H 0 int 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 The first order solution (2.11) 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 eq. (2.11) into eq. (2.8), and an analogous equation for anti-neutrinos [37], , which is not solved in this article since we assume no lepton asymmetry. Here H 0 free 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 0 free , N 0 αβ (t) and H 0 int , N 0 αβ (t) . Note that the differential equation (2.13) is not closed for both N αβ and N 0 αβ . To close and simplify the differential equation (2.13), 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 eq. (2.13) can be rewritten as (2.14) Thus, the time evolution of the expectation value of N 0 αβ (0) concerning the initial state, ρ p (0), is given by eq. (2.15) 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 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.

Neutrino oscillation in vacuum
The calculation of the first term in the RHS of eq. (2.15) is well established in the mass basis. The free Hamiltonian of neutrinos in the mass basis is given by where γ = (γ 1 , γ 2 , γ 3 ) are the gamma matrices. After substituting the free operators for left-handed neutrinos, the free Hamiltonian becomes The first term in the RHS of eq. (2.15) in the mass basis is written as where M 2 diag = diag(m 2 ν 1 , m 2 ν 2 , m 2 ν 3 ) 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 eq. (2.15) 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 eq. (2.15) , following the similar procedure in the mass basis, where M 2 = U PMNS M 2 diag U † PMNS is the neutrino mass matrix in the flavor basis. For antineutrinos, the corresponding term is obtained by adding a minus sign for the reverse indices in the anti-neutrino density matrix

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 where D Z µν (p) and D W µν (p) are the full propagator of Z 0 boson and W ± boson, 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 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 CC ∝ J µ CC . For the charged current interactions, the second term in the RHS of eq. (2.15), which represents forward scattering of neutrinos with e ± -background, is given by [37,38] i H 0 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 (2.27) where E e = p 2 + m 2 e . 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 eq. (2.27).
For the neutral current interactions, the second term in the RHS of eq. (2.15), which represents forward scattering of neutrinos via neutrino self-interactions, is given by [37,38] 28) 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 where we neglect neutrino masses since neutrinos are relativistic in neutrino decoupling. 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 antineutrinos follow the same evolution, ρ p =ρ T p =ρ p , in the case of no lepton asymmetry.

Collision term
Finally we discuss the third term in the RHS of eq. (2.15) 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 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 eν int + H ν int , while we ignore the part including only electrons and positrons, Here we have used the following Fierz transformation in the charged currents, where H ab↔cd 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 eq. (2.15), The other terms also denote forward scattering, which would give tiny modifications of eqs. (2.26) and (2.28). The first term in eq. (2.35) 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 antineutrinos.
In a schematic manner, the collision term for two-body reactions 1 + 2 ↔ 3 + 4 at tree level takes the following expressions, where ρ i (i = 1, 2, 3, 4) denote the neutrino density matrix, not the energy density, and , 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 eq. (2.36) can be reduced analytically to two integrals as in appendix B.
In the following, we rewrite the collision terms C[ρ p (t)] including eq. (2.35) 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].  (33) in eq. (2.32) [9].
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 0 We can calculate the corresponding collision terms, which are denoted as (2π Here f e ± (p) is the distribution function for electrons and positrons, respectively.

Continuity equation
In addition to the Boltzmann equations for the neutrino density matrix, the energy conservation law must be satisfied, 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, (2.49) The Hubble parameter in eq. (2.48) is calculated using the usual relation, 3H 2 m 2 Pl = 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.

Finite temperature QED corrections to m e , ρ and P up to
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 (2.59) at order O(e 2 G 2 F ) 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 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] where E p = p 2 + m 2 e and N F (p) is the sum of the distribution functions for e ± , The next-to-leading order of thermal corrections to ρ, P is O(e 3 ), not O(e 4 ). These nontrivial 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], Finally, we read the total energy density and the total pressure of electromagnetic plasma up to O(e 3 ) corrections as The thermal corrections to the e ± mass at O(e 2 ) is given by, through modifications of the e ± self energy [42], The last logarithmic terms in eqs. (2.51) and (2.56) 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.

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 eqs. (2.57)-(2.60) 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] 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, H 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 eq. (2.58) 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 eq. (2.55) (see also eqs. (2.49), (2.51) and (2.53) for the detail of eq. (2.55)).
The effective Hamiltonian for the neutrino oscillations in vacuum and the forward scattering of neutrinos in the e ± , ν,ν-background is given by 3 where G F is the Fermi coupling constant and m W , m Z are the W and Z boson masses, respectively. The first term in the RHS of eq. (2.59) denotes neutrino oscillations in vacuum and M 2 is the mass-squared matrix. In the flavor basis, we can write . 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 , where E e = p 2 + m 2 e + δm 2 e (p, T ) and f e ± (p) is the distribution function of e ± . δm 2 e (p, T ) is the QED finite temperature correction to m e , which is given by eq. (2.56) 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 CPconserving 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/3E 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 eq. (2.57) 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  [20,26]). Nine integrals in the collision term (2.36) can be reduced analytically to two integrals as in appendix B. We deal with both diagonal and off-diagonal collision terms in eqs. (2.37), (2.39) and (2.40) 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 eqs. (2.42) and (2.43) for the self-interactions of neutrinos, νν ↔ νν and νν ↔ νν, since the annihilations of electrons and positrons are important for the heating process of neutrinos while the selfinteractions of neutrinos less contribute to this heating process. We treat this collision term from neutrino self-interaction in eq. (A.13) 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 2 e + δm 2 e (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.

Computational method, initial conditions and values of neutrino masses and mixing
We solve kinetic equations for neutrinos of eqs. (2.57) and (2.58) with the following comoving variables instead of the cosmic time t, the momentum p, and the photon temperature T γ , 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 (2.57) 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 ρ in y i (x) are regarded as 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], 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 CPsymmetry, δ 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 2 31 > 0, not in the inverted ordering, ∆m 2 31 < 0, because the results do not change significantly.

Effective number of neutrino species N eff
To describe the process of neutrino decoupling, we first numerically solve a set of eqs. (2.57) and (2.58) 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 .

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 . 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.

Value of the effective number of neutrino species N eff
The effective number of neutrinos N eff can be rewritten, where δρ να = ρ να − ρ eq ν and ρ eq ν = d 3 p (2π) 3 pf eq . In tables 2 and 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 pf 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].   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.

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 eqs. (2.57) and (2.58) 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 2 F ) 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 2 F QED interactions also modify the weak interaction rates in the collision term C[ρ p (t)] and the Hamiltonian for the forward scattering of neutrinos (2.59) at order e 2 G 2 F 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 2 F ).
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 eq. (2.56). One can incorporate corrections (i) in the weak interaction rates by shifting m 2 e → m 2 e + δm 2 e(2) (p, T ), but it is numerically difficult to take into account the momentum-dependent part of δm 2 e(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. But 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 2 F ) 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.

Forward scattering of neutrinos via their self-interactions
In the Hamiltonian (2.59) in the Boltzmann equations (2.57), the forward scattering terms of neutrinos via their self-interactions correspond to Even in the case without lepton asymmetry, . 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, 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 .

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.

Relic neutrino spectra in the mass basis
We present relic neutrino spectra in the mass basis by solving a set of eqs. (2.57) and (2.58) in the mass basis directly. We can also obtain the same result by transforming relic neutrino spectra in the flavor basis through eq. (2.7).
In the mass basis, the neutral and charged currents including left-handed neutrino fields in eq. (2.23) are given by, using ν α = i=1,2,3 U αi ν i as in eq. (2.5), 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 H p → U † PMNS H p U PMNS analogous to M 2 diag = U † PMNS M 2 U PMNS in eq. (2.57) 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 * ei , as in eq. (4.2). 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 eqs. (2.57) and (2.58) 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.

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 T ν (t) is the effective neutrino temperature, which is ∝ a(t) −1 and normalized asT ν → T γ in high temperature limit. Under this definition ofT ν (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 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. (4.6) Neutrino number density and energy density per one degree of freedom in the current universe are also parametrized as whereñ 0 andρ 0 are given bỹ 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 ofT ν (t 0 ) ∆m 2 21 8.6 meV, |∆m 2 31 | 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.15T ν (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, 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 d ν i ≡ 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.  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 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.

Helicity of relic neutrinos -Majorana 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. 13 1.01 0.91 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. 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 nonrelativistic 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], where f ν i (p ν , t) is given by eq. (4.4) 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 where ν s i 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 righthanded 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 whereν s i denotes a sterile state of anti-neutrino. From eqs. (4.12) and (4.13), 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], (4.14) Then number density and energy density summed over helicity for Majorana and Dirac neutrinos also differ by a factor of two, 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.

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 where δn c ν i 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 c ν i , effects of e ± -annihilations into ν,ν during neutrino decoupling have not been taken into account simultaneously. For m ν i < 0.15eV, 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.15eV).
In the following, we discuss only the case where δn c ν i < 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, where δn d ν i is the enhancement factor by e ± -annihilations into ν andν during neutrino decoupling given in table 7.
0.53 50 12 100 50 200 300 Table 8: The enhancement factor, δn c ν i , due to neutrino clustering by our Galaxy and nearby galaxies for given values of neutrino masses [36].

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.

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, The total capture rate for the CνB in this process, Γ CνB , can be written 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 (5.5) where N T = M T /M3 H is the number of tritium, M T is the total tritium mass in the experimental setup, and M3 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 loc ν i (p ν , s ν ) is the local momentum distribution for relic cosmic neutrinos around the Earth, which satisfies n loc ν i (s ν ) = dp 3 ν (2π) 3 f loc ν i (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) where V ud 0.9740 is a component of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, m3 H 2808.921 MeV and m3 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παZE e /|p e | 1 − e −2παZEe/|pe| , 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) where K 0 end is the beta decay endpoint kinetic energy for massless neutrinos given by E ν i is so small compared to K 0 end and m e that we can safely neglect E ν i in eq. (5.8). Then we obtain Γ i with 1% precision substituting eq. (5.6) into eq. (5.5), where v ν i is the (unnormalized) average magnitude of velocity for ν i given by 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 eq. (5.10) with 1% precision. Here p 0 ∼ 3.15T ν (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 eq. (5.10) 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 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 helicitydependent part in Γ i is already suppressed by v ν i .

Majorana vs Dirac neutrinos
For non-relativistic neutrinos, i.e., v i 1, if we set v ν i = 0 in eq. (5.10), Γ 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 Γ M i is twice that for Dirac neutrinos [34], 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], Note again that the approximations in eqs. (5.12) and (5.13) 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 eq. (5.10) should be included precisely. In both neutrino mass ordering we take the following values of the PMNS matrix, |U e1 | 2 0.681, |U e2 | 2 0.297, |U e3 | 2 0.0222. (5.15) Note that neutrino squared-mass differences and neutrino mixing parameters currently include a few percent (about 10%) uncertainties even at 1σ (3σ) confidence level.

Values
In table 9, we show values of the capture rates on 100 grams of tritium in both the cases of NO and IO for Majorana and Dirac neutrinos with m lightest = 0. δΓ d i denotes the differences between the cases with and without effects of e ± -annihilation during neutrino decoupling and δΓ c i 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.  Table 9: Capture rates of relic cosmic neutrinos on 100 grams of tritium in unit of year −1 with m lightest = 0. δΓ d i is the differences between the cases with and without effects of e ±annihilation during neutrino decoupling and δΓ c i is the differences with and without gravitational clustering for relic neutrinos in nearby galaxies.

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, δΓ c,d i , and the estimated error of the capture rates for relic neutrinos on tritium in more detail.
To observe δΓ c,d i , we need a large number of events to satisfy typically 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 δΓ c,d i ∼ (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 ei | 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, δΓ c,d i for m lightest = 0 is correctly estimated since uncertainties of |U ei | are canceled out in δΓ c,d i . 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 GT 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. But, 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.

β-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, The β-decay spectrum and the capture rate for the β-decay process are given by [70] (see also appendix D) where 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), Then the maximal energy for the emitted electron in the β-decay process called the energy at β-decay endpoint is where m lightest is the lightest neutrino mass. We can see that the β-decay spectrum dΓ β /dE e vanishes for E e = E end e . Then the total tritium β-decay rate is obtained as 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 max e . 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) where E CνB,i e is the emitted electron energy from the CνB signal, ν i + 3 H → e − + 3 He given by eq. (5.8).
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 wouldbe 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, Substituting eq. (5.4) into eq. (5.24), the smeared spectrum of the emitted electron from the CνB signal can be written as eq. (5.26) is a Fredholm integral equation of the first kind and dΓ i dEe is a would-be observed quantity. After solving eq. (5.26) 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 e1 |. 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 e3 | while in the normal ordering that is composed of ν 1 with large |U e1 |.

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 the 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.

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 2 F ) 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 tables 4 and 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.

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 2 F ) in the homogeneous and isotropic universe, following refs. [9,39]. We consider the general form of the collision term in this case, 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: and decompose momentum integrations into the radial and angle components, For four-Fermi interaction processes at order of O(G 2 F ), all of |M| 2 have two kinds of forms, 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, In both cases of eqs. (B.6) and (B.7), we can perform all integrals for angle components in eq. (B.5) so that D(q 1 , q 2 , q 3 , q 4 ) in the case of eq. (B.6) reduces to while in the case of eq. (B.7), D(q 1 , q 2 , q 3 , q 4 ) is given by where D 1,2,3 are defined in eq. (A.12).
(B.13) D 3 is equal to that in eq. (B.11) 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.
C Kinematics for ν i + 3 H → e − + 3 He and 3 H → e − + 3 He+ν 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, m3 H and m3 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 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 ν +E3 He ), which yields From these conditions, the maximal energy of the electron for 3 H → e − + 3 He +ν i is given by The endpoint energy of the electron for the tritium β-decay is also given by 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 where we neglect the terms proportional to |p ν | 2 and |p ν ||p e | and leave the term proportional to E ν i m3 H because of m3 H |p e | |p ν |. For m3 H m e , the difference between E CνB,i e and E end is (C.8) Since E CνB,i e − 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. D Cross section for ν i + 3 H → e − + 3 He and decay rate for 3 H → e − + 3 He +ν 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Γ β /dE e .
D.1 Cross section for ν i + 3 H → e − + 3 He 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.1) where s = (p ν i + p3 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 = (m3 H + E ν i ) 2 − |p ν | 2 = m 2 3 H + 2m3 H E ν i + m 2 ν i , t = (E e − E ν i ) 2 − |p e − p ν | 2 (m e − m ν i ) 2 + 2|p e ||p ν | cos θ. (D.2) Using also dt/d cos θ = 2|p e ||p ν |, we obtain The matrix element for ν i + 3 H → e − + 3 He is effectively given by 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 GT √ 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 Using the completeness relations, we obtain the relation of Dirac spinors for 3 H, 3 He, and e − , s j =±1/2 u jūj = ( / p j + m j ), (D. 8) and for neutrinos with their helicity s ν , where S ν i is the spin vector for neutrinos given by (D.10) In the massless limit, the previous relation of the Dirac spinor for neutrinos becomes where we used mS µ = p µ and p µ S µ = 0. Using the above relations, we rewrite eq. (D.7) as T αβ 1 = 1 2 tr γ α 1 − γ 5 / p ν i + m ν i 1 + 2s ν γ 5 / S ν i γ β 1 − γ 5 / p e + m e , (D.12) T γδ 2 = tr γ γ F − Gγ 5 / p n + m n γ δ F − Gγ 5 / p p + m p . (D.13) Then we obtain T αβ 1 T 2αβ as T αβ 1 T 2αβ = 32 (G + F ) 2 [(p e · p3 He ) (p ν i · p3 H )] + (G − F ) 2 [(p e · p3 H ) (p ν i · p3 He )] + G 2 − F 2 m3 H m3 He (p e · p ν i ) − 64s ν m ν i (G + F ) 2 [(p e · p3 He ) (S ν i · p3 H )] + (G − F ) 2 [(p e · p3 H ) (S ν i · p3 He )] + G 2 − F 2 m3 H m3 He (p e · S ν i ) . (D.14) In the rest frame of 3 H, neglecting the momentum of 3-helium |p3 He |/m3 He ∼ (m3 H − m3 He )/m3 He ∼ O(10 −4 ), T αβ 1 T 2αβ is given by 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 (D. 16) D.2 Decay rate for 3 H → e − + 3 He +ν i The decay rate of the β-decay follows the standard formula at the rest frame of tritium, Γ β = 1 2 9 π 5 m3 H d 3 p e d 3 p ν i d 3 p3 He E e E ν i E3 He |M| 2 δ 4 (p3 H − p e − p ν i − p3 He ), = 1 2 6 π 4 m3 H dE e dE ν i |M β | 2 , (D. 17) 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, After similar calculations in appendix D.1, |M β | 2 for β-decays at rest of tritium is written as where we neglect the momentum of 3 He due to p3 He m3 He . In addition, we neglect the second term in eq. (D.23) since |p e | ∼ m3 H − m3 He E e . Thus, |M β | 2 approximately becomes (D.24) Plugging eq. (D.24) into eq. (D.22), we obtain 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 ), } is the endpoint energy of the tritium β-decay given by eq. (C.4) in appendix C.