MeV-scale reheating temperature and thermalization of oscillating neutrinos
by radiative and hadronic decays of massive particles
Takuya Hasegawa,1, 2 Nagisa Hiroshima,1, 3 Kazunori Kohri,1, 2, 4
Rasmus S. L. Hansen,5, 6 Thomas Tram,6, 7 and Steen Hannestad6
1
2

Theory Center, IPNS, KEK, Tsukuba 305-0801, JapanThe Graduate University of Advanced Studies (Sokendai), Tsukuba 305-0801, Japan
3

RIKEN Interdisciplinary Theoretical and MathematicalarXiv:1908.10189v1 [hep-ph] 21 Aug 2019

Sciences (iTHEMS), Wako, Saitama 351-0198, Japan
4

Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
5

Max-Planck-Institut f ür Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany
6

Department of Physics and Astronomy, University of Aarhus,
Ny Munkegade 120, DK8000 Aarhus C, Denmark
7

Aarhus Institute of Advanced Studies (AIAS),

Aarhus University, DK8000 Aarhus C, Denmark(Dated: August 28, 2019)
From a theoretical point of view, there is a strong motivation to consider an MeV-scalereheating temperature induced by long-lived massive particles with masses around the weak
scale, decaying only through gravitational interaction. In this study, we investigate lower
limits on the reheating temperature imposed by big-bang nucleosynthesis assuming bothradiative and hadronic decays of such massive particles. For the first time, effects of neutrino self-interactions and oscillations are taken into account in the neutrino thermalization calculations. By requiring consistency between theoretical and observational values oflight element abundances, we find that the reheating temperature should conservatively be
TRH & 1.8 MeV in the case of the 100% radiative decay, and TRH & 4 − 5 MeV in the caseof the 100% hadronic decays for particle masses in the range of 10 GeV to 100 TeV.

2I.

INTRODUCTION

In the standard big-bang cosmology it is normally assumed that radiation components (photons,
electrons/positrons, and neutrinos) were perfectly thermalized, and energy of radiation dominatedthe total energy density of the Universe well before the beginning of big-bang nucleosynthesis(BBN). In a modern picture of the early Universe, this radiation-dominated epoch is expected tobe realized after a decay of a massive particle such as the inflaton, the particle associated with the
inflaton field driving inflation, or another massive particle such as the curvaton. If such massive
particles abundantly existed in the early Universe, their non-relativistic energy could dominate
the total energy, and then an early matter-dominated epoch should have been realized before theradiation-dominated epoch. Therefore, particle production caused by their decays and subsequent
entropy production (called reheating) dramatically modify the thermal history of the Universe.
The Universe could experience the reheating more than once after inflation depending on thefundamental theory of particle physics. Since many theoretical models have been proposed as atheory beyond the standard model of particle physics, it is required to have some ways to find thetrue theory of nature. The one of the approaches is to investigate the possible value of “reheating
temperature” which is defined by the cosmic temperature when the radiation-dominated epoch
just started. This is because the reheating temperature is related to the property of the massive
particles, and we can constrain the theories through the observational bound on the reheating
temperature.
As for a candidate of inflaton field or curvaton field, a lot of unstable massive scalar fields, e.g.,
moduli, dilaton fields, are predicted in particle physics theories beyond the standard model such as
supergravity or superstring theory. They tend to dominate the total energy of the Universe during
their oscillation epochs. It is notable that they typically have masses at or above the weak scale anddecay only through gravitational interaction. This means that they have long lifetimes of O(0.1) sec– O(10) sec, and the reheating temperature after their decay is expected to be O(1) MeV. Since
neutrinos decoupled from the thermal plasma at around the cosmic temperature T ∼ O(1) MeV,
they would have suffered imperfect thermalization due to the late-time entropy production causedby their decay. Thus, we have a strong motivation to observationally test this kind of cosmologicalscenarios with decaying particles which induce the MeV-scale reheating temperature.
The theory of BBN based on the standard big-bang cosmology, i.e., standard BBN, successfully
explains observational light element abundances (see e.g., Ref. [1] and references therein), and thetheory say that the light elements are synthesized at around T ∼ O(0.01) MeV – O(0.1) MeV.

3
As we shall see in Sec. IV, BBN is highly sensitive to the neutrino abundances. Therefore, we can
examine the MeV-scale reheating scenarios by using BBN as a probe.
In Ref. [2], lower bounds on reheating temperature have been studied in terms of BBN for
the first time. They have shown that an incomplete thermalization of neutrinos gives the mostsignificant effect on BBN assuming that the 100% of the long-lived massive particles decay into
electromagnetic radiations such as photons or charged leptons. Because of a competition between
decreases and increases of the produced amount of 4 He by the imperfect thermalization of the neutrinos, we can constrain the reheating temperature. As a result, they have obtained a conservativelower bound on the reheating temperature TRH > 0.5 MeV – 0.7 MeV (95% C.L.).
Afterwards, in Ref. [3], they discussed hadronic decays of massive particles, i.e., direct decays
into quarks and/or gluons which immediately fragment into hadrons such as pions, kaons, or nucleons. The thermalization of radiations proceed in exactly the same way as in the case where
the 100% of the massive particles decay into electromagnetic radiations. This is because almost
all of the kinetic energy of hadrons are transferred into radiation through Coulomb scatteringwith background electrons/positrons or inverse-Compton like scattering with background photons,
and a neutral pion π 0 immediately decays into two photons. In the case of the hadronic decay,interconverting reactions between ambient protons and neutrons induced by emitted hadrons are
extraordinarily important because they increase the neutron to proton ratio which is a key parameter of resultant abundances of light elements. As a result, they obtained a lower bound TRH > 2.5– 4 MeV (95% C.L.) depending on the mass of the long-lived massive particles and their branchingratio into hadrons.
Subsequently, two- and three-flavor neutrino oscillations were respectively considered in thethermalization process of neutrinos in Refs. [4] and [5] where they obtained a lower bound TRH >2 MeV (95% C.L.) and TRH > 4.1 MeV (95% C.L.) assuming radiative decay of the massive
particles.
Some other cosmological probes other than BBN are also sensitive to neutrino abundances.
Here we briefly refer to the recent papers which focused on this topic. In Ref. [3], they discussed
possible effects of an incomplete thermalization of neutrinos on a temperature anisotropy andpolarization of Cosmic Microwave Background (CMB) and a galaxy power spectrum of Large Scale
Structure (LSS). Ref. [6] obtained a combined constraint TRH > 4 MeV (95% C.L.) by considering
BBN, CMB (WMAP) and LSS (2dF Galaxy Redshift Survey).
After that, authors in Ref. [7] have updated the CMB and LSS data by using WMAP threeyear data and SDSS luminous red galaxies data, and they obtained TRH > 2 MeV (95% C.L.).

4
Similar analysis have been done in Ref. [8] by using WMAP five-year data and SDSS luminous red
galaxies data where they obtained TRH > 2 MeV (95% C.L.) from CMB only and TRH > 3.2 MeV(95% C.L.) from CMB by including external priors from SDSS red luminous galaxy survey andthe constraint from the comic age. In addition, authors in Ref. [5] have reported a new constraint
TRH > 4.3 MeV (95% C.L.) from CMB by using Planck 2015 data.
In this paper, we extend the study of Ref. [3] by considering neutrino oscillation and neutrino
self-interaction in the calculation of the neutrino thermalization. We assume both radiative andhadronic decays of the massive particles and give an updated bound on the reheating temperature
set by BBN. This is the first study that consider effects of neutrino self-interactions on the neutrino
thermalization to constrain the reheating temperature.
The structure of this paper is as follows. In Sec. II, we introduce the formalism of the neutrino
thermalization assuming the MeV-scale reheating temperature. In Sec. III, we show the results ofneutrino thermalization in the reheating process, and describe how neutrino oscillation and neutrino
self-interaction affect the thermalization process. The results of BBN are shown in Sec. IV where
we discuss effects of both radiative and hadronic decays on light element abundances. Finally, wedraw our conclusion in Sec. V.

5II.

REHEATING PROCESS AND NEUTRINO THERMALIZATIONIn this section, we describe the neutrino thermalization in the low-reheating-temperature Universe and introduce the key equations.
As described in the previous section, there are some candidate particles in theories going beyond the Standard Model of particle physics which are weakly interacting and decay at around
BBN. Here, we call the long-lived massive particles just “massive particles” and label them X.
We assume the energy density of the massive particles dominates those of other particles at an
initial time and the Universe is completely matter-dominated before the massive particles start todecay. In this setting, the entropy production caused by the decay of X induces particle production
via electromagnetic, weak, and strong interactions. As a result, photons and charged leptons are
rapidly thermalized in the reheating via electromagnetic interactions, whereas the weakly interacting neutrinos are slowly created in the thermal bath of photons and charged leptons. 1 Since
neutrinos decouple from the thermal bath at around a temperature of O(1) MeV, neutrinos should
not be fully thermalized if TRH ∼ O(1) MeV. A degree of thermalization of neutrinos affects thelight element abundances [2]. For this reason, it is especially important to accurately calculate thethermalization of neutrinos in the reheating for TRH ∼ O(1) MeV. Therefore, next we look into
the dynamics of neutrino thermalization in the thermal plasma.
In the Universe with a temperature of O(1) MeV, electrons and positrons are the only charged
leptons which are abundant in the system since the abundances of muons and tau leptons arestrongly suppressed by Boltzmann factors. Therefore, neutrinos are mainly produced in the annihilation process of electrons and positrons, e− + e+ → να + ν̄α where α = e, µ, τ . Since electron
neutrino (νe ) is not only produced by the neutral-current weak interaction but also by the chargedcurrent one, it tends to be produced more than µ neutrinos (νµ ) and τ neutrinos (ντ ) when all
neutrinos are not fully thermalized. Consequently, neutrino oscillations play a role in equilibrating
neutrino abundances in this case, and we have to simultaneously consider the neutrino production
by collisions and flavor oscillations.
Our treatment of neutrino oscillation is consistent with that of Ref. [4]. That is, we adopt
the effective two-flavor mixing scheme which is a good description to approximately include full
three-flavor mixings when the collision rates of νµ and ντ are identical, and one mixing angle is1

In the case where the massive particles decay intro quarks and/or gluons, they fragment into mesons and baryons
after the hadronization, and almost all the kinetic energy of hadrons are transferred into radiation due to thereason mentioned before. Therefore, contributions of secondary neutrinos produced by the decay of hadrons canbe neglected. However, there is also another possibility that the massive particles directly decays into neutrinos,
e.g., X → να + ν̄α where α = e, µ, τ [6]. In this case, electromagnetic plasma is produced only from neutrinos via
weak interaction, which gives totally different results of the neutrino thermalization and BBN. We do not consider
the possibility in this paper.

6
predominantly important compared to others (see e.g., Ref. [9] for more details on the effective
two-flavor mixing scheme). The former condition is well satisfied because of the absence of muonsand tau leptons in the system with a temperature of O(1) MeV. On the other hand, the lattercondition is only approximately satisfied since the reactor neutrino mixing θ13 is known to be nonnegligible compared to other mixings, namely the solar neutrino mixing θ12 and the atmospheric
neutrino mixing θ23 [10]. In later sections, however, we show that the effect of θ13 on BBN isvery small compared to that of θ12 or neutrino self-interaction irrespective of the mass orderingof neutrinos. Therefore, the effective two-flavor mixing scheme (with the solar neutrino mixing)
gives a good description of the full three-flavor mixings at least for the current purpose. In thisscheme, we label the degenerate state of νµ and ντ as x neutrinos (νx ) and consider the flavor
mixing between νe and νx . Also, we label the other neutrino species which does not mix with otherflavors as spectator neutrino (νsp ) in this two-flavor treatment.
In general, neutrino states can be described by a one-particle irreducible density matrix ̺p ≡̺(p, t) where p ≡ |p| is the absolute momentum. 2 Since we focus on the effective two-flavor mixing,
the density matrix is expressed in terms of a 2 × 2 Hermitian matrix, and we label each element of
the density matrix as̺p = 

ρee ρexρ∗ex ρxx



.

(2.1)

The diagonal part of the matrix corresponds to the distribution function of mixed neutrinos (i.e., νeand νx ), that is, ρee = fνe and ρxx = fνx , while the off-diagonal elements represent the quantum
coherence among neutrinos with different flavors. In the current study, the chemical potentials
of neutrinos are set to be zero. Under this assumption, the density matrix of neutrinos is equal
to that of anti-neutrinos, i.e., ̺p = ̺¯p , and they have the same abundance. Therefore, it is notnecessary to follow the time evolution of anti-neutrinos separately from that of the corresponding
neutrinos.
We can obtain the time evolution of the neutrino density matrix by solving the momentumdependent Quantum Kinetic Equation [11, 12] which is formally written as
∂̺p
∂̺p
d̺p
=
−Hp
= −i [Hp , ̺p ] + C(̺p ) .
dt
∂t
∂p

(2.2)

In the above equation, the term including the Hubble parameter H corresponds to the effect of the2

Since we focus on the Universe with a temperature of O(1) MeV, we can neglect tiny neutrino masses which are
known to be sub-eV scale [1]. In this case, energy of neutrinos E is equal to its absolute momentum, i.e., E = p.

7expansion of the Universe, and C(̺p ) is the collision term of neutrinos expressed asC(̺p ) = 

Rν e

−Dρex

−Dρ∗ex

Rν x



,

(2.3)

where Rνe and Rνx are the production rates of νe and νx , respectively. Also, D is the collisionaldamping rate which breaks the flavor coherence among different flavors of neutrinos. In thispaper, we adopt a simplified treatment of the damping effects discussed in Ref. [5] and neglect theadditional contributions such as “damping-like terms” which appear in Ref. [11]. In the current
study, we consider the collisional processes a(k) + b(p) → c(k′ ) + d(p′ ) shown in Table. I including
those of neutrino self-interaction. In this case, the expressions of the repopulation and the damping
terms are [13],Rνα (k) = 2πZdΠk′ dΠp′ dΠp δE (kp|k′ p′ )
×

X
i


V 2 [να (k), ν̄α (p)|i(k′ ), ī(p′ )] fi (Ek′ )fī (Ep′ )(1 − fνα (k))(1 − fν̄α (p))

− fνα (k)fν̄α (p)(1 − fi (Ek′ ))(1 − fī (Ep′ ))X+V 2 [να (k), j(p)|να (k′ ), j(p′ )] fνα (k′ )fj (Ep′ )(1 − fνα (k))(1 − fj (Ep ))j

D(k) = π

XZ
α

×


−fνα (k)fj (Ep )(1 − fνα (k′ ))(1 − fj (Ep′ )) ,(2.4)

dΠk′ dΠp′ dΠp δE (kp|k′ p′ )XiV 2 [να (k), ν̄α (p)|i(k′ ), ī(p′ )] fi (Ek′ )fī (Ep′ )(1 − fν̄α (p))

+fν̄α (p)(1 − fī (Ep′ ))(1 − fi (Ek′ ))X+V 2 [να (k), j(p)|να (k′ ), j(p′ )] fνα (k′ )fj (Ep′ )(1 − fj (Ep ))

(2.5)

j


+fj (Ep )(1 − fj (Ep′ ))(1 − fνα (k′ )) ,where α = e, x, and k, p, k′ and p′ are absolute momenta of the particle a, b, c and d, respectively.
3

d p
′ ′
(1) (E + E − E ′ − E ′ ) is the 1D Dirac delta function
Also, dΠp ≡ (2π)
pp
3 , and δE (kp|k p ) ≡ δkkcorresponding to energy conservation for each process. The summation index i runs over electrons
and all flavors of neutrinos other than να (i.e., νβ where β 6= α), while j runs in addition overpositrons, να , and all flavors of anti-neutrinos. The expression of V 2 is written as

V 2 [a(p), b(k)|c(p′ ), d(k′ )] = (2π)3 δ(3) (k + p, k′ + p′ )Na2 Nb2 Nc2 Nd2 S|M |2 (a(p), b(k)|c(p′ ), d(k′ )) , (2.6)

8Process (α 6= β)

S |M |2
25 G2F [ (2 sin2 θW ± 1)2 (k · p′ )(p · k ′ ) + 4 sin4 θW (k · k ′ )(p · p′ )+2 m2e sin2 θW (2 sin2 θW ± 1)(k · p) ]25 G2F [ (2 sin2 θW ± 1)2 (k · k ′ )(p · p′ ) + 4 sin4 θW (k · p)(k ′ · p′ )−2 m2e sin2 θW (2 sin2 θW ± 1)(k · p′ ) ]2
5 2
2
2 GF [ (2 sin θW ± 1) (k · p)(k ′ · p′ ) + 4 sin4 θW (k · k ′ )(p · p′ )−2 m2e sin2 θW (2 sin2 θW ± 1)(k · p′ ) ]I.

e− + e+ → να + ν̄αII.

να + e + → να + e +III.

να + e − → να + e −IV.

να + να → να + ναV.

να + νβ → να + νβ

25 G2F (k · p)(k ′ · p′ )VI.

να + ν α → να + ν α

27 G2F (k · p′ )(p · k ′ )VII.

να + ν β → να + ν β

25 G2F (k · p′ )(p · k ′ )VIII.

να + ν α → νβ + ν β

25 G2F (k · p′ )(p · k ′ )26 G2F (k · p)(k ′ · p′ )Table. I: Collision process a(p) + b(k) → c(p′ ) + d(k′ ) which contributes to the thermalization ofneutrinos of each flavor να , νβ (α, β = e, µ, τ where α 6= β). The process I is the production
process of neutrinos due to electron annihilation, the processes II - III are the scattering processes
between neutrinos and electrons, and the processes IV - VIII are the self-interaction processes
among neutrinos. Here, θW is the Weinberg angle, GF is the Fermi-coupling constant, S is thesymmetry factor, and |M |2 is the squared scattering matrix element. The positive sign in theexpression is for νe and the minus sign for νµ or ντ (i.e., for νx and νsp ). The expressions of theprocess I, II, IV and V are also applied to the corresponding anti-neutrinos.
where S|M |2 (a(p), b(k)|c(p′ ), d(k′ )) is the squared scattering matrix element for the processes inTable. I summed over initial and final spins, and symmetrized over identical particles in the initialpand the final state. Also, Ni ≡ 1/2Ei where Ei is the energy of particle i (i = a, b, c, d) andδ(3) (k + p, k ′ + p′ ) ≡ δ(3) (k + p − k′ − p′ ) is the 3D Dirac delta function corresponding to themomentum conservation. As for the processes in Table. I, we analytically reduce the dimensionof momentum integrals in the above expressions from nine to two and calculate the full collision
terms without any simplifying assumptions in the same way as in Ref. [13].
In the expression of Eq. (2.2), Hp is neutrino Hamiltonian which is expressed as
√


Z
√
Eν
M2 8 2 GF p El
+ 2 + 2 GF dΠp′ (̺p′ − ̺¯∗p′ ) ,−Hp = Hp, vac + Hp, mat =
2p
3
m2W
mZ

(2.7)

where GF is the Fermi-coupling constant. Also, mW and mZ are the masses of W and Z bosons,respectively. In the above expression, the first term, Hp, vac , is the contribution which induces thevacuum oscillation where M is the mass matrix in flavor basis. The mass matrix M is related to theone in mass basis M as M2 = U M2 U † where U is the PMNS matrix. In the effective two-flavor

9
mixing scheme,
M2 =



m2
 1
0

0
m22





 , U =

cos θ

sin θ

− sin θ cos θ



,

(2.8)

where δm2 ≡ m22 − m21 is the squared-mass difference and θ is the mixing angle in vacuum between
νe and νx . Also, the second term in the Hamiltonian, Hp, mat , corresponds to the matter potentials
which arise from coherent scatterings between neutrinos and charged-leptons. In the term, El
corresponds to the total energy density of charged leptons, while Eν to that of neutrinos:




Z
ρe 0
ρ
ρ
ν
ν
ex
 , Eν = dΠp′ p′ (̺p′ + ̺¯∗p′ ) =  e,El = 
∗0 0ρνex ρνx(2.9)

where ρe = ρe− + ρe+ is the total energy density of electrons and positrons, while ρνe and ρνxR
are those of νe and νx , respectively. Also, we have defined ρνex ≡ dΠp′ p′ (ρex + ρ̄∗ex ) and ρ∗νex ≡RdΠp′ p′ (ρ∗ex + ρ̄ex ). In the expression of El , we have neglected the existence of muons or tau

leptons due to their large masses. 3 The asymmetric part of Hp, mat is often assumed to vanish
when ρ = ρ̄ for neutrinos and the number density of electrons and positrons are identical. For the
diagonal part of the Hamiltonian this is true, however the off-diagonal part gets a contribution
from the neutrinos as shown in Eq. (2.7) since ̺∗p 6= ̺p .

As for the oscillation parameters, we use the best fit values of the mass-squared differences andmixing angles reported in Ref. [10]: 4
δm212 = 7.55 × 10−5 eV2 , sin2 θ12 = 3.20 × 10−1 ,(2.10)

δm213 = 2.50 × 10−3 eV2 , sin2 θ13 = 2.160 × 10−2 (NO) ,(2.11)

δm213 = −2.42 × 10−3 eV2 , sin2 θ13 = 2.220 × 10−2 (IO) ,(2.12)

where “NO” (“IO”) means normal (inverted) mass ordering of neutrinos, respectively. For thenumerical calculation, we rewrite the 2 × 2 density matrix with polarization vectors (P0 , P ):ρee ρex = 1 [ P0 (p) σ0 + P (p) · σ ] ,(2.13)
̺p = 
2
ρ∗ ρ
ex

xx

where σ0 = 1 is the identity matrix and σ = (σx , σy , σz ) are the Pauli matrices. With P =(Px , Py , Pz ), distribution functions of mixed neutrinos can be written asf νe =
3

4

1
1(P0 + Pz ) , fνx = (P0 − Pz ) .
22(2.14)In the effective two-flavor mixing scheme, we need to treat both νµ and ντ in the same way. Therefore, we do notconsider background muons or tau leptons whose contribution is very small compared to that of electrons.
The atmospheric neutrino mixing (θ23 , δm223 ) is irrelevant to the oscillation between νe and νx in the effective
two-flavor mixing scheme. Therefore, we do not use the value in this paper.

10
In addition, we can rewrite the expression of Eq. (2.2) as follows:→−Ṗ = H × P − D (Px x + Py y) + (Rνe − Rνx ) z ,(2.15)

Ṗ0 = Rνe + Rνx ,(2.16)

which leads toṖνe = Hx Py − Hy Px + 2 Rνe ,(2.17)Ṗνx = Hy Px − Hx Py + 2 Rνx ,1
P˙x = Hy (Pνe − Pνx ) − Hz Py − D Px ,2
1P˙y = Hz Px − Hx (Pνe − Pνx ) − D Py ,2

(2.18)(2.19)(2.20)

where x, y and z are unit vectors, and we defined Pνe ≡ P0 + Pz and Pνx ≡ P0 − Pz . Each−→component of the neutrino potential H, i .e., Hi = Tr (Hp σi ) where i = x, y, z, is written as√Z
16 2 GF p
δm2
sin 2θ −Hx =dΠp′ p′ Px ,
2p
3 m2Z
Z
√
Hy = 2 2 GF dΠp′ Py ,Hz = −

δm2
cos 2θ + Hmat .
2p

The second term in Hz is the matter contribution which is explicitly written as√ρνe − ρνx
8 2
ρe
+
,
Hmat = −
GF p
3
m2W
m2Z
√

Z ∞
Ee
4 2
ge
dp′ p′2= − 2 GF p
2
3π
exp(Ee /Tγ ) + 1
mW 0

Z ∞
gν
′ ′3
dp p (fνe − fνx ) ,+ 2
mZ 0
where Tγ is the photon temperature, and Ee =p

(2.21)(2.22)(2.23)(2.24)(2.25)

p2 + m2e is the energy of electrons. Also, ge = 4

is statistical degree of freedom of electrons and gν = 2 is that of neutrinos of each flavor.
With the matter potential Hmat , the mass-squared difference and the mixing angle in vacuum
are modified in medium by MSW effect as follows [14, 15]:
s
2


δm2M
δm2 2 2δm2,=
sin 2θ + −cos 2θ + Hmat2p2p2p
 2 2δm
sin2 2θ
2p2
sin 2θM = 
2

2 ,
δm2
δm2
2
sin
2θ
+
−
cos
2θ
+
H
mat
2p
2p

(2.26)(2.27)

11
where δm2M and θM are the in-medium mass-squared difference and the mixing angle, respectively.
We note here that in Eqs. (2.26)-(2.27) we simplify the expressions by neglecting the small contributions of the off-diagonal part of the matter potential Hp, mat in Eq. (2.7) (see Ref. [9] for the exact
expressions of the MSW effect). As can be seen from Eqs. (2.26) and (2.27), the matter potential
2−1|Hmat | ∝ Tγ5 dominates the vacuum one |Hvac | ≡ | − δm
2p cos 2θ| ∝ Tγ , i.e., |Hmat | >> |Hvac |,at high temperature such as Tγ > O(10) MeV and θM ∼ 0 holds for most energy modes. 5 On
the other hand, the opposite hierarchy, i.e., |Hvac | >> |Hmat |, holds at low temperature and themixing parameters take the same values as those in vacuum: θM ∼ θ, δm2M ∼ δm2 . For neutrinos
with momentum p = hpi ∼ 3.15 Tγ where h·i means a thermal average, the level crossing between
these potentials, i.e., |Hmat | ∼ |Hvac |, occurs at the temperature Tc


1/6

δm212
 3 MeV−1/3
2.5×10−3 eV2
Tc ∼ GF (δm2 cos 2θ)1/6 ∼1/6 ,
δm213

 5 MeV2−5
7.5×10eV(2.28)

where, in the above evaluation, we have replaced p with hpi ∼ 3.15 Tγ and approximated(cos 2θ)1/6 ∼ 1 which is well satisfied for θ12 and θ13 . Therefore, neutrino oscillation becomes
effective at around a temperature of O(1) MeV for the solar neutrino mixing (θ12 , δm12 ) andthe reactor neutrino mixing (θ13 , δm13 ), which is the reason for taking its effect on the neutrino
thermalization into account.
On the other hand, since νsp decouple from flavor mixings of neutrinos, the time evolution ofthis neutrino species is just given by the classical Boltzmann equation:
dfνsp
∂fνsp
∂fνsp
=
−Hp
= C(fνsp ) ,
dt
∂t
∂p

(2.29)

where fνsp is the distribution function of νsp , and C(fνsp ) is the collision term whose expression isequal to that of νx , i.e., C(fνsp ) = Rνx (see Eq. (2.4)).
In order to calculate the thermalization process of neutrinos in the expanding Universe, we alsoneed to compute the energy conservation equation:dρ= −3H(ρ + P ) ,dt(2.30)

which can be expressed as the time evolution of the photon temperature Tγ :dρ−ΓX ρX + 4H(ργ + ρν ) + 3H(ρe + Pe ) + dtν
dTγ
=−
,
∂ργ
dt
|a(t) + ∂ρe |a(t)∂Tγ5(2.31)∂Tγ

We assume here that the neutrinos are thermalized with photons and have a temperature Tγ for simplicity.

12
where a(t) is the scale factor at the cosmic time t, ΓX is the decay rate of the massive particles,
whereas ρ and P are the total energy density and the total pressure, respectively:ρ = ργ + ρe + ρν + ρX ,
Z ∞
Ee
ge
π2 4
Tγ + 2dp′ p′2=15
2π 0exp(Ee /Tγ ) + 1
Z ∞
gν
+ 2
dp′ p′3 (fνe + fνx + fνsp ) + ρX ,2π 0
P = Pγ + Pe + Pν ,
Z ∞
1
p′4
ge
π2 4
Tγ + 2
dp′
=
45
6π 0
Ee exp(Ee /Tγ ) + 1
Z ∞
gν
+ 2
dp′ p′3 (fνe + fνx + fνsp ) .
6π 0(2.32)(2.33)

Here, ργ (Pγ ), ρe (Pe ), ρν (Pν ) and ρX mean the energy density (pressure) of photons, electrons,
neutrinos and the massive particles, respectively. The total energy density and the total pressureof neutrinos are a sum of three contributions: ρν = ρνe + ρνx + ρνsp , Pν = Pνe + Pνx + Pνsp . In
addition, the Hubble parameter H is obtained by solving the Friedmann equation:
r
8πGρ
ȧ
.
H≡ =
a
3

(2.34)In the above expression, we can obtain the time evolution of ρX by solving the Boltzmann equationof the massive particles X:dρX= −ΓX ρX − 3HρX ,dt(2.35)

which can be integrated analytically for non-relativistic particles X:
ρX = ρX, 0 e−ΓX t ,(2.36)

where ρX, 0 is the initial energy density of X which is assumed to be much larger than those ofother particles, i.e., ρrad ≡ ργ + ρe + ρν . In addition, ΓX is related to TRH through the Hubble
parameter H as follows:ΓX = 3H(TRH ) ,(2.37)

where TRH is the argument of the Hubble parameter. Since the Universe is in the radiationdominated epoch after most of the massive particles decayed and Tγ ∼ TRH is realized, we canapproximately write the Hubble parameter as
H(TRH ) =r

2g∗ π 2 TRH,
90 mpl(2.38)

13
where mpl ∼ 2.435×1018 GeV is the reduced Planck mass, and g∗ = 10.75 is the relativistic degrees
of freedom in the Universe with a temperature of O(1) MeV. Hence, the relation between TRH andthe decay rate of X is approximately written as
TRH ∼ 0.7



ΓX
sec−1

1/2

MeV .

(2.39)

From the above expression, we can see that TRH ∼ O(1) MeV corresponds to the lifetime of the

6 In order to obtain the neutrino distribution functions andmassive particles τX = Γ−1
X ∼ O(1) sec.

a degree of the neutrino thermalization in the reheating, we simultaneously solve the Eqs. (2.17),(2.18), (2.19), (2.20), (2.29), (2.31), (2.34) and (2.36) from the initial time t = 10−4 sec to thefinal time t = 107 sec corresponding to the cosmic time well before and after BBN, respectively.
We find that the final results are independent of the choice of the initial time as long as theinitial temperature of electromagnetic particles is much higher than TRH . To calculate neutrino
thermalization processes, we use a modified version of LASAGNA code [16, 17] which is, in theoriginal version, a solver of ordinary differential equations for calculating sterile neutrino productionin the early Universe.
In the next section, we show our numerical results of the neutrino thermalization and BBN in
the low-reheating-temperature Universe.

6

Since the actual value of g ∗ depends on the value of TRH , Eq. (2.37) just gives a rough estimate of when theradiation dominant epoch is realized.

14
III.

NUMERICAL RESULTS: NEUTRINO THERMALIZATION IN THE REHEATINGPROCESSIn this section, we show our numerical results of neutrino thermalization. In order to express thetime evolution of the neutrino thermalization, we define the effective number of neutrino species
Neff :
Neff = Neff, νe + Neff, νx + Neff, νsp =X

ρνα /ρνα, std ,(3.1)

α = e, x, sp

where Neff, να is the contribution for each neutrino species, and ρνα, std is the energy density of eachneutrino species in the standard big-bang cosmology. 7 The value of Neff is almost equal to theactual number of neutrino species when all neutrinos are fully thermalized.
Figure 1 shows the relation between TRH and Neff for the cases with and without neutrino
self-interaction. As shown in Fig. 1, the value of Neff increases as TRH becomes large, and thevalue is almost equal to 3.046 above TRH & 10 MeV which is the canonical value in the standard
big-bang cosmology with large TRH [18, 19]. The above threshold value of TRH arises from thefact that weak reaction processes which are responsible for the neutrino thermalization decouple at1−2
5 /(T 2 /m ) ∼ 1, i.e., T
around a temperature Tdec given by Γweak /H ∼ G2F Tdecpldec ∼ (GF mpl ) 3 ∼dec

O(1) MeV where Γweak is the thermal reaction rate of weakly-interacting particles. Therefore, if
TRH is larger than Tdec , neutrinos have enough time to be fully thermalized before decoupling. In
addition, it can be seen from Fig. 1 that both neutrino oscillation and neutrino self-interaction
increase the value of Neff . This is because the production rate of νe is larger than that of νx(Rνe > Rνx ), and thereby neutrino oscillation increases the total production rate of neutrinos
Rν, tot (≡ Rνe + Rνx ). To understand this effect more quantitatively, let us assume that all neutrino
species are almost thermalized. In this case, we can approximate the production rates of νe andνx as [20]Rνe ∼ Ce G2F Tγ5 (feq − fνe ) ,(3.2)

Rνx ∼ Cx G2F Tγ5 (feq − fνx ) ,(3.3)

where feq is the Fermi-Dirac distribution feq = 1/(1 + exp(p/Tγ )), and fνe and fνx are distributionfunctions of νe and νx , respectively. Also, Ce ∼ 1.27 is the collision coefficient for νe and Cx ∼ 0.92
is that for νx [21]. By denoting the effect of neutrino oscillation at a certain time by ∆f ≡7

The energy density of νe is slightly larger than those of νx and νsp after electron annihilation due to the larger
reaction rate of νe with electrons. Therefore, we discriminate among ρνα, std with different flavors.

15

ΓX (sec−1 )
10

100

Neff

3 3.046
No self-int

2

No osciwith osci: θ13 (NO)with osci: θ13 (IO)
with osci: θ12

1

Neff

3 3.046
with self-int

2

No osciwith osci: θ13 (NO)with osci: θ13 (IO)
with osci: θ12

1

Neff

3 3.046
2

with osci: θ12

1
0

No self-int
with self-int

1

10

TRH (MeV)
Fig. 1: Relations between TRH and Neff . The top- and middle panels respectively show the effectof neutrino oscillation in the case without and with neutrino self-interaction, while the bottom
one shows the effect of neutrino self-interaction in the case with (δm212 , θ12 ). The canonical value
Neff = 3.046 [18] is also plotted with the black dotted horizontal line.

fν |with osci −fν |no osci ≡ −∆fνe = ∆fνx , we can evaluate the effect of neutrino oscillation on the total
production rate of neutrinos ∆Rν, tot ≡ Rν, tot |with osci − Rν, tot |no osci = (Rνe |with osci + Rνx |with osci ) −(Rνe |no osci + Rνx |no osci ) = (Rνe |with osci − Rνe |no osci ) + (Rνx |with osci − Rνx |no osci ) ≡ ∆Rνe + ∆Rνx as

164
TRH = 10 MeV

Neff

3

TRH = 5 MeV

2
TRH = 2 MeV

1
TRH = 1 MeV

0

10

1

0.1

0.01

Tγ (MeV)Fig. 2: Time evolution of Neff for each value of TRH . The black solid line is for TRH = 10 MeV,
the blue long-dashed line is for TRH = 5 MeV, the red middle-dashed line is for TRH = 2 MeV,
and the green short-dashed line is for TRH = 1 MeV. Neutrino oscillation with (δm212 , θ12 ) andneutrino self-interaction are considered in the calculation.
follows:∆Rν, tot = ∆Rνe + ∆Rνx ,∼ G2F Tγ5 (−Ce ∆fνe − Cx ∆fνx ) ,= G2F Tγ5 (Ce − Cx )∆f .

(3.4)

As we can see from the expression, the quantity ∆Rν, tot is larger than zero when ∆f = fνe −fνx > 0
which holds if the reheating temperature is sufficiently low for neutrinos to be fully thermalized andthereby fνe > fνx . Consequently, we can see that neutrino oscillation increase the total production
rate of neutrinos unless all neutrinos are completely thermalized. Neutrino self-interaction plays
a role similar to neutrino oscillation. That is, they equilibrate abundances of neutrinos among
themselves and enhance the thermalization of neutrinos in the same way as neutrino oscillation.
We can also see from Fig. 1 that the effect of θ13 is much smaller than that of θ12 or neutrino
self-interaction. This is true independent of the neutrino mass ordering. The relative differences ofeffects among different mixings can be understood as follows. If the vacuum term (Hvac ) dominates
other matter terms (Hmat ) and neutrino oscillation occur adiabatically, the effective transition rateof neutrinos from one flavor to another (i.e., να → νβ where α 6= β) due to neutrino oscillation can

17
be written as [22]
Γtrans =1
sin 2θ Γcoll ,
4

(3.5)

where Γcoll is the collision rate of neutrinos. Therefore, the value of the mixing angle solely
determines how large oscillation happens in this case. As we can see from Eqs. (2.10)-(2.12), thevalue of sin2 θ12 is almost ten times larger than that of sin2 θ13 . That is the reason that the effectof θ12 on the neutrino thermalization is larger than θ13 in the case of normal mass ordering. In
the case of inverted mass ordering, neutrino oscillation proceeds via MSW resonance, and nonadiabatic effects can be important. As for this point, authors in Ref. [9] evaluated the adiabaticityof the MSW resonance and concluded that the non-adiabatic effects are negligible. Therefore, anefficient oscillation should occur when a large population of neutrinos go though the resonanceeven if we adopt the reactor neutrino mixing θ13 . Since the MSW resonance happens at around atemperature of Tc ∼ 5 MeV for neutrino with p = hpi ∼ 3.15 Tγ , we can expect larger oscillation
effects in the case of θ13 (IO) at TRH < Tc , which in fact can be seen in Fig. 1.
Figure 3 shows the time evolution of Neff , and Fig. 3 is the same as Fig. 2, but for the contribution
for each neutrino species Neff, να (α = e, x, sp). Since the final abundance of neutrinos does notdepend on the condition before reheating, we assume that neutrinos have thermal spectra (i.e.,
Fermi-Dirac distributions) at the initial time. We can see from Fig. 2 that the value of Neffdecreases until TRH is realized. This is due to the entropy production from decays of the massive
particles. The value of Neff then increases at Tγ < TRH until neutrinos are decoupled from otherparticles at around a few MeV. This corresponds to the upturn behavior in the evolution of Neff .
In addition, we can see from the evolution of Neff, να in Fig. 3 that neutrino oscillation becomes
effective at around a temperature of a few MeV, whereas neutrino self-interaction becomes effective
at higher temperature. The former is because, as we have discussed in the previous section, neutrino
oscillation with the solar neutrino mixing (δm212 , θ12 ) are effective when the photon temperature
is lower than Tc ∼ 3 MeV. The latter is due to the reason that the reaction rates of neutrino
self-interactions monotonically increase with the photon temperature.
The role of neutrino oscillation and neutrino self-interaction is shown in Fig. 4, where we plot
final energy spectra of neutrinos for the cases with and without neutrino oscillation or neutrino selfinteraction. The energy spectra are evaluated at Tγ ∼ 10−2 MeV which corresponds to the epochwell after electron annihilation. In the case of TRH = 2 MeV, both of these effects decrease thedifference in neutrino abundances. On the other hand, if TRH is large enough (e.g., TRH = 10 MeV),
neutrinos are almost completely thermalized well before decoupling. Therefore, neutrino oscillationNeff, νe

18

No osci, No self-int
No osci, with self-int
with osci, No self-int
with osci, with self-int

1

TRH = 2 MeV

Neff, νx

0
1

Neff, νsp

0
1

0

10

1

0.1

0.01

Tγ (MeV)Fig. 3: Time evolution of Neff, να for TRH = 2 MeV. Neff, νe (top panel) is for νe , Neff, νx (middle
panel) is for νx and Neff, νsp (bottom panel) is for νsp . In the red long-dashed and green solidlines, neutrino oscillation with (δm212 , θ12 ) is taken into account, and neutrino self-interaction isconsidered in the blue middle-dashed and green solid lines. In the figure of Neff, νsp , the black
short-dashed (blue middle-dashed) and red long-dashed (green solid) lines are overlapping.
or neutrino self-interaction plays no role in the final abundance of neutrinos.
Figure 5 shows the dependence of the mean energy of νe , νx and νsp (i.e., ρνα /nνα ) on TRH . The
quantity Rdist on the vertical axis was introduced to measure the distortion in the final energy
spectrum of neutrinos in Ref. [3], and it is defined as
Rdist =

1ρν,
3.15 Tν, eff nν

(3.6)
2

4π
where Tν, eff is the effective temperature of neutrinos Tν, eff = [ 3ζ(3)nν ]1/3 . As we can see from this19

p3 fν ( MeV3 )

5
4
3
2

Equilibrium
No self-int (νe )No self-int (νx )No self-int (νsp )with self-int (νe )with self-int (νx )with self-int (νsp )Equilibrium
No osci (νe )No osci (νx , νsp )with osci (νe )with osci (νx )with osci (νsp )4
3
2

TRH = 10 MeV

10

TRH = 2 MeV

TRH = 2 MeV

105

p3 fν ( MeV3 )Equilibrium
No self-int (νe )No self-int (νx )No self-int (νsp )with self-int (νe )with self-int (νx )with self-int (νsp )Equilibrium
No osci (νe )No osci (νx , νsp )with osci (νe )with osci (νx )with osci (νsp )

0

2

4

p/Tγ

6

8

TRH = 10 MeV

10 0

2

4

6

8

10

p/Tγ

Fig. 4: Final energy spectra of neutrinos. The horizontal axis is the neutrino energy p divided by
the photon temperature Tγ . The vertical axis is the differential energy spectrum of neutrinos. In
the figure, neutrino self-interaction is considered in the left column, whereas neutrino oscillation
with (δm212 , θ12 ) is considered in the right column. We consider (δm212 , θ12 ) in the case withneutrino oscillation in the left column. The thermal spectrum is plotted with the black
dot-dashed line. In the top-left panel, the thin red short-dashed and blue thick short-dashed lines
are overlapping, while the red thin long-dashed and blue thick long-bashed lines are overlappingin the top-right panel. Also, all plots are almost overlapping in the case of TRH = 10 MeV.

definition, Rdist = 1 corresponds to the thermal spectrum, and Rdist > 1 indicates a larger meanenergy. It can be seen from Fig. 5 that the value of Rdist increases as TRH becomes smaller. Thisis because neutrinos are only produced from the annihilation of electrons e− + e+ → να + ν̄α , andneutrinos in the final state therefore have an energy larger than twice the electron mass. Thus,
if neutrinos are mainly produced when the electron mass is not negligible, and the equilibration
process e± + να → e± + να is not effective as in the case of TRH . O(1) MeV, Rdist becomeslarger than unity. In addition, since νe scatter with electrons stronger than νx and νsp due to the20ΓX (sec−1 )Rdist, νsp

Rdist, νx

Rdist, νe

10100

No osci, No self-int
No osci, with self-int
with osci, No self-int
with osci, with self-int

3.0
2.0
1.0
3.0
2.0
1.0
3.0
2.0
1.0
1

10

TRH (MeV)
Fig. 5: Dependence of the mean energy of νe (top panel), νx (middle panel) and νsp (bottom
panel) on TRH . The vertical axis is the distortion parameter Rdist for each neutrino species. Rdist= 1 corresponds to the thermal spectrum. In the middle panel, the blue long-dashed and green
solid lines are overlapping. Also, in the bottom panel, the black dot-dashed (red short-dashed)
and blue long-dashed (green solid) lines are overlapping. We consider (δm212 , θ12 ) in the case withneutrino oscillation.

charged-current interaction, the energy distribution of νe is closer to the thermal spectrum. That is
the reason that the relation Rdist, νe < Rdist, νx , Rdist, νe < Rdist, νsp holds for a sufficiently small TRH .
Furthermore, we can see from Fig. 5 that both neutrino oscillation and neutrino self-interaction
increase Rdist, νe , while decrease Rdist, νx . This is because neutrino oscillation and neutrino selfinteraction equilibrate the neutrino abundances of different flavors as shown in Fig. 4. The reason

21
is as follows: If the final distribution function of neutrinos is changed by a factor of κ (i.e., fν 7→ κfνwhere κ < 1 for νe and κ > 1 for νx ) due to neutrino oscillation or neutrino self-interaction, thenthe distortion parameter should be also modified by
−1ρν /nν
κρν /κnνe= κ−1/3 ,Rdist /Rdist =3.15(κ1/3 Tν, eff ) 3.15Tν, eff(3.7)

edist and Rdist are distortion parameters for the cases with and without effects of neutrino
where R

oscillation or neutrino self-interaction, respectively. Therefore, these effects increase Rdist, νe anddecrease Rdist, νx as long as neutrino oscillation or neutrino self-interaction is effective. In addition,
since the reaction rate of neutrino self-interaction strongly depends on the number density ofneutrinos, its effect on Rdist becomes small more rapidly than neutrino oscillation as TRH decreases.

These effects on Rdist can be estimated by comparing neutrino distribution functions in Fig. 4 andis consistent with the results in Fig. 5. On the other hand, since νsp does not mix with other flavorof neutrinos, they are only affected by neutrino self-interaction.
In the next section, we discuss the light element abundances created in the process of BBN
taking our computed neutrino thermalization into account.

22
IV.

BIG BANG NUCLEOSYNTHESISAs mentioned in Sec. II, incomplete thermalization of neutrinos affects the dynamics of thestandard BBN. In this section, we explain the role of neutrinos in the production process of lightelements and show our results of BBN obtained by assuming TRH ∼ O(1) MeV.
A.

Formulation of BBN

We have seen in the previous section that the late-time entropy production due to decays of X
induces the incomplete thermalization of neutrinos before decoupling. Since neutrinos take part in
the weak reaction processes,n ↔ p + e− + ν̄e ,(4.1)

e+ + n ↔ p + ν̄e ,(4.2)

νe + n ↔ p + e− ,(4.3)

which interchange ambient neutrons and protons with each other, non-thermal spectra of neutrinos
significantly change the freeze-out value of the neutron-to-proton ratio (n/p)f ≡ (nn /np )T =Tf wherenn and np are the number density of neutrons and protons, respectively, whereas Tf is the freezeout temperature of the processes (4.1)-(4.3). As described later in this section, the theoreticalvalues of light element abundances are very sensitive to the neutron-to-proton ratio before BBN.
Therefore, theoretical predictions of the standard BBN should be modified in the Universe withsmall TRH . Since the predictions of standard BBN is well consistent with the observational values,
we can constrain TRH by requiring that the late-time entropy production does not spoil the current
success of the standard BBN.
In the case where the massive particles have a hadronic branching ratio, there are additional
neutron-proton interchanging processes other than (4.1)-(4.3) via strong interactions caused byinjected hadrons N + H ↔ N ′ + H ′ where N and N ′ are nucleons, and H and H ′ are mesons orbaryons. If the hadronic branching ratio is large enough, the hadronic processes dominantly affectthe neutron-to-proton ratio, which result in different light element abundances compared to thecase of the 100% radiative decays of the massive particles [3]. In the current study, we consider the
hadronic processes involving pions (π ± ) and nucleons (n, n̄, p, p̄) which are injected from hadronic
decays of the massive particles. The energetic hadrons produced in the decay of the massive particles are instantaneously stopped by Coulomb scattering with background electrons/positrons or23
inverse-Compton like scattering with background photons [23–25]. Therefore, the hadrons affecting
neutron-proton inter-conversions are thermalized, and we can use thermal cross sections for the
calculation. As for the hadronic cross sections, we adopt those given in Table.1 of Ref. [23] for themean values and assume 30% experimental error in each cross section for a conservative treatment(see also Refs. [3, 26]).
In order to follow the evolution of light element abundances, we solve the Boltzmann equationsof light elements using the Kawano code [27]. Since some of the nuclear reaction rates in the code
are already outdated, we replace them with the latest ones (see Ref. [28] for more information).
In addition, we rewrite some equations in the code to allow for the late-time entropy production
accompanied by the decays of X. Moreover, since the free neutron decay (i.e., the forward processof (4.1)) continues even after the other weak processes of (4.1)-(4.3) decoupled at Tγ ∼ Tf , thevalue of the neutron-to-proton ratio just before BBN depends on the lifetime of neutrons (seee.g., [29]). In the current study, we use the value of the neutron lifetime τn reported in Ref. [1]:τn = 880.2 ± 1.0 sec (68% C.L.) .

(4.4)

As for the observational values of light elements, we adopt the primordial mass fraction ofhelium 4 He, Yp , reported in Ref. [30]:Yp = 0.2449 ± 0.0040 (68% C.L.) ,(4.5)

whereas for the observational value of primordial abundance of deuterium D, we adopt the latest
value reported in Ref. [31]:D/H = (2.545 ± 0.025) × 10−5 (68% C.L.) .

B.

(4.6)

Results of BBN: Radiative decay

First we show the results of radiative decay, i.e., the hadronic branching ratio Br = 0. In thiscase, photons and charged leptons emitted from the decay of X are instantaneously thermalized
via electromagnetic force, and results of neutrino thermalization and BBN are independent of themass of X, mX . In Fig. 6 and Fig. 7, the relation between TRH and D/H and Yp are shown,respectively. We assume the 100% radiative decays of the massive particles in these figures. The
baryon-to-photon ratio ηB is the only free parameter in the standard BBN. In the low-reheatingtemperature Universe, a baryon number is diluted by the entropy production due to the decays

24

ΓX (sec−1 )10

No self-int

Yp

0.260
0.240

2 σ obs boundNo osciwith osci: θ13 (NO)with osci: θ13 (IO)
with osci: θ12

0.220
0.200

with self-int

Yp

0.260
0.240

2 σ obs boundNo osciwith osci: θ13 (NO)with osci: θ13 (IO)
with osci: θ12

0.220
0.200

with osci: θ12

0.260

Yp

100

0.240
0.220

2 σ obs boundNo self-int
with self-int

0.200
1

10

TRH (MeV)
Fig. 6: Relations between TRH and Yp in the case of the 100% radiative decays of X. We adopt
ηB = 6.13 × 10−10 in the figure. The top- and middle panels show the effect of neutrino oscillation
for the cases without and with neutrino self-interaction, respectively. The black dot-dashed line isfor the case without neutrino oscillation, the red short-dashed line is for the case with θ13 (NO),the blue long-dashed line is for the case with θ13 (IO), and the green solid line is for the case with
θ12 . The bottom panel shows the effect of neutrino self-interaction when we consider neutrino
oscillation with (δm212 , θ12 ). The black dashed- and solid lines are for the cases without and withneutrino self-interaction, respectively. The gray-shaded region corresponds to the 2σ
observational bound.

of the massive particles, and hence ηB is decreased by many orders of magnitude. Therefore, we
set the large initial value of ηB so that the final value of ηB is consistent with observations of lightelements. To plot Fig. 6 and Fig. 7, we fix the final value of ηB to the median value reported by25

ΓX (sec−1 )10100

D/H (×105 )

2.6
2.4

No self-int

2.2

2 σ obs boundNo osciwith osci: θ13 (NO)with osci: θ13 (IO)
with osci: θ12

2.0
1.8

D/H (×105 )

2.6
2.4

with self-int

2.2

2 σ obs boundNo osciwith osci: θ13 (NO)with osci: θ13 (IO)
with osci: θ12

2.0
1.8

D/H (×105 )

2.6
2.4

with osci: θ12

2.2
2 σ obs boundNo self-int
with self-int

2.0
1.8
1

10

TRH (MeV)
Fig. 7: Same as Fig. 6, but for D/H.
Planck collaboration [32]:ηB = 6.13 × 10−10 .

(4.7)

Since almost all neutrons are processed into 4 He which is the most stable among light elements, theprimordial mass fraction of 4 He can be written as Yp ≡ ρ4 He /ρB ∼ 2/{1 + (n/p)−1
BBN } ∼ 0.25, where(n/p)BBN ≡ (n/p)f e−t/τn is the neutron-to-proton ratio just before deuterium bottleneck opens(i.e., Tγ ∼ 0.08 MeV and t ∼ 200 sec), and the last approximation holds in the standard big-bang
cosmology where (n/p)BBN ∼ 1/7 [29]. Therefore, the value of (n/p)BBN almost entirely determines

the final abundance of 4 He. As for the final abundance of D, the value of Neff is also importantbecause it is related to the Hubble parameter (see Eqs. (2.32), (2.34) and (3.1)) and determines
when each light element departs from the nuclear statistical equilibrium [33, 34]. 8 In addition,
8

We can intuitively understand the dependence of the deuterium abundance on the expansion rate of the Universein the BBN epoch by comparing the binding energy of D and 4 He; i.e., BD ∼ 2.22 MeV and B4 He ∼ 28.3 MeV.

26
as we can see from Fig. 6 and Fig. 7, the influences of neutrino oscillation and self-interaction onlight element abundances are similar, and both of these effects increase Yp and D/H. In order tounderstand the numerical results on light element abundances, next we focus on the dynamics ofthe freeze-out of the neutron-to-proton ratio.
Since nucleons are always non-relativistic, (n/p)f can be expressed with the freeze-out temperature if TRH is MeV scale as (n/p)f ∼ exp(−Q/Tf ) where Q ≡ mn − mp ∼ 1.3 MeV is the massdifference of nucleons. We note that Tf is determined by the relative values of the neutron-protoninter-converting weak reaction rates Γnp and the Hubble parameter H and is roughly given byΓnp (Tf )/H(Tf ) ∼ 1. In the low-reheating-temperature Universe, the total energy density is smaller
than that in the standard big-bang cosmology under the same photon temperature due to theincomplete thermalization of neutrinos. 9 Therefore, the expansion rate of the Universe is also
small in the scenario, and this effect delays the decoupling of the processes (4.1)-(4.3) and thereby
increases (n/p)f .
The influence is not only in the Hubble parameter H but also in the reaction rates Γnp . Specifically, the reaction rate of the processes (4.1)-(4.3) can be written as [3]#"Z Q−me′2pp
1 − fνe (p′ ) ,(p′ − Q)2 − m2e (Q − p′ )dp′Γn→pe−ν̄e = K′
1 + e(p −Q)/Tγ
0
#
"
Z ∞
′2
p

p
1 − fνe (p′ ) ,dp′
Γne+ →pν̄e = K(p′ − Q)2 − m2e (p′ − Q) (pνe −Q )/Tγ + 1
e
Q+me
"
#
Z ∞
p
p′ 2
′
′
′
′
2
2
dp
Γnνe →pe− = K(p + Q) − me (p + Q)fνe (p ) ,′1 + e−(p +Q)/Tγ0"#
Z Q−me′2pp′′′(p′ − Q)2 − m2e (Q − p ) −(pdpfνe (p ) ,Γpe− ν̄e →n = Ke νe −Q )/Tγ + 1
0
"
#
Z ∞
′2
p

p
′
′
′
dp
Γpe− →nνe = K(p′ + Q)2 − m2e (Q + p ) (p1 − fνe (p ) ,e νe +Q )/Tγ + 1
0
"
#
Z ∞
′2
p
p
dp′
Γpν̄e →ne+ = K(p′ − Q)2 − m2e (Q − p′ )fνe (p′ ) ,1 + e−(pνe −Q )/Tγ
Q+me

where me is the electron mass, and K ∼ (1.636τn )−1 is a normalization factor whose value isdetermined by the neutron lifetime τn . Of these reaction rates, some depend on fνe and others
Since B4 He is much larger than BD , most of deuterium created in the process p + n → D + γ are processed into 4 He.
However, if the expansion rate is large, a process which is responsible for the 4 He production decouples earlier,and there remains a large abundance of deuterium. That is the reason that a large expansion rate in the BBN
9epoch leads to a large deuterium abundance and vice versa.
We recall the reader that Tγ determines when light element abundances are to be created since the reaction
process p + n → D + γ responsible for the deuterium production is the first step of BBN, and its backward reaction
rate depends on Tγ . For this reason, the values of Γnp and H should not be characterized by the cosmic time
but rather by Tγ . That is the reason that a larger value of Neff ∝ ρν /ργ leads to a larger expansion rate of the
Universe at the epoch of BBN.

27
on 1 − fνe . In the low-reheating-temperature Universe, the neutrino abundance of each flavor issmaller than the case of the standard big-bang cosmology. Therefore, by denoting the reductionof fνe due to the incomplete thermalization of neutrinos by ∆fνe , the following relation holds for
sufficiently small TRH (< Tdec ∼ O(1) MeV):

∆fνe
∆fνe
≪
(1 − fνe )f νe

for fνe ≪ 0.5 .

(4.8)

As a result, with such a small value of TRH , the total reaction rate Γnp ≡ Γn→pe− ν̄e + Γne+ →pν̄e +Γnνe →pe− + Γpe− ν̄e + Γpe− →nνe + Γpν̄e→ne+ becomes smaller than that of the standard big-bang
cosmology as written in Ref. [3]. Therefore, this effect accelerates the decoupling of the processes (4.1)-(4.3) and thereby decreases (n/p)f . Consequently, the relative magnitude of these two
opposite contributions determine the net effect of incomplete thermalization of neutrinos on (n/p)f .
As described in the previous section, neutrino oscillation and self-interaction slightly enhance
the neutrino thermalization and increase the total energy density of neutrinos and hence Neff . As aresult, the Hubble expansion rate increases due to these effects. In addition, since the νe abundanceis decreased by the conversion νe → νx , and only νe take part in the reaction processes (4.1)-(4.3),
Γnp decreases by considering these effect. Therefore, neutrino oscillation and self-interaction always
play a role in delaying the freeze-out of neutron-to-proton ratio and increasing Tf and (n/p)f , whichleads to larger values of Yp and D/H.
As is the case for Neff (see Fig. 1), we can also see from Fig. 6 and Fig. 7 that the impactof the solar neutrino mixing (δm212 , θ12 ) is much larger than that of the reactor neutrino mixing(δm213 , θ13 ) independent of the neutrino mass ordering. Therefore, the effective two-flavor mixing
with (δm212 , θ12 ) gives a good approximation to the full three-flavor neutrino mixings. For thisreason, we hereafter only consider (δm212 , θ12 ) in the case with neutrino oscillation.
To obtain the observational constraint on TRH , we perform a Monte-Carlo calculation of BBN
and χ2 analysis at each point on the grids of ηB and TRH assuming observational values for Yp(Ref. [30]) and D/H (Ref. [31]). 10 In the Monte-Carlo calculation, we assume that the reaction
rates in the standard BBN, the hadronic reaction rates and the neutron lifetime follow Gaussian
distribution and propagate their reported errors to obtain theoretical uncertainties on the lightelement abundances. Since an allowed region is defined by a parameter space where theoretical10

As written in e.g., [6, 35], it is technically incorrect to adopt the CMB bound ηB = (6.13 ± 0.04) × 10−10 reportedby the Planck collaboration [32] as a prior of BBN because the recombination process depends on the values of
Neff and Yp , and there are correlations between ηB and these quantities. In other words, CMB is not independent
from the neutrino thermalization and BBN. In Ref. [32], they adopt the canonical value Neff = 3.046 [18] and Yp
calculated by assuming the standard BBN, which are not necessarily realized in the low-reheating-temperature
Universe.

28abundances of light elements explain each observational value, we give the lower bound on TRH
combining χ2 values of both D/H and Yp : 11
χ2D/H + Yp ≡ χ2D/H + χ2Yp ={(D/H)th (ηB , TRH ) − (D/H)obs }2 {Yp, th (ηB , TRH ) − Yp, obs }2+ 2, (4.9)2
2
σD,
σYp , th (ηB , TRH ) + σY2 p , obs
th (ηB , TRH ) + σD, obs

where χ2D/H and χ2Yp are χ2 values of D/H and Yp , respectively. Also, σi, th and σi, obs where i =D/H and Yp are respectively the theoretical and observational 1σ variance of each light element
abundance.
Figure 8 shows the allowed region in the plane of ηB and TRH in the case of the 100% radiative
decays. In the current study, we assume that χ2D/H and χ2Yp follow Gaussian distribution. In thiscase, we can find the lower bound at 95% C.L. on TRH by requiring χ2D/H + Yp (ηB , TRH ) < 5.991:
TRH & 1.8 MeV ,(4.10)

in the case with both neutrino oscillation and self-interaction. Also, we depict in Fig. 9 thecomparison between cases with and without neutrino oscillation or neutrino self-interaction. As canbe seen from Fig. 9, we find TRH & 1.5 MeV in the case with neutrino oscillation and without selfinteraction, whereas TRH & 0.6 MeV in the case without neutrino oscillation and with neutrino selfinteraction. The BBN bound in the case with neutrino oscillation or self-interaction is tighter thanthat in the case without them. This is because, as we can see from Fig. 7, neutrino oscillation andself-interaction increase the value of Yp , and the discrepancy between theoretical and observational
values becomes large compared to the case without these effects.

C.

Results of BBN: Hadronic decay

As described above, if the massive particles have a branching ratio into hadrons, the constraint
on TRH imposed by BBN can be modified compared to when the decays of X are fully radiative(i.e., Br = 0). The effects of hadronic decays on light element abundances are shown in Fig. 10
where we plot the dependence of D and 4 He abundances on TRH for each value of mX and Br.
The case of Br = 0 in the figure corresponds to the 100% radiative decays of X which is plotted
for reference. In the figure, we assume that the massive particles have a non-negligible branchingratio into u ū quark pairs to calculate the number of hadrons produced in the decays of X with
11

There remains a long-standing problem in the standard BBN that the theoretical prediction of the 7 Li abundance isapproximately three times larger than that of the observational value if we input the value of the baryon-to-photon
ratio from CMB into the calculation of BBN (see e.g., [1]). Therefore, we refrain from using the 7 Li abundance toconstrain TRH in the current study.

29

10

10

ΓX (sec−1)

TRH (MeV)

100

95% C.L. (Yp + D/H)
68% C.L. (Yp + D/H)

1
5.0

6.0

7.0

ηB (×1010)

100

10

1
5.0

95% C.L. (D/H only)
68% C.L. (D/H only)
6.0

95% C.L. (Yp only)
68% C.L. (Yp only)
7.0 5.0ηB (×1010)

6.0

ΓX (sec−1)

TRH (MeV)

10

7.0

ηB (×1010)

Fig. 8: Allowed region in the (ηB , TRH ) plane in the case of the 100% radiative decays of X. The
95% (68%) C.L. contour is plotted with the blue solid (red dashed) line. The top panel shows theallowed region in terms of both Yp and D/H, whereas the bottom-left and bottom-right panels
show those of D/H and of Yp , respectively. Neutrino oscillation and neutrino self-interaction are
considered in the calculation.
Pythia 8.2 code. 12 Since the neutron-proton inter-converting processes caused by the injected
hadrons increase the neutron-to-proton ratio, both the 4 He and D abundances, which increase12

We have checked that the BBN bound does not depends on the quark flavor emitted from the massive particles if
the mass of the massive particles is much larger than the total mass of emitted quarks (i.e., mX >> mqα wheremqα is the quark mass of particular flavor α).

30

10
100

10
with osci, with self-int
with osci, No self-int
No osci, with self-int

1

4.0

5.0

6.0

ΓX (sec−1)

TRH (MeV)

95% C.L. (Yp + D/H)

7.0

ηB (×1010)
Fig. 9: Comparison of allowed regions in the (ηB , TRH ) plane in the cases with or without
neutrino oscillation and self-interaction. We assume the 100% radiative decays of X. The 95%
contour with the green solid line is for the case with both neutrino oscillation and self-interaction,
one with the red long-dashed line is for the case only with neutrino oscillation, and one with theshort-dashed line is for the case only with neutrino self-interaction.

with (n/p)f , get larger than those in the standard BBN case. In addition, we find that the effectof hadronic decays is large for a large Br or a small mX . Comparing the cases of mX = 100 TeV,
Br = 0.001 (red thin dashed) and mX = 10 GeV, Br = 1 (blue thick solid) in Fig. 10, we can seethat the discrepancy of D/H or Yp between these cases are of the same order or much larger thanthat of the 2σ observational error if TRH is a few MeV.
To understand the reason, we define the comoving variable for the initial abundance of themassive particles YX = nX /s where nX is the number density of the massive particles X, and s is
the total entropy density of the Universe. If we assume that X dominates the total energy at theinitial time and most of it is transferred to radiation components before the reheating process iscompleted, we can write the initial value of YX as follows:2( π g∗ T 4 )/mX
3 TRH
nX
∼ 30 2π2 RH
,
∼
YX =
3
∗
s
4 mX
gs TRH(4.11)

45

where g∗ and gs∗ are relativistic degrees of freedom defined by energy and entropy density respectively. In the standard big-bang cosmology, g ∗ ∼ gs∗ holds before electron-positron annihilation

31

−1

ΓX (sec )10100

D/H (×105)

3.2
2.8
2.4

2 σ obs boundBr = 02.0

Br = 1.0 (mX = 10 GeV)
Br = 0.001 (mX = 10 GeV)
Br = 1.0 (mX = 100 TeV)
Br = 0.001 (mX = 100 TeV)

1.6
0.380

Yp

0.340
0.300
0.260
0.220
1

10

TRH (MeV)
Fig. 10: D/H and Yp as a function of TRH in the case with hadronic decays of X. he blue thicksolid- and dashed lines are for the case of mX = 10 GeV, while the red thin solid- and dashedlines are for the case of mX = 100 TeV. In addition, the blue- and red solid lines are for the caseof Br = 1, while the blue- and red dashed lines are for the case of Br = 0.001. For comparison, we
also plot the case of Br = 0 with the black dot-dashed line. In the figure, we consider bothneutrino oscillation and self-interaction.

sets in. As we can see from the above expression, YX gets larger for smaller mX . In addition, thenumber of hadrons emitted from the decays of X is almost proportional to m0.4
X (see Ref. [25]), andtherefore the total number of hadrons emitted from X is almost proportional to m−0.6
X . Since theenergetic hadrons instantaneously lose their energy and are thermalized with background particles
before inter-converting ambient neutrons and protons, the number of emitted hadrons determine
the magnitude of the hadronic-decay effect on BBN. Therefore, the influence of hadronic decays
on BBN should be stronger for smaller mX .
For the purpose of showing the effects of neutrino oscillation and self-interaction on the light32

ΓX (sec−1 )ΓX (sec−1 )D/H (×105)

10
2.8
2.4

10

100

mX = 10 GeV

mX = 100 TeV

Br = 1.0

Br = 0.001

2σ obs boundwith osci, with self-int
with osci, No self-int
No osci, with self-int
No osci, No self-int

2.0

0.260

Yp

100

mX = 10 GeV

mX = 100 TeV

Br = 1.0

Br = 0.001

0.250
0.240
1

10 1

TRH (MeV)10

TRH (MeV)

Fig. 11: Effects of neutrino oscillation and self-interaction on D/H and Yp in the case withhadronic decays of X. We adopt ηB = 6.13 × 10−10 in the figure. The red solid- and blue
long-dashed lines are for the case with neutrino oscillation, while the green short-dashed andblack dot-dashed line are for the case with neutrino self-interaction.

element abundances in the case of hadronic decays, we plot in Fig. 11 the dependence of Yp andD/H on TRH for (mX , Br) = (10 GeV, 1.0) and (100 TeV, 0.001) where we expect large and small
effects of hadronic decays, respectively. As we can see from Fig. 11, if TRH is a few MeV, neutrino
oscillation and self-interaction affect light element abundances at the level of O(10)% for D/H andO(1)% for Yp when mX = 10 GeV and Br = 1.0, whereas the correction is O(10)% for both casesof D and Yp when mX = 100 TeV and Br = 0.001. Since we give the observational bound on TRH
by summing up the χ2 values of D/H and Yp , the constraint on TRH should be changed by O(1)%when mX = 10 GeV and Br = 1.0 and by O(10)% when mX = 100 TeV and Br = 0.001.
We show in Fig. 12 the allowed region in the same plane as Fig. 8, but in the case when hadronic
decays are included. In the figure, we show four representative cases of (mX , Br) = (10 GeV, 1.0),(10 GeV, 0.001), (100 TeV, 1.0) and (100 TeV, 0.001).

33

100

95% C.L.
68% C.L.

10mX = 100 TeV, Br = 1.0

mX = 10 GeV, Br = 1.0

ΓX (sec−1)

TRH (MeV)

10

1
100

10
mX = 10 GeV, Br = 0.001

1
5.0

6.0

mX = 100 TeV, Br = 0.001

7.0 5.0ηB (×1010)

6.0

ΓX (sec−1)

TRH (MeV)

10

7.0

ηB (×1010)

Fig. 12: Same as Fig. 8, but for the case of hadronic decays of the massive particles.
A possible minimum value of the reheating temperature TRH, min in terms of BBN is shown in
Fig. 13 as a function of mX . We can see from the figure that the BBN bound is tighter in the caseof a small mX or a large Br. Consequently, we obtain the lower bound on TRH at 95% C.L.:
TRH & 4.1 − 4.9 MeV

for mX = 10 GeV − 100 TeV ,(4.12)

when the hadronic branching ratio Br = 1.0, whereas
TRH & 2.1 − 3.7 MeV

for mX = 10 GeV − 100 TeV ,(4.13)

when Br = 0.001 in the case with both neutrino oscillation and neutrino self-interaction. In
addition, we find neutrino oscillation and neutrino self-interaction can change the value of TRH, minat the level of O(1)% for most of the range of mX in the case of hadronic decays.

34

TRH, min (MeV)

6
Br = 1.0

4 Br = 0.001
2 Radiative decay

00.01

95% C.L. (Yp + D/H)

0.1

1

10

100

mX (TeV)Fig. 13: Lower bound on the reheating temperature TRH, min at 95% C.L. as a function of mX in
the case with both neutrino oscillation and neutrino self-interaction. We plot the results in thecases of Br = 0.001 (blue dashed) and 1.0 (red solid). The case of the 100% radiative decay isalso plotted with the black dashed line.

35
V.

CONCLUSION

In this paper, we have investigated the possibility that the reheating temperature of the Universeis O(1) MeV motivated by long-lived massive particles which often appear in the particle physics
theory beyond the standard model and induce a late-time entropy production by their decays. In
this scenario, neutrinos are not necessarily thermalized well before the beginning of BBN. Hence,the expansion rate of the Universe and weak reaction processes are significantly altered, whichchanges the freezeout value of the neutron to proton ratio. We have calculated the thermalization
process of neutrinos including effects of both neutrino oscillation and neutrino self-interaction(Figs. 1-3), and obtained a lower bound on the reheating temperature TRH & 1.8 MeV (95% C.L.)(Fig. 8) in the case of the 100% radiative decay.
On the other hand, if the massive particles also decay into hadrons, there is an additional effecton BBN via inter-conversion of ambient neutron and proton through the scatterings of the hadrons.
In this case, the constraint becomes tighter than that of the 100% radiative decay (Fig. 10). Then,
we obtained the lower bound TRH & 2 MeV – 5 MeV (95% C.L.) depending on the masses of themassive particles (10 GeV – 100 TeV) and the hadronic branching ratio (Figs. 12 - 13).
In addition, we found that neutrino oscillation and neutrino self-interaction increase the efficiency of neutrino thermalization (Figs. 1 and 3) and decrease the exchange rate between neutronsand protons, thereby enhancing the theoretically expected abundances of helium, Yp (Fig. 6), anddeuterium, D/H (Fig. 7). These effects increase the minimum value of the reheating temperature
at the level of O(10)% in the case of the 100% radiative decays (Fig. 8) and O(1)% in most casesof hadronic decays (Fig. 11).
Finally, let us comment on the future prospects of this study. This time, we only focused on
BBN to constrain TRH . On the other hand, as described in Sec. I, CMB and LSS also depend on
the expansion rate in the Universe, and therefore have a sensitivity to the neutrino thermalization.
In addition, the recombination history depends on the 4 He abundance which is strongly affected by
the hadronic decay effects. Therefore, theoretical results of CMB and LSS should be different fromthose in the case of radiative decay. We will discuss observational constraints on TRH from CMB
and LSS in addition to BBN assuming hadronic decays of the massive particles in a forthcoming
paper [36].

36
ACKNOWLEDGMENTS

Numerical computations were carried out on PC clusters at Center for Computational Astrophysics, National Astronomical Observatory of Japan (NAOJ) and Computing Research Center, High Energy Accelerator Research Organization (KEK). KK is supported by JSPS KAKENHI Grants No. JP17H01131, MEXT Grant-in-Aid for Scientific Research on Innovative Areas
Nos. JP15H05889, JP18H04594, JP19H05114, and by WPI, MEXT, Japan. The work of RSLH,
TT, and SH is supported by the Villum Foundation.

37[1] M. Tanabashi, et al. Review of Particle Physics. Phys. Rev., Vol. D98, No. 3, p. 030001, 2018.
[2] M. Kawasaki, K. Kohri, and N. Sugiyama. Cosmological constraints on late time entropy production.
Phys. Rev. Lett., Vol. 82, p. 4168, 1999.
[3] M. Kawasaki, K. Kohri, and N. Sugiyama. MeV scale reheating temperature and thermalization ofneutrino background. Phys. Rev., Vol. D62, p. 023506, 2000.
[4] K. Ichikawa, M. Kawasaki, and F. Takahashi. The Oscillation effects on thermalization of the neutrinosin the Universe with low reheating temperature. Phys. Rev., Vol. D72, p. 043522, 2005.
[5] P. F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor, and O. Pisanti. Bounds on very low
reheating scenarios after Planck. Phys. Rev., Vol. D92, No. 12, p. 123534, 2015.
[6] S. Hannestad. What is the lowest possible reheating temperature? Phys. Rev., Vol. D70, p. 043506,
2004.
[7] K. Ichikawa, M. Kawasaki, and F. Takahashi. Constraint on the Effective Number of Neutrino Species
from the WMAP and SDSS LRG Power Spectra. JCAP, Vol. 0705, p. 007, 2007.
[8] F. De Bernardis, L. Pagano, and A. Melchiorri. New constraints on the reheating temperature of theuniverse after WMAP-5. Astropart. Phys., Vol. 30, pp. 192–195, 2008.
[9] L. Johns, M. Mina, V. Cirigliano, M. W. Paris, and G. M. Fuller. Neutrino flavor transformation in
the lepton-asymmetric universe. Phys. Rev., Vol. D94, No. 8, p. 083505, 2016.
[10] P. F. de Salas, D. V. Forero, C. A. Ternes, M. Tortola, and J. W. F. Valle. Status of neutrino oscillations
2018: 3σ hint for normal mass ordering and improved CP sensitivity. Phys. Lett., Vol. B782, pp. 633–
640, 2018.
[11] B. H. J. McKellar and M. J. Thomson. Oscillating doublet neutrinos in the early universe. Phys. Rev.,
Vol. D49, pp. 2710–2728, 1994.
[12] G. Sigl and G. Raffelt. General kinetic description of relativistic mixed neutrinos. Nucl. Phys., Vol.
B406, pp. 423–451, 1993.
[13] S. Hannestad, R. S. Hansen, T. Tram, and Y. Y. Y. Wong. Active-sterile neutrino oscillations in theearly Universe with full collision terms. JCAP, Vol. 1508, No. 08, p. 019, 2015.
[14] L. Wolfenstein. Effects of Matter on Neutrino Oscillations. AIP Conf. Proc., Vol. 52, pp. 108–112,
1979.
[15] S. P. Mikheyev and A. Y. Smirnov. Resonant neutrino oscillations in matter. Prog. Part. Nucl. Phys.,Vol. 23, pp. 41–136, 1989.
[16] S. Hannestad, I. Tamborra, and T. Tram. Thermalisation of light sterile neutrinos in the early universe.
JCAP, Vol. 1207, p. 025, 2012.
[17] S. Hannestad, R. S. Hansen, and T. Tram. Can active-sterile neutrino oscillations lead to chaotic
behavior of the cosmological lepton asymmetry? JCAP, Vol. 1304, p. 032, 2013.
[18] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico. Relic neutrino decoupling

38
including flavor oscillations. Nucl. Phys., Vol. B729, pp. 221–234, 2005.
[19] P. F. de Salas and S. Pastor. Relic neutrino decoupling with flavour oscillations revisited. JCAP, Vol.
1607, No. 07, p. 051, 2016.
[20] N. F. Bell, R. R. Volkas, and Y. Y. Y. Wong. Relic neutrino asymmetry evolution from first principles.
Phys. Rev., Vol. D59, p. 113001, 1999.
[21] K. Enqvist, K. Kainulainen, and M. J. Thomson. Stringent cosmological bounds on inert neutrino
mixing. Nucl. Phys., Vol. B373, pp. 498–528, 1992.
[22] R. Foot and R. R. Volkas. The Exact parity symmetric model and big bang nucleosynthesis. Astropart.
Phys., Vol. 7, pp. 283–295, 1997.
[23] M. H. Reno and D. Seckel. Primordial Nucleosynthesis: The Effects of Injecting Hadrons. Phys. Rev.,
Vol. D37, p. 3441, 1988.
[24] K. Kohri. Primordial nucleosynthesis and hadronic decay of a massive particle with a relatively short
lifetime. Phys. Rev., Vol. D64, p. 043515, 2001.
[25] M. Kawasaki, K. Kohri, and T. Moroi. Big-Bang nucleosynthesis and hadronic decay of long-lived
massive particles. Phys. Rev., Vol. D71, p. 083502, 2005.
[26] M. Pospelov and J. Pradler. Metastable GeV-scale particles as a solution to the cosmological lithiumproblem. Phys. Rev., Vol. D82, p. 103514, 2010.
[27] L. Kawano. Let’s go: Early universe. 2. Primordial nucleosynthesis: The Computer way. 1992.
[28] M. Kawasaki, K. Kohri, T. Moroi, and Y. Takaesu. Revisiting Big-Bang Nucleosynthesis Constraints
on Long-Lived Decaying Particles. Phys. Rev., Vol. D97, No. 2, p. 023502, 2018.
[29] G. Steigman. Primordial Nucleosynthesis in the Precision Cosmology Era. Ann. Rev. Nucl. Part. Sci.,
Vol. 57, pp. 463–491, 2007.
[30] E. Aver, K. A. Olive, and E. D. Skillman. The effects of He I λ10830 on helium abundance determinations. JCAP, Vol. 1507, No. 07, p. 011, 2015.
[31] E. O. Zavarygin, J. K. Webb, S. Riemer-Srensen, and V. Dumont. Primordial deuterium abundance atzabs = 2:504 towards Q1009+2956. J. Phys. Conf. Ser., Vol. 1038, No. 1, p. 012012, 2018.
[32] N. Aghanim, et al. Planck 2018 results. VI. Cosmological parameters. 2018.
[33] R. Esmailzadeh, G. D. Starkman, and S. Dimopoulos. Primordial nucleosynthesis without a computer.
Astrophys. J., Vol. 378, pp. 504–518, 1991.
[34] M. S. Smith, L. H. Kawano, and R. A. Malaney. Experimental, computational, and observational
analysis of primordial nucleosynthesis. Astrophys. J. Suppl., Vol. 85, pp. 219–247, 1993.
[35] J. Hamann, S. Hannestad, G. G. Raffelt, and Y. Y. Y. Wong. Sterile neutrinos with eV masses incosmology: How disfavoured exactly? JCAP, Vol. 1109, p. 034, 2011.
[36] T. Hasegawa, N. Hiroshima, K. Kohri, S. Wang, R. S. Hansen, T. Tram, and S. Hannestad. In
preparation, 2019.