Precision Early Universe Thermodynamics made simple: N eff and Neutrino Decoupling in the Standard Model and beyond
KKCL-2019-85
Prepared for submission to JCAP
Precision Early UniverseThermodynamics made simple: N eff and Neutrino Decoupling in theStandard Model and beyond Miguel Escudero Abenza
Theoretical Particle Physics and Cosmology GroupDepartment of Physics, King’s College London, Strand, London WC2R 2LS, UKE-mail: [email protected]
Abstract.
Precision measurements of the number of effective relativistic neutrino species andthe primordial element abundances require accurate theoretical predictions for early Universeobservables in the Standard Model and beyond. Given the complexity of accurately modellingthe thermal history of the early Universe, in this work, we extend a previous method presentedby the author in [1] to obtain simple, fast and accurate early Universe thermodynamics.The method is based upon the approximation that all relevant species can be described bythermal equilibrium distribution functions characterized by a temperature and a chemicalpotential. We apply the method to neutrino decoupling in the Standard Model and find N SMeff = 3 .
045 – a result in excellent agreement with previous state-of-the-art calculations.We apply the method to study the thermal history of the Universe in the presence of a verylight (1 eV < m φ < λ (cid:46) − ) neutrinophilic scalar. We find ourresults to be in excellent agreement with the solution to the exact Liouville equation. Finally,we release a code: NUDEC BSM (available in both Mathematica and
Python formats), withwhich neutrino decoupling can be accurately and efficiently solved in the Standard Modeland beyond: https://github.com/MiguelEA/nudec_BSM . ORCID: http://orcid.org/0000-0002-4487-8742 a r X i v : . [ h e p - ph ] M a y ontents A.1 SM ν - e and ν - ν interaction rates 29A.2 Comparison with previous literature on the SM: supplementary material 31A.3 N eff in the Standard Model including neutrino chemical potentials 36A.4 Detailed neutrinophilic scalar thermodynamics 37A.5 Evolution equations in the Maxwell-Boltzmann approximation 39A.6 Thermodynamic formulae 42A.7 Number and energy density transfer rates in the MB approximation 44– 1 – Introduction
The Planck collaboration [2] reports unprecedented precision measurements of the number ofeffective relativistic neutrino species, N eff . Within the framework of ΛCDM, Planck legacydata analyses yield: N eff = 2 . ± .
34 at 95% CL [3]. Additionally, upcoming CMB ex-periments like SPT-3G [4] and the Simons Observatory [5] will soon improve upon Planck’sprecision on N eff , and future experiments such as CMB-S4 [6], PICO [7], CORE [8] or CMB-HD [9] are expected to deliver a 1% precision determination of N eff . Likewise, the primordialhelium and deuterium abundances are now measured with 1% precision [10]; and in thefuture, the primordial deuterium abundance could be measured with 0.1% precision [11, 12].Precision determinations of the primordial element abundances and N eff simultaneouslyserve as confirmation of the thermal history in the Standard Model (SM) and a strongconstraint on many scenarios beyond, see e.g. [13–27]. However, precision measurements alsorepresent a challenge in finding accurate theoretical predictions for early Universe observables.In fact, obtaining an accurate prediction for N eff in the Standard Model is already a non-trivial task [28, 29] (see also [30–39]).The challenge in obtaining accurate predictions for N eff arises mainly as a result of thefact that neutrinos decouple from electrons and positrons in the early Universe at tempera-tures T dec ν ∼ N SMeff >
3. However, this is not the only relevant effect since finite temperature correctionsalter N eff at the level of ∆ N SMeff ∼ .
01 [40–42], and neutrinos start to oscillate prior to neu-trino decoupling [43, 44]. Thus, state-of-the-art treatments of neutrino decoupling accountfor non-thermal neutrino spectral distortions, finite temperature corrections and neutrinooscillations, leading to the Standard Model prediction of N SMeff = 3 .
045 [28, 29] . To obtainsuch an accurate number, state-of-the-art calculations resort to the density matrix formalismand solve a system of hundreds of stiff integro-differential equations for the neutrino distri-butions – something which is computationally challenging. Although this approach is veryaccurate, its complexity represents a drawback in studying generalized scenarios Beyond theStandard Model (BSM).In this work – building upon [1] – we propose a simplified, accurate and fast method tocalculate early Universe BSM thermodynamics. The method is based on the approximationthat any species can be described by thermal equilibrium distribution functions with evolvingtemperature and chemical potential. Within this approximation, simple ordinary differentialequations for the time evolution of the temperature and chemical potential of any givenspecies can be derived. These equations account for all relevant interactions, are easy tosolve, and track accurately the thermodynamics.This study extends the approach of [1] by allowing for non-negligible chemical potentials.This is important since, unlike in the SM, chemical potentials cannot be neglected in manyBSM theories. In addition, we include the next-to-leading order finite temperature correctionsto the electromagnetic plasma from [42], and spin-statistics and the electron mass in the SMreaction rates. Including these effects is required in order to obtain N eff in the StandardModel with an accuracy of 0 . Next to leading order finite temperature corrections are expected to shift this number by -0.001 [42]. – 2 –hen applying the method to neutrino decoupling in the Standard Model, and byaccounting for finite temperature corrections, and the electron mass and spin-statistics in the ν - ν and ν - e interaction rates, we find N SMeff = 3 . < m φ < λ ∼ O (10 − ) neutrinophilic scalar. The phenomenology and cosmology ofthis scenario was first highlighted in [45]. Recently, using the methods developed in thisstudy and by analyzing Planck legacy data [2, 3], strong constraints on the parameter spacewere derived in [46]. In this work, we compute all relevant thermodynamic quantities andcosmological parameters and compare them to the results obtained by solving the exactLiouville equation for the neutrino-scalar system. We find excellent agreement between thetwo computations and we therefore claim that the proposed method can be used to find fastand accurate thermodynamic quantities in BSM scenarios. We note that the method can beapplied to other frameworks not necessarily restricted to neutrino decoupling.Finally, to facilitate the implementation of the method to the interested reader, we pub-licly release a Mathematica and
Python code to solve for neutrino decoupling: NUDEC BSM.The typical execution time of the code in an average computer is O (10) s. Thanks to its sim-plicity, speed and accuracy we believe it can represent a useful tool to study BSM earlyUniverse thermodynamics. With NUDEC BSM , neutrino decoupling can be solved in the Stan-dard Model, in the presence of dark radiation, with MeV-scale species in thermal equilibriumduring neutrino decoupling, and in the presence of a light and weakly coupled neutrinophilicscalar. Given these examples, other scenarios should be easy to implement.
Structure of this work:
This work is organized as follows: In Section 2, we present the main method to solvefor the thermodynamic evolution of any species in the early Universe. We discuss the ap-proximations upon which the method is based and provide a first-principles derivation ofthe relevant evolution equations for the temperature and chemical potential describing thethermodynamics. In Section 3, we apply the method to neutrino decoupling in the StandardModel and compare with previous treatments. In Section 4, we solve for the thermal historyof the Universe in the presence of a very light and weakly coupled neutrinophilic scalar. Wepresent a detailed comparison between the solution using our method and that of solvingthe exact Liouville equation. In Section 5, we discuss the results, and argue on theoreticalgrounds the reason the proposed method is very accurate in many scenarios. We also discussthe limitations of our approach. In Section 6, we outline a recipe to model the early Universethermodynamics in generic extensions of the Standard Model. We conclude in Section 7.The practitioner is also referred to the Appendices A. There we outline various details,comparison between the proposed method against solutions to the Liouville equation, anduseful thermodynamic formulae. – 3 –
Fast and Precise Thermodynamics in the Early Universe
In this section we develop a simplified method to track the thermodynamic (in and out ofequilibrium) evolution of any species in the Standard Model and beyond. We first start byreviewing the Liouville equation that governs the time evolution of the distribution functionof a given species. Then, we assume that such distribution function is well-described bya thermal equilibrium distribution. Starting from the Liouville equation, we find ordinarydifferential equations for the time evolution of the temperature and a chemical potentialdescribing the thermodynamics. Finally, we provide useful formulae for the number andenergy density transfer rates for decays, annihilations and scatterings.
The distribution function f determines the thermodynamics of any given species in theearly Universe. In a fully homogeneous and isotropic Universe, the time evolution of thedistribution function f is governed by the Liouville equation [30, 47, 48]: ∂f∂t − Hp ∂f∂p = C [ f ] , (2.1)where p is the momentum, H = (8 πρ T / (3 m )) / is the Hubble rate, ρ T is the total energydensity of the Universe, m Pl = 1 . × GeV is the Planck mass, and C [ f ] is the collisionterm. In full generality, the collision term for a particle ψ is defined as [30]: C [ f ψ ] ≡ − E ψ (cid:88) X,Y (cid:90) (cid:89) i d Π X i (cid:89) j d Π Y j (2 π ) δ ( p ψ + p X − p Y ) × (2.2) |M| ψ + X → Y f ψ (cid:89) i f X i (cid:89) j (cid:2) ± f Y j (cid:3) − |M| Y → ψ + X (cid:89) j f Y j [1 ± f ψ ] (cid:89) i [1 ± f X i ] . The collision term accounts for any process ψ + X → Y where X ≡ (cid:80) i X i and Y ≡ (cid:80) j Y j represent the initial and final states accompanying ψ , and the indices i, j label particlenumber. Here, M ψ + X → Y is the probability amplitude for the process ψ + X → Y . The sumis taken over all possible initial ψ + X and final states Y , d Π X i ≡ π ) d p Xi E Xi , and the + signscorrespond to bosons, while the − signs correspond to fermions. The actual solution to the Liouville equation (2.1) strongly depends upon the processes thatthe given species is undergoing in the early Universe. It is well-known (see e.g. [48]), that wheninteractions between a particle ψ and the rest of the plasma are efficient, the distributionfunction of the ψ species is well described by equilibrium distribution functions. Namely,fermions and bosons follow the usual Fermi-Dirac (FD) and Bose-Einstein (BE) distributionfunctions: f FD ( E ) = 11 + e ( E − µ ) /T , f BE ( E ) = 1 − e ( E − µ ) /T , (2.3)respectively. Here, T is the temperature, µ is the chemical potential, and E = (cid:112) p + m isthe energy with m being the mass of the particle. Note that for bosons: µ < Upon integrating Equation (2.1) with the measures gd p/ (2 π ) and gEd p/ (2 π ) , for a par-ticle with g internal degrees of freedom, we find: dndt + 3 Hn = δnδt = (cid:90) g d p (2 π ) C [ f ] , (2.4) dρdt + 3 H ( ρ + p ) = δρδt = (cid:90) g E d p (2 π ) C [ f ] , (2.5)where n , ρ and p are the number, energy and pressure densities the species, respectively. δn/δt and δρ/δt represent the number and energy transfer rate between a given particle andthe rest of the plasma.By summing Equation (2.5) over all species in the Universe, one finds the usual conti-nuity equation: dρ T dt = − H ( ρ T + p T ) . (2.6)From a microphysics perspective, this continuity equation arises as a result of energy conser-vation in each process in the plasma; while from the fluid perspective, it simply results fromthe conservation of the total stress-energy tensor ∇ µ T µν = 0.By applying the chain rule to (2.4) and (2.5), we find: dTdt = − (cid:18) dndt ∂ρ∂µ − dρdt ∂n∂µ (cid:19) (cid:30) (cid:18) ∂n∂µ ∂ρ∂T − ∂n∂T ∂ρ∂µ (cid:19) , (2.7a) dµdt = − (cid:18) dρdt ∂n∂T − dndt ∂ρ∂T (cid:19) (cid:30) (cid:18) ∂n∂µ ∂ρ∂T − ∂n∂T ∂ρ∂µ (cid:19) . (2.7b)This set of equations can be considerably simplified if chemical potentials are neglected. Thistypically occurs as a result of some efficient interactions. In the Standard Model, for example,efficient e + e − ↔ γγ and e + e − ↔ γγγ interactions allow one to set µ e ( t ) = µ γ ( t ) = 0. If, dµ/dt = 0 then dndt = dρdt ∂n∂T / ∂ρ∂T and the previous equations (2.7) are simplified to: dTdt = dρdt (cid:30) ∂ρ∂T = (cid:20) − H ( ρ + p ) + δρδt (cid:21) (cid:30) ∂ρ∂T . (2.8)Therefore, Equation (2.8) can be used to track the thermodynamics of a species when chemicalpotentials are negligible (which occurs in many scenarios).– 5 –he set of Equations (2.7) can be rewritten in terms of the Hubble parameter and theenergy and number density transfer rates. Explicitly, they read: dTdt = 1 ∂n∂µ ∂ρ∂T − ∂n∂T ∂ρ∂µ (cid:20) − H (cid:18) ( p + ρ ) ∂n∂µ − n ∂ρ∂µ (cid:19) + ∂n∂µ δρδt − ∂ρ∂µ δnδt (cid:21) , (2.9a) dµdt = − ∂n∂µ ∂ρ∂T − ∂n∂T ∂ρ∂µ (cid:20) − H (cid:18) ( p + ρ ) ∂n∂T − n ∂ρ∂T (cid:19) + ∂n∂T δρδt − ∂ρ∂T δnδt (cid:21) . (2.9b)In general, n , ρ , p and their derivatives cannot be written analytically unless the given speciesfollows Maxwell-Boltzmann (MB) statistics and m = 0, or unless µ = m = 0 for FD, BEand MB statistics. Hence, in general, one requires numerical integration in order to obtainthermodynamic quantities. In Appendix A.6 we list all the relevant thermodynamic formulaewhile in Appendix A.5 we particularize the evolution equations for the case of MB statistics. In order to solve for the early Universe thermodynamics of a given system, the last ingredientneeded is to calculate the expressions for δn/δt and δρ/δt that enter Equations (2.9) and (2.8).This step is not generic and has to be carried out for each particular scenario. Here, in orderto facilitate the procedure, we write down relevant formulae for decays, annihilations andscattering processes that can be useful in many cases. We write down all the expressionsin the Maxwell-Boltzmann approximation. Namely, by considering f MB = e − ( E − µ ) /T . Thespin-statistics correction to the rates typically represents a (cid:46)
20% correction to the Maxwell-Boltzmann rate, but the advantage of the rates in the Maxwell-Boltzmann approximation isthat they can be written analytically for 1 ↔ ↔ Decay rates
The collision term for 1 ↔ a , characterized by T a and µ a , decays into two identical and massless statescharacterized by T (cid:48) and µ (cid:48) . In the Maxwell-Boltzmann approximation, the rates read asfollows (see Appendix A.7): δn a δt = Γ a g a m a π (cid:20) T (cid:48) e µ (cid:48) T (cid:48) K (cid:16) m a T (cid:48) (cid:17) − T a e µaTa K (cid:18) m a T a (cid:19)(cid:21) , (2.10) δρ a δt = Γ a g a m a π (cid:20) T (cid:48) e µ (cid:48) T (cid:48) K (cid:16) m a T (cid:48) (cid:17) − T a e µaTa K (cid:18) m a T a (cid:19)(cid:21) , (2.11)where Γ a is the decay rate at rest, g a is the number of internal degrees of freedom of the a particle, m a is its mass, and K j are modified Bessel functions of the j kind.– 6 – nnihilation rates For the process 1 + 2 ↔ T = T = T , µ = µ = µ , T = T = T (cid:48) , and µ = µ = µ (cid:48) , the number and energy density transfer rate for the 1 particle read (seeAppendix A.7): δn δt = g g π (cid:90) ∞ s min ds p √ s σ ( s ) (cid:20) T (cid:48) e µ (cid:48) T (cid:48) K (cid:18) √ sT (cid:48) (cid:19) − T e µT K (cid:18) √ sT (cid:19)(cid:21) , (2.12) δρ δt = g g π (cid:90) ∞ s min ds p (cid:0) s + m − m (cid:1) σ ( s ) (cid:20) T (cid:48) e µ (cid:48) T (cid:48) K (cid:18) √ sT (cid:48) (cid:19) − T e µT K (cid:18) √ sT (cid:19)(cid:21) , (2.13)where s is the Mandelstam variable, s min = min[( m + m ) , ( m + m ) ], p = [ s − ( m + m ) ] / [ s − ( m − m ) ] / / (2 √ s ), and σ ( s ) is the usual cross section for the process (namely,the cross section summed over final spin states and averaged over initial spins).These equations can be further simplified if all species are massless. Yielding: δn δt = g g π (cid:90) ∞ ds s / σ ( s ) (cid:20) T (cid:48) e µ (cid:48) T (cid:48) K (cid:18) √ sT (cid:48) (cid:19) − T e µT K (cid:18) √ sT (cid:19)(cid:21) , (2.14) δρ δt = g g π (cid:90) ∞ ds s σ ( s ) (cid:20) T (cid:48) e µ (cid:48) T (cid:48) K (cid:18) √ sT (cid:48) (cid:19) − T e µT K (cid:18) √ sT (cid:19)(cid:21) . (2.15)One must be careful with the spin factors. For example, for the energy transfer rate of asingle neutrino species and for the process ν α ¯ ν α ↔ e + e − , g = g = 1. On the other hand, ifone considers the energy transfer rate of an electron and the same process, then g = g = 2.Of course, the spin averaged/summed cross section takes into account this factor of 4 so that δρ e /δt = − δρ ν /δt .We exemplify the use of Equations (2.14) and (2.15) for typical annihilation cross sections: σ ( s ) = A/s −→ δn δt = A g g π T (cid:48) e µ (cid:48) T (cid:48) , δρ δt = A g g π T (cid:48) e µ (cid:48) T (cid:48) , (2.16) σ ( s ) = A sM (cid:63) −→ δn δt = AM (cid:63) g g π T (cid:48) e µ (cid:48) T (cid:48) , δρ δt = AM (cid:63) g g π T (cid:48) e µ (cid:48) T (cid:48) , (2.17)where we only show the 3 + 4 → δn /δt and δρ /δt for concreteness. Scattering rates
Scattering processes by definition yield δn/δt = 0 and if annihilation or decay interactionsare efficient the energy exchange from scatterings is typically subdominant. For example,in the Standard Model, electron-positron annihilations to neutrinos are a factor of 5 timesmore efficient transferring energy than electron-neutrino scatterings [1]. However, even ifscattering interactions can typically be neglected in the energy transfer rates, their effect isvery important since efficient scatterings distribute energy within the relevant species and leadto distribution functions that resemble equilibrium ones – and this is the key approximationthat leads to a simple evolution of the thermodynamics, see Section 2.2.Scattering rates are substantially more complicated to work out than decay or annihi-lation rates and simple expressions for them can only be found in very few cases. One suchcase is the scattering between two massless particles that interact via a four-point interac-tion. We perform the exact phase space integration and report the rates for such process inAppendix A.7 in the Maxwell-Boltzmann approximation. We also provide spin-statistics cor-rections for them. Finally, we point the reader to [52–54] where, in the context of WIMP darkmatter, scatterings of non-relativistic particles against massless ones have been investigatedin detail. – 7 –
Neutrino Decoupling in the Standard Model
In this section we solve for neutrino decoupling in the Standard Model. To this end, byusing Equation (2.8) we write down differential equations describing the time evolution ofthe temperature of the electromagnetic plasma T γ and the neutrino fluid T ν . This calculationis therefore fully analogous to the one performed by the author in [1], but here we update itby accounting for next-to-leading order finite temperature corrections [42], and spin-statisticsand the electron mass in the ν - ν and ν - e interaction rates. We solve for neutrino decoupling in the Standard Model by using the following well-justifiedapproximations:1.
Neutrinos follow perfect Fermi-Dirac distributions . Non-thermal corrections to an in-stantaneously decoupled neutrino distribution function are present in the StandardModel, but they encode less than 1% of the total energy density carried by the neu-trinos [28, 29]. Here we will show (see Section 3.5) that, in fact, thermal distributionfunctions can account very accurately for the energy density carried out by neutrinos.2.
Neglect neutrino oscillations . Neutrino oscillations are active at the time of neutrinodecoupling [28, 29, 43, 44]. However, we neglect them on the grounds that their impacton N eff in the Standard Model is small, ∆ N SMeff = 0 . T ν = T ν e = T ν µ = T ν τ .3. Neglect chemical potentials . The number of photons is not conserved in the earlyUniverse and the electron chemical potential is negligible given the very small baryon-to-photon ratio, µ e /T γ ∼ − . This therefore justifies setting µ e = µ γ = 0. In addition,in the Standard Model, neutrino chemical potentials can be neglected since interactions¯ νν ↔ e + e − ↔ γγγ are highly efficient in the relevant temperature range. Nonetheless,we have explicitly checked that allowing for non-negligible neutrino chemical potentialsdoes not render any significant change on any of the results presented in this Section,see Appendix A.3 for details. At the time of neutrino decoupling, electrons and positrons are tightly coupled to photons viaefficient annihilations and scatterings. Therefore we can model the electromagnetic sectorby the photon temperature T γ = T e . Since chemical potentials are negligible, by usingEquation (2.8) we can write the time evolution of the photon temperature in the SM as: dT γ dt = − Hρ γ + 3 H ( ρ e + p e ) + δρ νe δt + 2 δρ νµ δt∂ρ γ ∂T γ + ∂ρ e ∂T γ . (3.1)Here, δρ ν /δt accounts for the energy transfer between one neutrino species (including theantineutrino) and the rest of the electromagnetic plasma. δρ ν /δt enters this equation becauseenergy conservation ensures δρ e /δt = − δρ ν e /δt − ρ ν µ /δt .– 8 –e account for the next-to-leading order (NLO) Quantum Electrodynamics (QED)finite temperature correction to the electromagnetic energy and pressure densities from [42].Across the entire paper, however, when comparing with [28, 29] we shall use the leading order(LO) corrections from [40, 41] as used in [28, 29] to render a fair comparison. We denotethe finite temperature corrections to the electromagnetic pressure and energy densities as: P int and ρ int = − P int + T γ dP int dT γ , respectively. Where, in the notation of [42], dP int /dT γ ≡ (2 / T γ [ G ( m e T γ ) − m e T γ G ( m e T γ )] and d P int /dT γ ≡ T γ G ( m e T γ ). Expressions for P int and G , can be found in [42]. In the code we simply use a tabulated version of P int and its derivatives.By using Equation (2.8), we find the photon temperature evolution including finitetemperature corrections to be: dT γ dt = − Hρ γ + 3 H ( ρ e + p e ) + 3 HT γ dP int dT γ + δρ νe δt + 2 δρ νµ δt∂ρ γ ∂T γ + ∂ρ e ∂T γ + T γ d P int dT γ . (3.2)The neutrino temperature evolution can also be obtained from (2.8) and simply reads: dT ν α dt = − H T ν α + δρ ν α δt (cid:30) ∂ρ ν α ∂T ν α . (3.3)Note that the differential equations we have written for the neutrino temperature can beapplied to each neutrino flavor α or to the collective neutrino fluid. Applying them to theentire neutrino fluid is a good way of taking into account the fact that neutrino oscillationsbecome active prior to neutrino decoupling. In order to solve the relevant system of differential equations, we need to calculate theenergy transfer rates δρ ν α /δt . We take into account all relevant processes: e + e − ↔ ¯ ν α ν α , e ± ν α ↔ e ± ν α , e ± ¯ ν α ↔ e ± ¯ ν α , ν α ν β ↔ ν α ν β , ν α ¯ ν β ↔ ν α ¯ ν β , and ¯ ν α ν α ↔ ¯ ν β ν β . We consider thematrix elements reported in Table 2 of [30], and we follow the integration method developedin [32] in order to reduce the collision term to just two dimensions . We write the energytransfer rate as a function of T γ and T ν by making corrections to the analytical Maxwell-Boltzmann rates (see Appendix A.1 for details): δρ ν e δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = G F π (cid:2) (cid:0) g eL + g eR (cid:1) F ( T γ , T ν e ) + 2 F ( T ν µ , T ν e ) (cid:3) , (3.4a) δρ ν µ δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = G F π (cid:2) (cid:0) g µL + g µR (cid:1) F ( T γ , T ν µ ) − F ( T ν µ , T ν e ) (cid:3) , (3.4b)where, G F is Fermi’s constant, and as relevant for the energies of interest ( E (cid:28) M Z ) [10, 57] : g eL = 0 . , g eR = 0 . , g µL = − . , g µR = 0 . , (3.5)and where F ( T , T ) = 32 f FD a ( T − T ) + 56 f FD s T T ( T − T ) . (3.6)with f FD a = 0 . f FD s = 0 . m e . The reader is referred to the appendices of [55] and [56] for recent and detailed descriptions of the method.The exact phase space reduced matrix elements are available from the author upon request. In the first version of the paper we used the tree level values for g L and g R with s W = 1 − M W /M Z = 0 . ν − e rates in this version are 1 .
8% larger. – 9 – eutrino Decoupling BSM UCL 27-03-19Miguel Escudero (KCL) 1 e ± ⌫ W/Z
20 10 5 3 2 1 0.8 0.5 0.1 0.07 0.01 T γ (MeV)1.01.11.21.31.4 T γ / T ν N eu t r i no D e c oup li ng p ↔ n F r ee z e - ou t H e li u m S y n t he s i s e + e − → γ γ Fast, N SMeff = 3.045Instantaneous, N eff = 3Evolution in the Standard Model Figure 1 . Neutrino temperature evolution in the Standard Model. Key cosmological events arehighlighted: neutrino decoupling, proton-to-neutron freeze-out, and helium synthesis. The red line isthe result obtained in this work that leads to N eff = 3 .
045 and the black dashed line corresponds tothe case of instantaneous decoupling, leading to N eff = 3. The upper cartoons represent snapshotsof the relative number density of neutrinos (red), electrons (purple) and photons (blue). A data filecontaining the thermodynamic evolution in the Standard Model can be found on ArXiv. Here we present the results of solving the system of Equations (3.2) and (3.3). We solve themstarting from a sufficiently high temperature so that neutrino-electron interactions are highlyefficient. In particular, we start the integration from T γ = T ν = 10 MeV, corresponding to t = 1 / (2 H ) | T =10 MeV . We track the system until T γ ∼ .
01 MeV, at which point neutrinoshave fully decoupled from the plasma and electron-positron annihilation has taken place.We use 10 − as the relative and absolute accuracies for the integrator which yields a highnumerical accuracy. Given these settings, the typical CPU time used to integrate theseequations in Mathematica is ∼
20 s and in
Python (cid:46)
10 s. We have explicitly checked thatthe continuity equation is fulfilled at each integration step to a level of 10 − or better.In Figure 1 we show the neutrino temperature evolution as a function of T γ . We highlightsome key cosmological events and also compare with instantaneous neutrino decoupling.The central parameter to this study is the number of effective relativistic neutrinospecies as relevant for CMB observations: N eff . It is explicitly defined as: N eff ≡ (cid:18) (cid:19) / (cid:18) ρ rad − ρ γ ρ γ (cid:19) = 3 (cid:18) (cid:19) / (cid:18) T ν T γ (cid:19) , (3.7)where in the last step we have assumed that the neutrinos follow perfect Fermi-Dirac distri-bution functions with negligible chemical potentials.In Table 1, we report the resulting values of N eff and T γ /T ν for T γ (cid:28) m e in theSM under various approximations and considering both a common neutrino temperature( T ν = T ν e = T ν µ = T ν τ ) and by evolving separately T ν e and T ν µ, τ = T ν µ = T ν τ .– 10 –eutrino Decoupling in the SM T ν e = T ν µ, τ T ν e (cid:54) = T ν µ, τ Scenario T γ /T ν N eff T γ /T ν e T γ /T ν µ N eff Instantaneous decoupling 1.4010 3.000 1.4010 1.4010 3.000Instantaneous decoupling + LO-QED 1.3997 3.011 1.3997 1.3997 3.011Instantaneous decoupling + NLO-QED 1.3998 3.010 1.3998 1.3998 3.010MB collision term 1.3961 3.043 1.3946 1.3970 3.042MB collision term + NLO-QED 1.3950 3.052 1.3935 1.3959 3.051FD collision term 1.3965 3.039 1.3951 1.3973 3.038FD collision term + NLO-QED 1.3954 3.049 1.3941 1.3962 3.048FD+ m e collision term 1.3969 3.036 1.3957 1.3976 3.035FD+ m e collision term + LO-QED 1.39568 3.046 1.3945 1.3964 3.045 FD+ m e collision term + NLO-QED Table 1 . N eff and T γ /T ν as relevant for CMB observations in the SM by taking different approxima-tions and neglecting neutrino chemical potentials. The last row shows the case in which we accountfor both quantum statistics and the electron mass in the relevant collision terms, for which we find N SMeff = 3 . − . One of our main results is that – by considering spin-statistics and m e in the e - ν and ν - ν energy transfer rates – we find N SMeff = 3 . N eff in the SM [28, 29] that account for non-thermal correctionsto the neutrino distribution functions and neutrino oscillations. The reader is referred toAppendix A.2 for a comparison in terms of the neutrino distribution functions at T (cid:28) m e .Thus, we have found that although describing the neutrino populations by a Fermi-Dirac distribution functions is just an approximation to the actual scenario, it suffices forthe purpose of computing N SMeff with a remarkable accuracy.
We compare our results for early and late Universe observables from our calculation of neu-trino decoupling with state-of-the-art calculations in the Standard Model [28, 29, 58, 59].These studies account for non-thermal neutrino distribution functions, finite temperaturecorrections, and neutrino oscillations in the primordial Universe. We compare our results interms of N eff , the energy density of degenerate non-relativistic neutrinos (Ω ν h ), the effec-tive number of species contributing to entropy density ( g (cid:63)s ), and the primordial abundancesof helium ( Y P ) and deuterium (D/H P ). To obtain the relative differences in terms of Y P and D/H P we have modified the BBN code PArthENoPE [58, 59]. We refer the reader toAppendix A.2 for details.In Table 2 we outline our main results and comparison with previous state-of-the-artliterature. We find an agreement of better than 0.1% for any cosmological parameter. Theaccuracy of our approach in the Standard Model is well within the sensitivity of future CMBexperiments to N eff [4–9] and of future measurements of the light element abundances [11, 12].Finally, in addition to accuracy, we stress that the two other key features of our approachare simplicity and speed. One needs to solve for a handful of ordinary differential equationsand the typical execution time of NUDEC BSM is ∼
20 s in
Mathematica and ∼
10 s in
Python .Thus, we believe this approach has all necessary features to be used to model early UniverseBSM thermodynamics. This is the subject of study of the next sections.– 11 –eutrino Decoupling in the Standard Model: Key Parameters and ObservablesParameter N eff Y P D/H | P g (cid:63) s (cid:80) m ν / Ω ν h This work 3.045 - - 3.931 93.05 eVDifference w.r.t. instantaneous ν -dec 1.5 % 0.1 % 0.4 % 0.6 % 1.2 %Difference w.r.t. [28, 29, 58, 59] 0.03 % 0.008 % 0.08 % 0.05 % 0.09 %Current precision [3, 10] 5-6 % 1.2 % 1.1 % - -Future precision [5, 6, 11, 12] 1-2 % < Table 2 . Standard Model results as obtained in this work (NUDEC BSM) for relevant cosmologicalparameters by solving for the time evolution of T ν and T γ neglecting chemical potentials. The values of Y P and D/H P have been calculated using a modified version of PArthENoPE [58, 59]. We compare ourresults in terms of N eff , g (cid:63) s and (cid:80) m ν / Ω ν h against state-of-the-art calculations in the SM [28, 29].We compare the results in terms of Y P and D/H | P against the default mode of PArthENoPE [58, 59]that includes non-instantaneous neutrino decoupling as obtained in [29]. Current precision on N eff isfrom Planck [3] and on Y P and D/H | P is from the PDG [10]. Expected future precision on N eff is fromthe Simons Observatory [5] and CMB-S4 [6]. Expected precision on Y P and D/H | P is from [11, 12]. In this section we study the thermal history of the Universe in the presence of a very light(1 eV < m φ < λ < − ) neutrinophilic scalar: φ . This isprototypically the case of majorons [60–63] where the very small coupling strengths areassociated to the small neutrino masses, and where sub-MeV φ masses are consistent withthe explicit breaking of global symmetries by gravity. We will assume that the neutrinophilicscalar posses the same couplings to all neutrino flavors. Furthermore, we note that even if thescalar does not couple with the same strength to each neutrino flavor, neutrino oscillationsin the early Universe [43, 44] will render T ν e (cid:39) T ν µ (cid:39) T ν τ regardless.We begin by explicitly defining the model we consider. We follow by posing and solvingthe relevant system of differential equations for the temperature and chemical potentials forthe neutrinos and the φ scalar. We solve this system of equations and briefly comment onthe phenomenology in the region of parameter space where neutrino-scalar interactions arehighly efficient in the early Universe. We finally write down the Liouville equation for theneutrino-scalar system. We solve it and compare with the results as output by the methodproposed in this work (NUDEC BSM). We find an excellent agreement between the twoapproaches. The interaction Lagrangian that describes this model is L int = λ φ (cid:88) i ¯ ν i γ ν i , (4.1)where λ is a coupling constant, ν i corresponds to a neutrino mass eigenstate, we have assumedthat neutrinos are Majorana particles, and we will restrict ourselves to 1 eV < m φ < φ particles and we shall focus on the regime in which the neutrino- φ interaction strength isvery small λ < − . In this regime, 2 ↔ − − T ν / m φ − − − − h Γ ¯ νν → φ i / H Γ eff ¯ ⌫⌫ Figure 2 . Left:
Ratio of the thermally averaged neutrino inverse decay rate (¯ νν → φ ) to the Hubbleparameter as a function of temperature for different neutrino- φ interaction strengths, Γ eff (4.4). Right: φ → ¯ νν decay diagram. Therefore, in this scenario, the only cosmologically relevant processes are decays of φ into neutrinos and neutrino inverse decays – as depicted in Figure 2. The decay width atrest of φ reads: Γ φ = λ π m φ (cid:88) i (cid:115) − m ν i m φ (cid:39) λ π m φ , (4.2)where in the last step we have neglected neutrino masses.Given the decay width at rest, the number and energy density transfer rates as a result of φ ↔ ¯ νν interactions, in the Maxwell-Boltzmann approximation, are given by Equations (2.10)and (2.11) respectively so that: δn φ δt = Γ φ m φ π (cid:20) T ν e µνTν K (cid:18) m φ T ν (cid:19) − T φ e µφTφ K (cid:18) m φ T φ (cid:19)(cid:21) , (4.3a) δρ φ δt = Γ φ m φ π (cid:20) T ν e µνTν K (cid:18) m φ T ν (cid:19) − T φ e µφTφ K (cid:18) m φ T φ (cid:19)(cid:21) . (4.3b)Since the relevant process is 1 ↔
2, the collective neutrino fluid interaction rates fulfil: δρ ν /δt = − δρ φ /δt and δn ν /δt = − δn φ /δt .In order to efficiently explore the parameter space and to understand the extent towhich departures from thermal equilibrium occur, it is convenient to introduce the effectiveinteraction strength rate:Γ eff ≡ (cid:18) λ × − (cid:19) (cid:18) keV m φ (cid:19) (cid:39) (cid:104) Γ ¯ νν → φ (cid:105) H = 1 H ( T ν = m φ /
3) 1 n ν δn ν δt (cid:12)(cid:12)(cid:12)(cid:12) ¯ νν → φ . (4.4)Figure 2 shows the temperature dependence of the thermally averaged interaction rate ¯ νν → φ , and we can clearly appreciate that if Γ eff > φ and ν species. – 13 – .2 Evolution equations In order to describe the early Universe thermodynamics of the scenario described abovewe follow the method outlined in Section 2. We can simply use two sets of Equations (2.9)tracking the scalar and neutrino temperature and chemical potentials ( T φ , µ φ , T ν , µ ν ) sourcedwith the number and energy density transfer rates induced by φ ↔ ¯ νν interactions (4.3).For concreteness, we shall assume that φ ↔ ¯ νν interactions become relevant after electron-positron annihilation (although including the effect of e + e − annihilations is straightforward).Explicitly, the evolution equations describing the thermodynamics read: dT φ dt = 1 ∂n φ ∂µ φ ∂ρ φ ∂T φ − ∂n φ ∂T φ ∂ρ φ ∂µ φ (cid:20) − H (cid:18) ( p φ + ρ φ ) ∂n φ ∂µ φ − n φ ∂ρ φ ∂µ φ (cid:19) + ∂n φ ∂µ φ δρ φ δt − ∂ρ φ ∂µ φ δn φ δt (cid:21) , (4.5a) dµ φ dt = − ∂n φ ∂µ φ ∂ρ φ ∂T φ − ∂n φ ∂T φ ∂ρ φ ∂µ φ (cid:20) − H (cid:18) ( p φ + ρ φ ) ∂n φ ∂T φ − n φ ∂ρ φ ∂T φ (cid:19) + ∂n φ ∂T φ δρ φ δt − ∂ρ φ ∂T φ δn φ δt (cid:21) , (4.5b) dT ν dt = 1 ∂n ν ∂µ ν ∂ρ ν ∂T ν − ∂n ν ∂T ν ∂ρ ν ∂µ ν (cid:20) − H (cid:18) ( p ν + ρ ν ) ∂n ν ∂µ ν − n ν ∂ρ ν ∂µ ν (cid:19) − ∂n ν ∂µ ν δρ φ δt + 2 ∂ρ ν ∂µ ν δn φ δt (cid:21) , (4.5c) dµ ν dt = − ∂n ν ∂µ ν ∂ρ ν ∂T ν − ∂n ν ∂T ν ∂ρ ν ∂µ ν (cid:20) − H (cid:18) ( p ν + ρ ν ) ∂n ν ∂T ν − n ν ∂ρ ν ∂T ν (cid:19) − ∂n ν ∂T ν δρ φ δt + 2 ∂ρ ν ∂T ν δn φ δt (cid:21) , (4.5d) dT γ dt = − H T γ , (4.5e)where δn φ /δt and δρ φ /δt are given by (4.3) and the thermodynamic quantities of the neutrinoscorrespond to all neutrinos and antineutrinos, namely, to a Fermi-Dirac fluid with g ν = 6. Initial Conditions:
We solve this set of equations with initial conditions: T ν /m φ = 100 , T γ /T ν = 1 . , µ ν /T ν = − − , T φ /T ν = 10 − , µ φ /T ν = − − , (4.6)and t = 1 / (2 H ). The initial condition for T γ /T ν corresponds to the temperature ratio asobtained after e + e − annihilation in the Standard Model. The initial conditions for T φ and µ φ ensure that the integration starts with a tiny abundance of φ particles, ρ φ /ρ ν | t < − .In practice, when solving the system of Equations (4.5) we split the integration in two.We integrate from T ν /m φ = 100 until T ν /m φ = 20 setting m φ = 0 in the thermodynamicformulae since that allows for a great simplification and it is clearly a good approximation inthat regime. Then, we ran from T ν /m φ = 20 using thermodynamic formulae with m φ (cid:54) = 0. Parameter Space:
We solve the set of Equations (4.5) for the range 10 − < Γ eff < . Since we consider1 eV < m φ < − < λ < − . Notethat the results presented here will only depend upon Γ eff and apply to the entire range1 eV < m φ < φ particles becomecosmologically relevant in a radiation dominated Universe and the evolution is only dependentupon ratios of temperature/chemical potential to m φ . Thus, for one single value of m φ , theresults can be mapped into the entire range of φ masses for a given Γ eff . For m φ < .3 Thermodynamics in the strongly interacting regime If Γ eff (cid:29)
1, then φ ↔ ¯ νν interactions are highly efficient in the early Universe and thermalequilibrium is reached between the ν and φ populations. In this regime, one can find thethermodynamic evolution of the joint ν - φ system independently of Γ eff and without the actualneed to solve for Equations (4.5) .When Γ eff (cid:29)
1, the φ scalars thermalize with neutrinos while relativistic, T (cid:29) m φ .Since at such temperatures all particles can be regarded as massless, entropy is conservedand a population of φ particles is produced at the expense of neutrinos. In this regime, energyand number density conservation can be used to find an expression for the temperature andchemical potential of the joint ν - φ system after equilibration occurs. By denoting T ν as thetemperature before equilibration, and T eq and µ eq the temperature and chemical potentialof the neutrinos after equilibration, we can express the previous conditions as: ρ ν ( T ν ,
0) = ρ ν ( T eq , µ eq ) + ρ φ ( T eq , µ eq ) , (4.7) n ν ( T ν ,
0) = n ν ( T eq , µ eq ) + 2 n φ ( T eq , µ eq ) , (4.8)where µ φ eq = 2 µ ν eq as required from the fact that φ ↔ ¯ νν processes are highly efficient. Wecan easily solve these equations to find: T eq = 1 . T ν , µ eq = − . T ν . (4.9)This result bounds the maximum energy density of the φ species to be: ρ φ ( T eq , µ eq ) ρ ν ( T eq , µ eq ) + ρ φ ( T eq , µ eq ) (cid:39) .
09 and equivalently ρ φ T γ < . . (4.10)Therefore, the energy density in φ ’s represents less than 10% of the total ν - φ system.Once thermal equilibrium is reached the thermodynamic evolution of the system is fixed.This is a result of the fact that in thermal equilibrium entropy is conserved in the joint ν - φ system, and we know that for every φ that decays a pair of ¯ νν are produced. This allows usto write down two conservation equations for the entropy and number densities: (cid:2) s ν ( T (cid:48) , µ (cid:48) ) + s φ ( T (cid:48) , µ (cid:48) ) (cid:3) a (cid:48) = [ s ν ( T, µ ) + s φ ( T, µ )] a , (4.11a) (cid:2) n ν ( T (cid:48) , µ (cid:48) ) + 2 n φ ( T (cid:48) , µ (cid:48) ) (cid:3) a (cid:48) = [ n ν ( T, µ ) + 2 n φ ( T, µ )] a . (4.11b)And therefore, this set of equations – when solved simultaneously – provide the evolutionof T (cid:48) and µ (cid:48) as a function of the scale factor a (cid:48) . In this case: T (cid:48) γ a (cid:48) = T γ a , and we can findthe evolution as a function of T γ . We can solve the system starting from the equilibriumconditions when T (cid:29) m φ until T (cid:28) m φ , so that the φ population has disappeared from theplasma. By following this procedure we find that after the φ particles decay, the temperatureand neutrino chemical potentials are: T γ /T ν = 1 . , T ν /µ ν = − . , ( T (cid:28) m φ , Γ eff (cid:29) , (4.12)which implies ∆ N eff = N eff − N SMeff = 3 . − .
045 = 0 . (cid:39) . . (4.13)This is in excellent agreement with the results of [45]. Note that so long as Γ eff (cid:29)
1, theseexpressions are independent of the actual coupling strength. This is analogous to what happens in the SM with electron-positron annihilation. The net heating of thephoton bath does not depend upon the interaction strength of electromagnetic interactions but only upon thenumber of internal degrees of freedom of the electrons and photons and their fermionic and bosonic nature. – 15 – ! ¯ ⌫ ⌫ ¯ ⌫ ⌫ ! / / / /
10 1 / T ν / m φ ρ / T γ νφ Γ eff = 1 ↑ ∆ N eff = 0.11 Figure 3 . Energy density evolution of a neutrinophilic scalar φ (blue) coupled to neutrinos (red) ina region of parameter space in which φ -decays to neutrinos and neutrino inverse decays are the onlyrelevant processes. The interaction strength has been chosen to be Γ eff = 1, which corresponds toa neutrino- φ coupling strength of λ = 4 × − (cid:112) m φ / keV. The sketch illustrates an approximatesnapshot in comoving space of the joint thermodynamic system. At very high temperatures, the¯ νν → φ process is thermally suppressed as a result of the high boost of neutrinos which suppressesthe phase space for producing very light φ particles. When T ν ∼ m φ , ¯ νν → φ processes start to beefficient. As soon as T ν ∼ m φ / φ particles feel they are non-relativistic and decay to neutrinos.However, by then, the φ population is non-relativistic and the decay products from φ → ¯ νν are moreenergetic than the typical neutrino in the plasma, thereby yielding ∆ N eff = 0 .
11 as relevant for CMBobservations for 1 eV (cid:46) m φ (cid:46) MeV. The reader is referred to [45] and [46] for the cosmologicalimplications of this scenario.
In this section we present the solution to the thermodynamic evolution of the early Universein the presence of a light and weakly coupled neutrinophilic scalar.First, in Figure 3 we show the evolution of the energy density in neutrinos and φ speciestogether with an sketch of the joint thermodynamic system. We can appreciate that themain effect of efficient φ ↔ ¯ νν interactions (Γ eff >
1) is to render a siezable population of φ particles at T (cid:38) m φ . As a result of the fact that these particles have a mass, when thetemperature of the Universe is T (cid:46) m φ , they do not redshift as pure radiation and upondecay they render a more energetic neutrino population (relative to the photons) than theone at T (cid:29) m φ .In Figure 4 we show the resulting value of ∆ N eff as a function of Γ eff . We can clearlyappreciate that when Γ eff (cid:29)
1, ∆ N eff asymptotically reaches the value obtained in (4.13).We see that wide regions of parameter space in this scenario are within the reach of nextgeneration of CMB experiments [5, 6]. We note that for all the cases presented in this study,we have explicitly checked that the continuity equation (2.6) is fulfilled at each integrationstep with a relative accuracy of better than 10 − . Therefore the resulting values of ∆ N eff are numerically accurate to a level of ∼ . − − − − Γ eff ∆ N e ff Simons Observatory (1 σ )CMB Stage-IV (1 σ ) Figure 4 . ∆ N eff as a function of Γ eff as relevant for CMB observations and as obtained by solvingEquations (4.5). The 1 σ sensitivities for the Simons Observatory [5] and CMB-S4 [6] are shownfor illustration. We appreciate that when Γ eff (cid:29) N eff matches the perfect thermal equilibriumprediction in (4.13). − − T ν / m φ × ρ φ / T γ Γ eff − − T ν / m φ ( ρ φ + ρ ν ) / T γ Γ eff Thermal Equilibrium
Figure 5 . Left panel:
Evolution of ρ φ as a function of temperature for various choices of Γ eff . Notethat ρ φ never exceeds the thermal equilibrium prediction of (4.10). Right panel:
Evolution of ρ φ + ρ ν as a function of temperature. The dashed line indicates the evolution as dictated by Equation (4.11). In the left panel of Figure 5, we show the evolution of the energy density of φ speciesfor representative values of the interaction strength Γ eff . We can clearly appreciate thatfor Γ eff (cid:38) φ scalar at T (cid:46) m φ is the same – as dictated by entropyconservation. Similarly, in the right panel of Figure 5, we show the evolution of ρ φ + ρ ν and appreciate that, as expected, for Γ eff (cid:29) T γ /T ν and T ν /µ ν for T (cid:28) m φ as a function of Γ eff .We clearly see that, for Γ eff >
1, the values of T γ /T ν and T ν /m ν asymptotically reach thevalues of (4.12). – 17 – igure 6 . Values of T ν and µ ν after the φ population has completely decayed away, i.e. T ν (cid:28) m φ .We can appreciate that for Γ eff (cid:38)
10 both T ν and µ ν match those derived assuming perfect thermalequilibrium (dashed lines), as given by Equation (4.12). The aim of this section is to solve for the exact background thermodynamics of the scenariooutlined in Section 4.1. We do this by means of solving the Liouville equation for theneutrino and φ distribution functions. As discussed in the introduction, solving the Liouvilleequation requires one to solve a system of hundreds of stiff integro-differential equationsthat is computationally challenging. We therefore only find solutions for some representativevalues of Γ eff . We finally compare the results of solving the exact Liouville equation with theapproach presented here and find they are in excellent agreement. The Liouville equation
In the region of parameter space in which the only relevant interactions are φ ↔ ¯ νν , andafter the neutrinos have decoupled ( T γ < φ distribution functionsare characterized by the following system of Boltzmann equations [50, 64]: df φ dt − H p φ ∂f φ ∂p φ = C φ ↔ ¯ ννφ [ f φ ] = − m φ Γ φ E φ p φ (cid:90) Eφ + pφ Eφ − pφ dE ν F dec ( E φ , E ν , E φ − E ν ) , (4.14a) df ν dt − H p ν ∂f ν ∂p ν = C φ ↔ ¯ ννν [ f ν ] = m φ Γ φ E ν p ν (cid:90) ∞| ( m φ / p ν ) − p ν | dp φ p φ E φ F dec ( E φ , E ν , E φ − E ν ) , (4.14b)where f ν corresponds to the distribution function of one single neutrino species, and F dec is: F dec ( E φ , E ν , E ν (cid:48) ) = f φ ( E φ ) [1 − f ν ( E ν )] [1 − f ν ( E ν (cid:48) )] − f ν ( E ν ) f ν ( E ν (cid:48) ) [1 + f φ ( E φ )] . (4.15) Numerics to solve for f ν and f φ . We solve the system of integro-differential equa-tions (4.14) by binning the neutrino and φ distribution functions in comoving momentumin the range y ≡ p a/m φ = [0 . ,
20] in 200 bins each. This represents a high accuracysetting [65]. The system thus consists of 400 stiff integro-differential equations that we solveusing backward differentiation formulas from vode in Python . We use the default settingsfor the absolute and relative accuracies of the integrator.– 18 –
PU time usage.
The typical time needed to solve the exact Liouville equation (4.5)is roughly t CPU / day = 0 . . . . This is in sharp contrast with t CPU ∼ C or Fortran implementation of the solver could potentially yield an order ofmagnitude improvement in the speed but still would be orders of magnitude far from 1 min.
Initial Conditions and Integration Range.
We use as initial conditions for f ν and f φ Fermi-Dirac and Bose-Einstein formulas respectively with temperature and chemical poten-tials as in (4.6). We run the integrator until the largest time between T ν = m φ /
15 and t = 5 / Γ φ ( T ν ∼ m φ / (cid:112) Γ eff / . φ population has completelydecayed away by then. Namely, ρ φ /ρ ν < − . Comparison
Here we compare the results of solving the Liuoville equation (4.14) and the system of equa-tions (4.5) that describe the same system: a light neutrinophilic scalar where φ ↔ ¯ νν in-teractions are cosmologically relevant. We shall denote the solution to the exact Liouvilleequation as “Full” while we denote the solution to the system of differential equations (4.5)as “Fast”. We explicitly compare the solution of the two approaches in terms of the timeevolution of ρ φ + ρ ν and ∆ N eff . The reader is referred to Appendix A.4 for a comparison interms of ρ φ , ρ ν , n φ , n ν , f ν and f φ where very good agreement is also found.We focus the comparison for interaction strengths in the range 10 − < Γ eff <
20. ForΓ eff >
20 solving the Liouville equation is too time consuming ( t CPU ∼
25 (Γ eff / . days),while for Γ eff < − the cosmological consequences of the very light neutrinophilic scalar weconsider are negligible (as can be seen from Figure 4).In the right panel of Figure 7 we show the total energy density of the system, ρ φ + ρ ν .We appreciate an excellent agreement between the two approaches. For 10 − < Γ eff <
20 therelative difference between the two approaches in terms of ρ φ + ρ ν at any given temperatureis always better than 0 . N eff = N eff − .
045 in the left panel of Figure 7.Very similar results are found between the two approaches and only when Γ eff (cid:38) N eff . This difference is well below anyexpected future sensitivity from CMB observations. In fact, we believe that in the regimeΓ eff (cid:38) N eff is that given by the Fast approach and not the Full one.The reason is twofold: i) the Fast approach does converge to the value of ∆ N eff predictedby assuming pure thermal equilibrium, that must hold for Γ eff (cid:29) ii) the largerΓ eff is, the more numerically unstable the Liouville equation becomes. On the one hand,as Γ eff gets larger, f ν and f φ are more accurately described by perfect thermal equilibriumdistribution functions. On the other hand, if f ν and f φ resemble equilibrium distributionfunctions, Equation (4.15) tends to 0. Thus, we believe that for Γ eff (cid:38)
1, small numericalinstabilities prevent perfect thermalization of f ν and f φ leading to a slightly underestimatedvalue of ∆ N eff by the Full solution.Thus, we have demonstrated that the set of differential equations (2.9) describe verywell the thermodynamic evolution of the scalar-neutrino system considered in this work.Importantly, the total energy density of the system ( ρ φ + ρ ν ) agrees between the Full andFast solutions to better than 0.4% at any given temperature. In addition, the asymptoticvalues of ∆ N eff agree within 0.01 or better. Therefore, the Fast approach developed in thisstudy describes the thermodynamics of the system very accurately and can be used to studythe reach of ultrasensitive CMB experiments [6–9] to this particular scenario.– 19 – − − − − Γ eff ∆ N e ff Simons Observatory (1 σ )CMB Stage-IV (1 σ ) N eff -Fast N eff -Full ( t CPU / min ≃ t CPU / day ≃ Γ ) − T ν / m φ − − | − ( ρ φ + ρ ν ) | F AS T ( ρ φ + ρ ν ) | F U LL | ( % ) Γ eff Figure 7 . Left panel: ∆ N eff as a function of Γ eff . In black (Fast) the results from solving (4.5)and in red (Full) those by solving (4.14). The agreement is very good and the small differences arewell below any expected CMB experiment sensitivity. Right panel:
Relative difference between theFast and Full solution for ρ ν + ρ φ . The agreement between the two is better than 0.4% in all casesconsidered in this study. In this work we have shown that thermal equilibrium distributions can accurately track thethermal history of the Universe in various scenarios. We have explicitly demonstrated thisstatement in two specific cases: i) neutrino decoupling in the Standard Model (see Section 3)and ii) a SM extension featuring a very light and weakly coupled neutrinophilic scalar (seeSection 4). This poses the question of why this occurs given that out-of-equilibrium processesare relevant in the two cases. We believe that the reason is twofold:1. The evolution equations that dictate the temperature and chemical potential evolu-tion (2.9) describing each species arise from conservation of energy and number density– Equations (2.4) and (2.5) respectively. Given that they result from conservationequations for ρ and n , the main thermodynamic quantities should be well described.2. Departures from thermal equilibrium are not too severe in either of the two scenarios.(a) Neutrino Decoupling in the Standard Model . The out-of-equilibrium injection ofneutrinos from e + e − annihilations mainly occurs at T ∼ m e . The interactionrate for this process is only mildly sub-Hubble: Γ /H ∼ ( m e /T dec ν ) ∼ .
02, andthe neutrinos produced from e + e − annihilations have an energy E ν ∼ T γ + m e which is similar to the mean neutrino energy (cid:104) E ν (cid:105) ∼ T ν . Hence, a Fermi-Dirac distribution function for the neutrinos describes well the thermodynamics– provided that the temperature evolution accounts for the relevant interactions.(b) Light and Weakly Coupled Neutrinophilic Scalar . The interaction rate is a freeparameter in this scenario and we have found an excellent description of thethermodynamics for all situations in which decays and inverse decays fulfilledΓ eff ∼ Γ /H > − . In this case, Γ eff also controls the typical energy of theneutrinos injected from φ → ¯ νν decays. For Γ eff (cid:38) φ ↔ ¯ νν interactions have a typical energy E ν ∼ T ν and hence a thermal equilibrium distribution function characterizes very wellthe system. For Γ eff (cid:28) E φ → ¯ ννν ∼ T ν (cid:112) . / Γ eff , and therefore for Γ eff (cid:38) − the neutrinos injected have an energy sufficiently similar to (cid:104) E ν (cid:105) ∼ T ν andperfect thermal distribution functions can characterize the system accurately.– 20 –hus, given these examples, we believe that the thermodynamics of any system where de-partures from thermal equilibrium are moderate can be accurately described by the methodpresented in Section 2. Furthermore, the larger the interaction rates are, the better thesystem will be described by this approach. For very small interaction rates (Γ /H (cid:46) − )a given scenario should still be approximately described by our system of equations (2.9)but the accuracy of the approach could be substantially reduced. Of course, in this regime,solving the exact Liouville equation is not as challenging given the smallness of the interac-tion rates. This occurs in relevant cases such as sterile neutrino dark matter [66, 67], seee.g. [68], or freeze-in dark matter [69], see e.g. [70]. We finish by mentioning that, in fact,very out-of-equilibrum processes are particularly simple to solve and can be worked out atthe level of the energy and number density of particles, see for example Chapter 5.3 of [48]. Here we provide a four steps recipe to model the early Universe thermodynamics in genericBSM scenarios. The steps are the following:1.
Identify what are the relevant species and group them in sectors.
If interactions within a sector are strong, then one can write down evolution equationsfor the joint system rather than for each species. For example, this is the case in theStandard Model neutrino decoupling for electrons, positrons and photons.2.
Determine whether chemical potentials can be neglected or not.
Chemical potentials can be neglected if processes that do not conserve particle numberare active and there is no primordial asymmetry between particles and antiparticles.This leads to a considerable simplification of the evolution equations.3.
Include all relevant interactions among sectors.
Calculate the relevant energy and number density transfer rates from decay, annihila-tion and scattering processes. The expressions for decay and annihilation processes inthe Maxwell-Boltzmann approximation can be found in Section 2.4.4.
Write down evolution equations for T and µ for each sector within the scenario. Use Equations (2.9) if chemical potentials cannot be neglected and Equation (2.8)otherwise. In order to implement numerically the recipe, one can use as a startingpoint NUDEC BSM where a suite of BSM scenarios have already been coded up.
An example: Massless Neutrinophilic Scalar Thermalization during BBN
We illustrate the use of this recipe by considering the following scenario: a massless neu-trinophilic scalar that thermalizes during BBN. For concreteness, the scenario is describedby the same interaction Lagrangian as in Equation (4.1): L int = λ/ φ (cid:80) i ¯ ν i γ ν i . Given thisinteraction Lagrangian and since m φ = 0, kinematically, the only possible processes in which φ can participate are ¯ νν ↔ φφ , νφ ↔ νφ and ¯ νφ ↔ ¯ νφ .1. Identify what are the relevant species and group them in sectors.
There are three different sectors in this example: 1) e + , e − , γ , 2) ν and 3) φ . This is, if interactions are present, when the different species interact efficiently enough (Γ /H (cid:38) − )and when particles are produced with energies that are similar to the mean energy of the fluid ( (cid:104) E (cid:105) ∼ T ). – 21 –. Determine whether chemical potentials can be neglected or not.
Chemical potentials in the electromagnetic sector of the plasma can safely be ignoredas a result of the small baryon-to-photon ratio. Since we consider no primordial asym-metry for neutrinos we can also set µ ν = µ ¯ ν = 0. We can also safely set µ φ = 0 sincethe rest of the particles have µ = 0 and none of the processes φ participates in cangenerate a net chemical potential.3. Include all relevant interactions among sectors.
We should calculate the relevant annihilation and scattering energy density transferrates. The rates for ν - e interactions in the Standard Model can be found in (3.4). The φφ ↔ ¯ νν rate for this scenario can be obtained from Equation (2.15) given σ ( s ) = λ πs log (cid:0) s/m ν (cid:1) . The energy transfer rate for νφ ↔ νφ processes cannot be easilycalculated. However, it can be neglected on the grounds that scattering interactionstransfer substantially less energy than annihilation interactions, see Section 2.4.4. Write down evolution equations for T and µ for each sector within the scenario. Since chemical potentials are neglected for all the species, the temperature evolution ofeach sector is given by Equation (2.8) sourced with the relevant energy transfer rates.The actual system of differential equations for T γ , T ν and T φ can be found in [71]. Precision measurements of key cosmological observables such as N eff and the primordial ele-ment abundances represent a confirmation of the thermal history of the Standard Model anda stringent constraint on many of its extensions. In this work – motivated by the complexityof accurately modelling early Universe thermodynamics and by building upon [1] – we havedeveloped a simple, accurate and fast approach to calculate early Universe thermodynamicsin the Standard Model and beyond.In Section 2, we have detailed the approximations upon which the method is basedon and developed the relevant equations that track the early Universe thermodynamics ofany given system. In Section 3, the method was applied to study neutrino decoupling inthe Standard Model. In Section 4, we have used this approach to model the early Universethermodynamics of a BSM scenario featuring a light neutrinophilic boson. We have foundexcellent agreement between our approach and previous literature on the SM, and betweenour BSM results and those obtained by solving the exact Liouville equation. In Section 5, wehave discussed theoretical arguments supporting why the approach presented in this studyis precise and also have discussed its limitations.The main results obtained in this study are: • Neutrino decoupling in the Standard Model can be accurately described with perfectFermi-Dirac distribution functions for the neutrinos that evolve according to (3.3). Byaccounting for finite temperature corrections, and for spin-statistics and the electronmass in the ν - ν and ν - e interaction rates we find N SMeff = 3 . • By solving the exact Liouville equation we have explicitly demonstrated that thethermodynamic evolution of a very light (1 eV < m φ < λ (cid:46) − ) neutrinophilic scalar can be accurately described by Equations (2.9).– 22 –hese results allow us to draw the main conclusion of this work: • The early Universe thermodynamics of any given system in which departures from ther-mal equilibrium are not severe can be accurately modelled by a few simple differentialequations for the temperature and chemical potentials describing each of the relevantspecies. These equations are (2.9). One can (and should) use Equation (2.8) if chemicalpotentials can be neglected. These equations are simple and fast to solve, and includinginteractions among the various particles is straightforward as described in Section 2. Arecipe to implement the method in generic BSM scenarios is outlined in Section 6.To summarize, the method developed in this study and in [1] can be used to accurately trackthe early Universe evolution in the Standard Model and beyond. Thereby enabling one tofind precise predictions for N eff and the primordial element abundances in BSM theories. Wenote that the method has already been applied to study: i) the BBN effect of MeV-scalethermal dark sectors [1, 72], ii) the impact on BBN of invisible neutrino decays [71], iii) thecosmological consequences of light and weakly coupled flavorful gauge bosons [73], iv) N eff constraints on dark photons [74], v) variations of Newton’s constant at the time of BBN [75],and vi) the CMB implications of sub-MeV neutrinophilic scalars [46].We expect the proposed method to prove useful to study many other extensions of theStandard Model and in other contexts not necessarily restricted to neutrino decoupling orBBN. For example, the method could be used to study: late dark matter decoupling [76],generic dark sectors [77, 78], dark sectors equilibrating during BBN [79], or the BBN era inlow-scale Baryogenesis set-ups [80–83].We conclude by highlighting that we publicly release a Mathematica and
Python code:NUDEC BSM containing the codes developed in this study. The code is fast ( t CPU = O (10) s)and we believe that it could serve as a useful tool for particle phenomenologists and cosmol-ogists interested in calculating the early Universe thermodynamics of BSM scenarios. Acknowledgments
I am grateful to Sam Witte, Chris McCabe, Sergio Pastor, Stefano Gariazzo and Pablo F. deSalas for their very helpful comments and suggestions over a draft version of this paper, andto Toni Pich for useful correspondence. This work is supported by the European ResearchCouncil under the European Union’s Horizon 2020 program (ERC Grant Agreement No648680 DARKHORIZONS). v2 Bonus: N SMeff with NLO e + e − ↔ ¯ νν rates In our N SMeff calculation in Section 3 we have accounted for the e correction to the energyand pressure density of the electromagnetic plasma. The reader may be wondering what isthe impact of radiative corrections to the ν − e interaction rates. The weak corrections to therates are tiny since G F E = G F (10 MeV) < − . However, at finite temperature and fortemperatures relevant for neutrino decoupling, QED corrections have been shown to reducethe e + e − ↔ ¯ νν rates by ∼ −
4% [84]. We estimate the impact on N eff of such corrections byusing NUDEC BSM and by approximating the NLO results of [84] with a simple formula (andextrapolating for T (cid:38) N SMeff including such corrections isexpected to be shifted by (cid:39) − . δρδt (cid:12)(cid:12)(cid:12)(cid:12) NLO e + e − ↔ ¯ νν = δρδt (cid:12)(cid:12)(cid:12)(cid:12) LO e + e − ↔ ¯ νν (cid:20) − .
038 log (cid:18) T γ .
11 MeV (cid:19)(cid:21) −→ N SMeff (cid:12)(cid:12)
NLO − N SMeff (cid:12)(cid:12) LO (cid:39) − . , and hence, it seems that using LO ν − e rates limits the accuracy of N SMeff to be (cid:39) . eferences [1] M. Escudero, Neutrino decoupling beyond the Standard Model: CMB constraints on the DarkMatter mass with a fast and precise N eff evaluation , JCAP (2019) 007 [ ].[2]
Planck collaboration, Y. Akrami et al.,
Planck 2018 results. I. Overview and thecosmological legacy of Planck , .[3] Planck collaboration, N. Aghanim et al.,
Planck 2018 results. VI. Cosmological parameters , .[4] SPT-3G collaboration, B. A. Benson et al.,
SPT-3G: A Next-Generation Cosmic MicrowaveBackground Polarization Experiment on the South Pole Telescope , Proc. SPIE Int. Soc. Opt.Eng. (2014) 91531P [ ].[5]
Simons Observatory collaboration, J. Aguirre et al.,
The Simons Observatory: Sciencegoals and forecasts , JCAP (2019) 056 [ ].[6] K. Abazajian et al.,
CMB-S4 Science Case, Reference Design, and Project Plan , .[7] NASA PICO collaboration, S. Hanany et al.,
PICO: Probe of Inflation and Cosmic Origins , .[8] CORE collaboration, E. Di Valentino et al.,
Exploring cosmic origins with CORE:Cosmological parameters , JCAP (2018) 017 [ ].[9] N. Sehgal et al.,
CMB-HD: An Ultra-Deep, High-Resolution Millimeter-Wave Survey OverHalf the Sky , .[10] ParticleDataGroup collaboration, M. Tanabashi et al.,
Review of Particle Physics , Phys.Rev.
D98 (2018) 030001.[11] R. J. Cooke, M. Pettini, K. M. Nollett and R. Jorgenson,
The primordial deuterium abundanceof the most metal-poor damped Ly α system , Astrophys. J. (2016) 148 [ ].[12] E. B. Grohs, J. R. Bond, R. J. Cooke, G. M. Fuller, J. Meyers and M. W. Paris,
Big BangNucleosynthesis and Neutrino Cosmology , .[13] M. Pospelov and J. Pradler, Big Bang Nucleosynthesis as a Probe of New Physics , Ann. Rev.Nucl. Part. Sci. (2010) 539 [ ].[14] F. Iocco, G. Mangano, G. Miele, O. Pisanti and P. D. Serpico, Primordial Nucleosynthesis:from precision cosmology to fundamental physics , Phys. Rept. (2009) 1 [ ].[15] S. Sarkar,
Big bang nucleosynthesis and physics beyond the standard model , Rept. Prog. Phys. (1996) 1493 [ hep-ph/9602260 ].[16] A. D. Dolgov, S. H. Hansen, G. Raffelt and D. V. Semikoz, Heavy sterile neutrinos: Boundsfrom big bang nucleosynthesis and SN1987A , Nucl. Phys.
B590 (2000) 562 [ hep-ph/0008138 ].[17] M. Kawasaki, K. Kohri, T. Moroi and Y. Takaesu,
Revisiting Big-Bang NucleosynthesisConstraints on Long-Lived Decaying Particles , Phys. Rev.
D97 (2018) 023502 [ ].[18] T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram and S. Hannestad,
MeV-scale reheating temperature and thermalization of oscillating neutrinos by radiative andhadronic decays of massive particles , JCAP (2019) 012 [ ].[19] D. Cadamuro, S. Hannestad, G. Raffelt and J. Redondo,
Cosmological bounds on sub-MeVmass axions , JCAP (2011) 003 [ ].[20] M. Cicoli, J. P. Conlon and F. Quevedo,
Dark radiation in LARGE volume models , Phys.Rev.
D87 (2013) 043520 [ ].[21] C. Brust, D. E. Kaplan and M. T. Walters,
New Light Species and the CMB , JHEP (2013)058 [ ]. – 24 –
22] H. Vogel and J. Redondo,
Dark Radiation constraints on minicharged particles in models witha hidden photon , JCAP (2014) 029 [ ].[23] C. Boehm, M. J. Dolan and C. McCabe,
A Lower Bound on the Mass of Cold Thermal DarkMatter from Planck , JCAP (2013) 041 [ ].[24] R. J. Wilkinson, A. C. Vincent, C. Bhm and C. McCabe,
Ruling out the light weaklyinteracting massive particle explanation of the Galactic 511 keV line , Phys. Rev.
D94 (2016)103525 [ ].[25] S. Weinberg,
Goldstone Bosons as Fractional Cosmic Neutrinos , Phys. Rev. Lett. (2013)241301 [ ].[26] P. F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti,
Bounds on verylow reheating scenarios after Planck , Phys. Rev.
D92 (2015) 123534 [ ].[27] S. Gariazzo, P. F. de Salas and S. Pastor,
Thermalisation of sterile neutrinos in the earlyUniverse in the 3+1 scheme with full mixing matrix , JCAP (2019) 014 [ ].[28] P. F. de Salas and S. Pastor,
Relic neutrino decoupling with flavour oscillations revisited , JCAP (2016) 051 [ ].[29] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti and P. D. Serpico,
Relic neutrinodecoupling including flavor oscillations , Nucl. Phys.
B729 (2005) 221 [ hep-ph/0506164 ].[30] A. D. Dolgov,
Neutrinos in cosmology , Phys. Rept. (2002) 333 [ hep-ph/0202122 ].[31] D. A. Dicus, E. W. Kolb, A. M. Gleeson, E. C. G. Sudarshan, V. L. Teplitz and M. S. Turner,
Primordial Nucleosynthesis Including Radiative, Coulomb, and Finite TemperatureCorrections to Weak Rates , Phys. Rev.
D26 (1982) 2694.[32] S. Hannestad and J. Madsen,
Neutrino decoupling in the early universe , Phys. Rev.
D52 (1995) 1764 [ astro-ph/9506015 ].[33] S. Dodelson and M. S. Turner,
Nonequilibrium neutrino statistical mechanics in the expandinguniverse , Phys. Rev.
D46 (1992) 3372.[34] A. D. Dolgov, S. H. Hansen and D. V. Semikoz,
Nonequilibrium corrections to the spectra ofmassless neutrinos in the early universe , Nucl. Phys.
B503 (1997) 426 [ hep-ph/9703315 ].[35] S. Esposito, G. Miele, S. Pastor, M. Peloso and O. Pisanti,
Nonequilibrium spectra ofdegenerate relic neutrinos , Nucl. Phys.
B590 (2000) 539 [ astro-ph/0005573 ].[36] G. Mangano, G. Miele, S. Pastor and M. Peloso,
A Precision calculation of the effectivenumber of cosmological neutrinos , Phys. Lett.
B534 (2002) 8 [ astro-ph/0111408 ].[37] J. Birrell, C.-T. Yang and J. Rafelski,
Relic Neutrino Freeze-out: Dependence on NaturalConstants , Nucl. Phys.
B890 (2014) 481 [ ].[38] E. Grohs, G. M. Fuller, C. T. Kishimoto, M. W. Paris and A. Vlasenko,
Neutrino energytransport in weak decoupling and big bang nucleosynthesis , Phys. Rev.
D93 (2016) 083522[ ].[39] J. Froustey and C. Pitrou,
Incomplete neutrino decoupling effect on big bang nucleosynthesis , Phys. Rev. D (2020) 043524 [ ].[40] A. F. Heckler,
Astrophysical applications of quantum corrections to the equation of state of aplasma , Phys. Rev.
D49 (1994) 611.[41] N. Fornengo, C. W. Kim and J. Song,
Finite temperature effects on the neutrino decoupling inthe early universe , Phys. Rev.
D56 (1997) 5123 [ hep-ph/9702324 ].[42] J. J. Bennett, G. Buldgen, M. Drewes and Y. Y. Y. Wong,
Towards a precision calculation ofthe effective number of neutrinos N eff in the Standard Model I: The QED equation of state , . – 25 –
43] S. Hannestad,
Oscillation effects on neutrino decoupling in the early universe , Phys. Rev.
D65 (2002) 083006 [ astro-ph/0111423 ].[44] A. D. Dolgov, S. H. Hansen, S. Pastor, S. T. Petcov, G. G. Raffelt and D. V. Semikoz,
Cosmological bounds on neutrino degeneracy improved by flavor oscillations , Nucl. Phys.
B632 (2002) 363 [ hep-ph/0201287 ].[45] Z. Chacko, L. J. Hall, T. Okui and S. J. Oliver,
CMB signals of neutrino mass generation , Phys. Rev.
D70 (2004) 085008 [ hep-ph/0312267 ].[46] M. Escudero and S. J. Witte,
A CMB Search for the Neutrino Mass Mechanism and itsRelation to the H Tension , Eur. Phys. J.
C80 (2020) 294 [ ].[47] J. Bernstein,
Kinetic Theory in the Expanding Universe . Cambridge University Press,Cambridge, U.K., 1988, 10.1017/CBO9780511564185.[48] E. W. Kolb and M. S. Turner,
The Early Universe , Front. Phys. (1990) 1.[49] L. D. Landau and E. M. Lifshitz, Statistical Physics, Part 1 , vol. 5 of
Course of TheoreticalPhysics . Butterworth-Heinemann, Oxford, 1980.[50] M. Kawasaki, G. Steigman and H.-S. Kang,
Cosmological evolution of an early decayingparticle , Nucl. Phys.
B403 (1993) 671.[51] G. D. Starkman, N. Kaiser and R. A. Malaney,
Mixed dark matter from neutrino lasing , Astrophys. J. (1994) 12 [ astro-ph/9312020 ].[52] T. Bringmann and S. Hofmann,
Thermal decoupling of WIMPs from first principles , JCAP (2007) 016 [ hep-ph/0612238 ].[53] T. Bringmann,
Particle Models and the Small-Scale Structure of Dark Matter , New J. Phys. (2009) 105027 [ ].[54] T. Binder, L. Covi, A. Kamada, H. Murayama, T. Takahashi and N. Yoshida, Matter PowerSpectrum in Hidden Neutrino Interacting Dark Matter Models: A Closer Look at the CollisionTerm , JCAP (2016) 043 [ ].[55] A. Fradette, M. Pospelov, J. Pradler and A. Ritz,
Cosmological beam dump: constraints ondark scalars mixed with the Higgs boson , Phys. Rev.
D99 (2019) 075004 [ ].[56] C. D. Kreisch, F.-Y. Cyr-Racine and O. Dore,
The Neutrino Puzzle: Anomalies, Interactions,and Cosmological Tensions , .[57] J. Erler and S. Su, The Weak Neutral Current , Prog. Part. Nucl. Phys. (2013) 119[ ].[58] O. Pisanti, A. Cirillo, S. Esposito, F. Iocco, G. Mangano, G. Miele et al., PArthENoPE:Public Algorithm Evaluating the Nucleosynthesis of Primordial Elements , Comput. Phys.Commun. (2008) 956 [ ].[59] R. Consiglio, P. F. de Salas, G. Mangano, G. Miele, S. Pastor and O. Pisanti,
PArthENoPEreloaded , Comput. Phys. Commun. (2018) 237 [ ].[60] Y. Chikashige, R. N. Mohapatra and R. D. Peccei,
Are There Real Goldstone BosonsAssociated with Broken Lepton Number? , Phys. Lett. (1981) 265.[61] G. B. Gelmini and M. Roncadelli,
Left-Handed Neutrino Mass Scale and SpontaneouslyBroken Lepton Number , Phys. Lett. (1981) 411.[62] J. Schechter and J. W. F. Valle,
Neutrino Decay and Spontaneous Violation of LeptonNumber , Phys. Rev.
D25 (1982) 774.[63] H. M. Georgi, S. L. Glashow and S. Nussinov,
Unconventional Model of Neutrino Masses , Nucl. Phys.
B193 (1981) 297. – 26 –
64] A. D. Dolgov, S. H. Hansen, S. Pastor and D. V. Semikoz,
Unstable massive tau neutrinosand primordial nucleosynthesis , Nucl. Phys.
B548 (1999) 385 [ hep-ph/9809598 ].[65] A. D. Dolgov, S. H. Hansen and D. V. Semikoz,
Nonequilibrium corrections to the spectra ofmassless neutrinos in the early universe: Addendum , Nucl. Phys.
B543 (1999) 269[ hep-ph/9805467 ].[66] S. Dodelson and L. M. Widrow,
Sterile-neutrinos as dark matter , Phys. Rev. Lett. (1994)17 [ hep-ph/9303287 ].[67] X.-D. Shi and G. M. Fuller, A New dark matter candidate: Nonthermal sterile neutrinos , Phys. Rev. Lett. (1999) 2832 [ astro-ph/9810076 ].[68] T. Venumadhav, F.-Y. Cyr-Racine, K. N. Abazajian and C. M. Hirata, Sterile neutrino darkmatter: Weak interactions in the strong coupling epoch , Phys. Rev.
D94 (2016) 043515[ ].[69] L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West,
Freeze-In Production of FIMPDark Matter , JHEP (2010) 080 [ ].[70] C. Dvorkin, T. Lin and K. Schutz, Making dark matter out of light: freeze-in from plasmaeffects , Phys. Rev.
D99 (2019) 115009 [ ].[71] M. Escudero and M. Fairbairn,
Cosmological Constraints on Invisible Neutrino DecaysRevisited , Phys. Rev.
D100 (2019) 103531 [ ].[72] N. Sabti, J. Alvey, M. Escudero, M. Fairbairn and D. Blas,
Refined Bounds on MeV-scaleThermal Dark Sectors from BBN and the CMB , JCAP (2020) 004 [ ].[73] M. Escudero, D. Hooper, G. Krnjaic and M. Pierre,
Cosmology with a very light L µ - L τ gauge boson , JHEP (2019) 071 [ ].[74] M. Ibe, S. Kobayashi, Y. Nakayama and S. Shirai, Cosmological Constraint on Dark Photonfrom N eff , .[75] J. Alvey, N. Sabti, M. Escudero and M. Fairbairn, Improved BBN Constraints on theVariation of the Gravitational Constant , Eur. Phys. J. C (2020) 148 [ ].[76] J. A. D. Diacoumis and Y. Y. Y. Wong, Trading kinetic energy: how late kinetic decoupling ofdark matter changes N eff , JCAP (2019) 001 [ ].[77] M. A. Buen-Abad, G. Marques-Tavares and M. Schmaltz,
Non-Abelian dark matter and darkradiation , Phys. Rev.
D92 (2015) 023531 [ ].[78] Z. Chacko, Y. Cui, S. Hong and T. Okui,
Hidden dark matter sector, dark radiation, and theCMB , Phys. Rev.
D92 (2015) 055033 [ ].[79] A. Berlin, N. Blinov and S. W. Li,
Dark Sector Equilibration During Nucleosynthesis , Phys.Rev.
D100 (2019) 015038 [ ].[80] K. Aitken, D. McKeen, T. Neder and A. E. Nelson,
Baryogenesis from Oscillations ofCharmed or Beautiful Baryons , Phys. Rev.
D96 (2017) 075009 [ ].[81] G. Elor, M. Escudero and A. Nelson,
Baryogenesis and Dark Matter from B Mesons , Phys.Rev.
D99 (2019) 035031 [ ].[82] A. E. Nelson and H. Xiao,
Baryogenesis from B Meson Oscillations , Phys. Rev.
D100 (2019)075002 [ ].[83] G. Kane and M. W. Winkler,
Baryogenesis from a Modulus Dominated Universe , .[84] S. Esposito, G. Mangano, G. Miele, I. Picardi and O. Pisanti, Neutrino energy loss rate in astellar plasma , Nucl. Phys. B (2003) 217 [ astro-ph/0301438 ].[85] K. Enqvist, K. Kainulainen and V. Semikoz,
Neutrino annihilation in hot plasma , Nucl. Phys.
B374 (1992) 392. – 27 –
86] D. J. Fixsen,
The Temperature of the Cosmic Microwave Background , Astrophys. J. (2009) 916 [ ].[87] S. Dodelson,
Modern Cosmology . Academic Press, Amsterdam, 2003.[88] C. Pitrou, A. Coc, J.-P. Uzan and E. Vangioni,
Precision big bang nucleosynthesis withimproved Helium-4 predictions , Phys. Rept. (2018) 1 [ ].[89] B. D. Fields, S. Dodelson and M. S. Turner,
Effect of neutrino heating on primordialnucleosynthesis , Phys. Rev.
D47 (1993) 4309 [ astro-ph/9210007 ].[90] P. D. Serpico, S. Esposito, F. Iocco, G. Mangano, G. Miele and O. Pisanti,
Nuclear reactionnetwork for primordial nucleosynthesis: A Detailed analysis of rates, uncertainties and lightnuclei yields , JCAP (2004) 010 [ astro-ph/0408076 ].[91] A. Cuoco, J. Lesgourgues, G. Mangano and S. Pastor,
Do observations prove that cosmologicalneutrinos are thermally distributed? , Phys. Rev.
D71 (2005) 123501 [ astro-ph/0502465 ].[92] P. J. E. Peebles,
Recombination of the Primeval Plasma , Astrophys. J. (1968) 1.[93] R. Weymann,
Diffusion Approximation for a Photon Gas Interacting with a Plasma via theCompton Effect , Physics of Fluids (1965) 2112.[94] W. T. Hu, Wandering in the Background: A CMB Explorer , Ph.D. thesis, UC, Berkeley, 1995. astro-ph/9508126 .[95] C.-P. Ma and E. Bertschinger,
Cosmological perturbation theory in the synchronous andconformal Newtonian gauges , Astrophys. J. (1995) 7 [ astro-ph/9506072 ].[96] A. Arbey, J. Auffinger, K. P. Hickerson and E. S. Jenssen,
AlterBBN v2: A public code forcalculating Big-Bang nucleosynthesis constraints in alternative cosmologies , .[97] A. Baha Balantekin and B. Kayser, On the Properties of Neutrinos , Ann. Rev. Nucl. Part.Sci. (2018) 313 [ ].[98] P. Gondolo and G. Gelmini, Cosmic abundances of stable particles: Improved analysis , Nucl.Phys.
B360 (1991) 145.[99] J. Edsjo and P. Gondolo,
Neutralino relic density including coannihilations , Phys. Rev.
D56 (1997) 1879 [ hep-ph/9704361 ].[100] M. Kawasaki, K. Kohri and N. Sugiyama,
MeV scale reheating temperature andthermalization of neutrino background , Phys. Rev.
D62 (2000) 023506 [ astro-ph/0002127 ]. – 28 – Appendices
In these appendices we provide: details about Standard Model interaction rates includingspin-statistics and the electron mass A.1. A detailed comparison between our SM resultsand previous SM calculations A.2. The relevant equations and results from solving neutrinodecoupling in the SM including neutrino chemical potentials A.3. A comparison between thesolution to the Liouville equation and the proposed method for a very light and weakly cou-pled neutrinophilic scalar A.4. The temperature and chemical potential evolution equationsin the Maxwell-Boltzmann approximation A.5. Thermodynamic formulae A.6. An explicitcalculation of the number and energy density transfer rates in the Maxwell-Boltzmann ap-proximation for decay, annihilation and scattering processes together with spin-statisticscorrections to them A.7.
A.1 SM ν - e and ν - ν interaction rates In this appendix we outline the results of the neutrino-electron and neutrino-neutrino en-ergy and number transfer rates including spin-statistics and a non-negligible electron mass.For that purpose, we have calculated the exact collision terms following the phase spaceintegration method of [32] (see also the appendices of [55] and [56]).
Interaction rates in the Maxwell-Boltzmann approximation and m e = 0We begin by outlining the results obtained in [1], the reader is referred to Appendix A.2 ofthat reference for more details. In the Maxwell-Boltzmann approximation and by neglectingthe electron mass, the neutrino energy density transfer rates read: δρ ν e δt (cid:12)(cid:12)(cid:12)(cid:12) MBSM = G F π (cid:2) (cid:0) g eL + g eR (cid:1) F MB ( T γ , T ν e ) + 2 F MB ( T ν µ , T ν e ) (cid:3) , (A.1a) δρ ν µ δt (cid:12)(cid:12)(cid:12)(cid:12) MBSM = G F π (cid:2) (cid:0) g µL + g µR (cid:1) F MB ( T γ , T ν µ ) − F MB ( T ν µ , T ν e ) (cid:3) , (A.1b)where g eL , g eR , g µL , g µR are defined in Equation (3.5) and where we have defined: F MB ( T , T ) = 32 ( T − T ) + 56 T T ( T − T ) . (A.2)Similarly, the number density transfer rates read: δn ν e δt (cid:12)(cid:12)(cid:12)(cid:12) MBSM = 8 G F π (cid:104) (cid:0) g eL + g eR (cid:1) ( T γ − T ν e ) + 2 ( T ν µ − T ν e ) (cid:105) , (A.3a) δn ν µ δt (cid:12)(cid:12)(cid:12)(cid:12) MBSM = 8 G F π (cid:104) (cid:0) g µL + g µR (cid:1) ( T γ − T ν µ ) − ( T ν µ − T ν e ) (cid:105) , (A.3b)where δn i /δt ≡ (cid:80) j (cid:104) σv (cid:105) ij ( n i − n j ). Electron-neutrino rates with quantum statistics
By numerically evaluating the energy transfer rates between neutrinos and electrons we havefound that the rates including quantum statistics are: δρ ν e δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = G F π (cid:2) (cid:0) g eL + g eR (cid:1) F ( T γ , T ν e ) + 2 F ( T ν µ , T ν e ) (cid:3) , (A.4a) δρ ν µ δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = G F π (cid:2) (cid:0) g µL + g µR (cid:1) F ( T γ , T ν µ ) − F ( T ν µ , T ν e ) (cid:3) , (A.4b)– 29 – igure 8 . Neutrino-electron energy ( left panel ) and number ( right panel ) density transfer ratestaking into account spin statistics in the collision term and the electron mass as compared to the onesobtained in the Maxwell-Boltzmann approximation with m e = 0 (A.1). The dotted lines correspondto scatterings and the dashed lines correspond to annihilations. Clearly, the impact of the electronmass is small for the relevant range of temperatures. Yet, we take it into account in NUDEC BSM. where we have defined: F ( T , T ) = 32 f FD a ( T − T ) + 56 f FD s T T ( T − T ) , (A.5)where f FD a = 0 . f FD s = 0 . f FD a and f FD s account for thePauli blocking effects in the rates for annihilations and scatterings processes, respectively.For the number density transfer rates we find: δn ν δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = 0 . × δn ν δt (cid:12)(cid:12)(cid:12)(cid:12) MBSM . (A.6)The annihilation cross section taking into account Fermi-Dirac statistics was first calculatedin [85]. Since (cid:104) σv (cid:105) ν = ( (cid:104) σv (cid:105) ν e + 2 (cid:104) σv (cid:105) ν µ ) ≡ n ν ( T γ ) δnδt (cid:12)(cid:12)(cid:12) T ν =0 , we find the following neutrinoannihilation cross sections: (cid:104) σv (cid:105) MB¯ νν ↔ e + e − = 3 . G F T , (A.7) (cid:104) σv (cid:105) FD¯ νν ↔ e + e − = 2 . G F T . (A.8)In the MB approximation we find an annihilation cross section that is 6% higher than the onereported in [85]. The annihilation cross section taking into account Pauli blocking preciselymatches the one reported in [85].The Fermi-Dirac suppression factors for the energy transfer rates in (3.6) of 0.884 and0.829 have been calculated assuming that ( T γ − T ν ) /T γ ≡ ∆ T /T = 0 .
01. However, we haveexplicitly checked that even if ∆
T /T = ± . T (cid:38) ν - e and ν - ν rates in any realistic BSM temperature evolution scenario.– 30 – lectron-neutrino rates with quantum statistics and the electron mass We have numerically calculated the effect of the electron mass in the energy and numberdensity transfer rates. The effect of the electron mass is small when the e - ν interactions areefficient in the early Universe, namely for T > N SMeff , a non-negligible electron mass leads to a decrease of (cid:39) − . T γ − T ν ) /T γ ≡ ∆ T /T = 0 .
01. We have explicitly checked that the resulting rates consideringthe actual T ν ( T γ ) evolution in the early Universe does not lead to relevant changes in N eff . A.2 Comparison with previous literature on the SM: supplementary material
In this appendix we provide detailed comparison between our calculation of neutrino decou-pling with previous accurate calculations. We compare at the level of the energy carried byeach neutrino species, neutrino number density, entropy density of the Universe, yields of theprimordial elements, the resulting frozen neutrino number and energy contributions, and thetime evolution of the entropy released from out of equilibrium e - ν interactions.In order to make contact with previous literature, it is useful not to work in terms of T γ and T ν but to do so in terms of comoving dimensionless ratios. We define z γ ≡ a T γ /m , z ν ≡ a T ν /m , (A.9)where a is the dimensionless scale factor, and m is a normalization factor with dimensionsof energy. For concreteness, we shall choose m = m e . Given these definitions we can easilyfind their scale factor evolution: dzda = za (cid:20) H T dTdt (cid:21) . (A.10)By plugging Equations (3.2) and (3.3) into (A.10) we can explicitly write down theevolution equation for z γ and z ν in the Standard Model: dz ν α da = z ν α a δρ να δt Hρ ν α , (A.11a) dz γ da = z γ a (cid:104) T γ dρ e dT γ − p e + ρ e ) (cid:105) + T γ (cid:104) T γ d P int dT γ − dP int dT γ (cid:105) − H δρ ν δt T γ (cid:16) T γ d P int dT γ + dρ γ dT γ + dρ e dT γ (cid:17) . (A.11b)We solve these equations in the SM by taking into account all relevant interactions and NLOfinite temperature corrections. We use the same initial conditions as in Section 3.4 and findthat once neutrinos have completely decoupled and electrons and positrons annihilated away: T ν e (cid:54) = T ν µ,τ → z γ = 1 . , z ν e = 1 . , z ν µ,τ = 1 . , → N eff = 3 . , (A.12) T ν e = T ν µ,τ → z γ = 1 . , z ν e = 1 . , z ν µ,τ = 1 . , → N eff = 3 . . (A.13) N eff In Table 3 we compare the results from our solution to neutrino decoupling to those of thelatest and most accurate determination in the SM by de Salas & Pastor [28], and also to thoseof Dolgov et. al. [65] and Grohs et. al. [38] that do not account for neutrino oscillations or– 31 –eutrino Decoupling in the Standard ModelNo oscillations, no QED N eff z γ δρ ν e /ρ ν e (%) δρ ν µ /ρ ν µ (%) δρ ν τ /ρ ν τ (%)This work 3.0350 1.39903 0.971 0.407 0.407Dolgov et. al. [65] 3.0340 1.39910 0.946 0.398 0.398Grohs et. al. [38] 3.0340 1.39904 0.928 0.377 0.377No oscillations, LO-QED N eff z γ δρ ν e /ρ ν e (%) δρ ν µ /ρ ν µ (%) δρ ν τ /ρ ν τ (%)This work 3.0453 1.39782 0.959 0.401 0.401de Salas & Pastor [28] 3.0446 1.39784 0.920 0.392 0.392With oscillations, LO-QED N eff z γ δρ ν /ρ ν (%) δρ ν /ρ ν (%) δρ ν /ρ ν (%)This work 3.0462 1.39777 0.603 0.603 0.603NH, de Salas & Pastor [28] 3.0453 1.39779 0.636 0.574 0.519IH, de Salas & Pastor [28] 3.0453 1.39779 0.635 0.574 0.520 Table 3 . N SMeff as obtained in this work and compared with the calculations of de Salas & Pastor [28]that includes neutrino oscillations and LO-finite temperature corrections, and that of Dolgov et.al. [65] and of Grohs et. al. [38] that do not. When neutrino oscillations are neglected, the differencein N SMeff is 0 . .
001 as compared to [38, 65]. In our case, δρ ν α /ρ ν α ≡ z ν α − T ν e = T ν µ = T ν τ ) the difference in terms of N SMeff is 0 . finite temperature corrections. We can appreciate that the resulting values of N eff only differby at most 0 . . .
03% as compared to the results of Ref. [65].This appears to be a consequence of the starting point of the derivation of our temperaturetime evolution equations: energy conservation (2.5). We therefore estimate the accuracy ofour reported value of N SMeff to be 0 . ν h A parameter that is of relevance to cosmology is the energy density encoded in non-relativisticneutrino species: Ω ν h . Since we describe neutrinos with a perfect Fermi-Dirac distributionfunction, we find that for degenerate and non-relativistic neutrinos:Ω ν h ≡ (cid:80) ν m ν n ν ρ c /h = (cid:80) ν m ν ρ c /h ζ (3) π T ν = (cid:80) ν m ν .
05 eV , (A.14)where we have used T γ /T ν = 1 . T γ = 2 . ν h = (cid:80) ν m ν / (93 .
12 eV). An instantaneous neutrino decouplingwould correspond to Ω ν h = (cid:80) ν m ν / (94 .
11 eV) [87].When neutrino chemical potentials are allowed to vary (see Appendix A.3), we find thatthe resulting number density encoded in non-relativistic degenerate neutrinos is:Ω ν h = (cid:80) ν m ν .
14 eV , (A.15)and hence in excellent agreement with the results of [29].– 32 –rimordial nuclei abundances in the Standard ModelRelative uncertainty ∆ rel Y p ∆ rel (D/H) ∆ rel ( He / H) ∆ rel ( Li / H)Neutrino evolution, T ν e = T ν µ × − × − × − − × − Neutrino evolution, T ν e (cid:54) = T ν µ × − × − × − − × − Measurements, PDG [10] 1 . × − − - 1 . × − Nuclear rates,
PRIMAT [88] 6 . × − . × − . × − . × − Table 4 . ∆ rel X ≡ σ X /X . The first two rows correspond to the relative change in the primordialnuclei abundances using the thermodynamic evolution obtained in this work to that compared tothe default neutrino evolution used in PArthENoPE which is based on the SM calculation of [29]. Weappreciate that the relative difference between the evolution used in
PArthENoPE and that obtainedhere yields differences in the primordial element abundances that are at least one order of magnitudesmaller than the error associated with current measurements or with nuclear uncertainties.
Primordial nuclei abundances
Neutrino decoupling occurs at T ∼ T ∼ . i) neutrinos participate in proton-to-neutronreactions and hence the neutrino number density enters the proton-to-neutron conversionrates, ii) neutrinos affect the expansion rate of the Universe, and iii) residual electron-neutrino interactions generate a net amount of entropy in the Universe.State-of-the-art BBN codes, such as PArthENoPE [58, 59] and now also
PRIMAT [88],account for each of these effects by taking the relevant neutrino evolution in the SM asobtained in [29].Here, we compare the relative changes to the nuclei abundances of some relevant nucleigiven our neutrino evolution by implementing it in
PArthENoPE . In Table 4 we compare therelative changes in the nuclei abundances from our neutrino evolution to those obtained inthe default mode of
PArthENoPE . We can clearly appreciate that the relative changes to theprimordial nuclei abundances from our evolution in the SM as compared to the default modeof
PArthENoPE are smaller – at least by an order of magnitude – to current observational [10]and nuclear reaction [88] uncertainties. Finally, we note that we provide a file with therelevant Standard Model evolution that could be used in BBN codes.
Entropy
As a result of the out-of-equilibrium energy exchange between electrons and neutrinos atthe time of neutrino decoupling, a small amount of entropy is generated in the Universe,see e.g. [38]. We can parametrize the entropy of the Universe in terms of the usual numberof degrees of freedom contributing to entropy density: g (cid:63)s ≡ s / (2 π T γ ), where s is theentropy density. We obtain: g (cid:63)s = 2 + 214 (cid:18) T ν T γ (cid:19) = 3 . . (A.16)If neutrinos decouple instantaneously [87]: g (cid:63)s = 3 . g (cid:63)s = 3 . g (cid:63)s is accurate at the per mille levelas compared to that obtained from [28]. – 33 – y ν n i n s ν h dn ν d y ν − dn i n s ν d y ν i ( % ) Solid: de Salas & PastorDashed: this work ν µ , τ ν e Number Density, without Oscillations 0 2 4 6 8 10 12 y ν n i n s ν h dn ν d y ν − dn i n s ν d y ν i ( % ) Solid: de Salas & PastorDashed: this work ν ν ν ν Number Density, with Oscillations0 2 4 6 8 10 12 y ν ρ i n s ν h d ρ ν d y ν − d ρ i n s ν d y ν i ( % ) Solid: de Salas & PastorDashed: this work ν µ , τ ν e Energy Density, without Oscillations 0 2 4 6 8 10 12 y ν ρ i n s ν h d ρ ν d y ν − d ρ i n s ν d y ν i ( % ) Solid: de Salas & PastorDashed: this work ν ν ν ν Energy Density, with Oscillations
Figure 9 . Comparison between the results of de Salas & Pastor [28] (solid) and those obtainedhere (dashed). y ν is the comoving momentum that is related to the photon temperature and thephysical neutrino momentum as y ν ≡ z γ × p ν /T γ , where z γ = 1 . f ins ν ( y ) ≡ (1 + e y ) − isthe neutrino distribution function of an instantaneously decoupled neutrino. dn ν /dy ∝ y f ν and dρ ν /dy ∝ y f ν . Our results in the left panels are obtained by evolving separately T ν e and T ν µ, τ , whileour results for the right panels are obtained by solving for a common neutrino temperature, T ν . Theupper panels correspond to the contribution to the number density in neutrino species from residualelectron-positron annihilations. The lower frames correspond to the resulting change in the energydensity. dρ ν /dp ν and dn ν /dp ν Figure 9 shows the neutrino distribution functions at
T < m e /
20 in the case obtained hereand compared to that obtained by solving the full Liouville equation as found in Ref. [28].We appreciate that the resulting non-instantaneous neutrino decoupling corrections to thenumber density in [28] peak at y ν ≡ p ν /T ν ∼
4, while ours do at y ν ∼
3. In the case of dρ ν /dy ν , those from [28] peak at y ν ∼ y ν ∼
4. Our results for dρ ν /dy ν and dn ν /dy ν do not precisely match those of Ref. [28] (as is to be expected, because thedistribution functions that we consider are pure Fermi-Dirac distributions), however, theyare accurate to the 0.05% level for any momenta y ν . In addition, the agreement between thetwo calculations to the total neutrino energy density in each neutrino species is remarkable(see the first five rows of Table 3). – 34 – T γ (MeV) − − − − − [ − ρ ν / ρ PA r t h E N o PE ν ] ( % ) T ν e = T ν µ , τ T ν e = T ν µ , τ Neutrino energy comparison 10 5 3 2 1 0.7 0.4 0.2 T γ (MeV) − − − [ − N / N PA r t h E N o PE ] ( % ) T ν e = T ν µ , τ T ν e = T ν µ , τ Entropy release comparison
Figure 10 . Left panel: neutrino energy evolution in the Standard Model as obtained in this workand compared with that used in the BBN code
PArthENoPE [58, 59]. The difference between the twois at most 0.25%.
Right panel: entropy release as a result of residual electron-positron annihilationsinto neutrinos in the Standard Model. This is encoded in terms of the function N that controls thetime dependence of the energy exchange between electrons and neutrinos in the early Universe. Thedifference between the two is at the level of (cid:39) N in thecase T ν e (cid:54) = T ν µ, τ agree to the ∼
3% level with N as recently calculated in [39] by solving for the fullLiouville equation but neglecting neutrino oscillations. Evolution of entropy release: implications for BBN
We make contact with previous SM neutrino evolution results by considering the variable N , introduced in Ref. [58] ( PArthENoPE ), that accounts for the entropy release of electron-positron annihilation to the neutrino sector.Given our definitions (see Equation A.9), N can be written as: N ( a ) ≡ z γ (cid:32) a dda (cid:34) (cid:18) am (cid:19) ρ ν (cid:35)(cid:33) , (A.17)where ρ ν = ρ ν e + 2 ρ ν µ . We can work out this expression to write it as: N ( a ) = 4 z γ π a (cid:88) α z ν α dz ν α da , (A.18)and hence, it is clear that if δρ ν /δt = 0 then N = 0 (see Equation (A.11)). Ref. [58] providesa fitting function for N as a function of m e /T γ . In our case, m e /T γ = m e m az γ = az γ .The right panels of Figure 10 show the comparison between N as obtained from ourapproach and that used by default in PArthENoPE . We can appreciate that the difference isat the level of 10% or less. The left panels of Figure 10 show the comparison in terms of theneutrino energy density. We see that difference in terms of ρ ν is at most 0.25%. The collectiveeffect of these two variables in terms of the prediction for the primordial abundances is shownin Table 4 where one can see that the change in the relevant nuclei abundances is smaller –at least by an order of magnitude – to current observational and nuclear rates uncertainties.– 35 – .3 N eff in the Standard Model including neutrino chemical potentials In this appendix we consider the case of neutrino decoupling in the Standard Model by al-lowing chemical potentials to vary. We do so by solving equations (2.9a) and (2.9b) for theneutrino temperature and chemical potential respectively. For the electromagnetic compo-nent of the plasma we still neglect chemical potentials since they are highly suppressed by thesmall baryon-to-photon ratio µ e /T ∼ − . Therefore, the photon temperature still evolvesaccording to (3.2). In addition, we consider that there is no primordial lepton asymmetry,so that µ ν = µ ¯ ν . Neutrino-neutrino and electron-neutrino interactions
We include the effect of the neutrino chemical potentials in the neutrino-neutrino and neutrino-electron energy and number density transfer rates, that allowing for neutrino chemical po-tentials explicitly read: δρ ν e δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = G F π (cid:2) (cid:0) g eL + g eR (cid:1) G ( T γ , , T ν e , µ ν e ) + 2 G ( T ν µ , µ ν µ , T ν e , µ ν e ) (cid:3) , (A.19) δρ ν µ δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = G F π (cid:2) (cid:0) g µL + g µR (cid:1) G ( T γ , , T ν µ , µ ν µ ) − G ( T ν µ , µ ν µ , T ν e , µ ν e ) (cid:3) , (A.20) δn ν e δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = 8 f FD n G F π (cid:20) (cid:0) g eL + g eR (cid:1) ( T γ − T ν e e µνeTνe )+2( T ν µ e µνµTνµ − T ν e e µνeTνe ) (cid:21) , (A.21) δn ν µ δt (cid:12)(cid:12)(cid:12)(cid:12) FDSM = 8 f FD n G F π (cid:20) (cid:0) g µL + g µR (cid:1) ( T γ − T ν µ e µνµTνµ ) − ( T ν µ e µνµTνµ − T ν e e µνeTνe ) (cid:21) , (A.22)where f FD n = 0 .
852 and G ( T , µ , T , µ ) = 32 f FD a ( T e µ T − T e µ T ) + 56 f FD s e µ T e µ T T T ( T − T ) . (A.23)The effect of a non-negligible electron mass is also included as in (3.3). Results
We solve Equations (2.9) for the neutrinos and Equation (3.2) for the electromagnetic plasmausing as initial conditions T ν = T γ = 10 MeV and µ ν /T ν = − − . Such choice of initialvalue for the neutrino chemical potentials does not affect any result, since the energy densityin neutrinos scales as ρ ν ∼ e µ ν /T ν .Our results are outlined in Table 5. We can appreciate that when neutrino chemicalpotentials are allowed to vary and their value is dictated by the dynamics, they are stillnegligible: T ν /µ ν ∼ − ρ ν ∼ e µ ν /T ν and1 − e − / ∼ . N eff is small. This thereforejustifies solving for neutrino decoupling assuming vanishing neutrino chemical potentials.When comparing the results in terms of N eff , g (cid:63)S , Ω ν h with previous state-of-the-artcalculations in the Standard Model [28, 29] we find an excellent agreement. Specially, wenotice that allowing chemical potentials to vary yields an even better value of the neutrinonumber density and therefore of Ω ν h . This was to be expected, since neutrino chemicalpotentials arise as a result of the freeze-out of number changing processes.– 36 –eutrino Decoupling in the SM including neutrino chemical potentialsCase T γ /T ν T ν /µ ν N eff g (cid:63)s (cid:80) m ν / Ω ν h Instantaneous 1.40102 ∞ µ ν = 0, LO-QED 1.39567 ∞ µ ν (cid:54) = 0, LO-QED 1.39434 -236.5 3.045 3.931 93.12 eVRefs [28, 29] (LO-QED) - - 3.045 3.933 93.12 eV µ ν = 0, NLO-QED 1.39578 ∞ µ ν (cid:54) = 0, NLO-QED 1.39445 -236.5 Table 5 . The resulting neutrino temperature and chemical potentials as obtained in this work withand without allowing for non-negligible neutrino chemical potentials. We also show the resultingvalues of N eff , g (cid:63)s and (cid:80) m ν / Ω ν h for degenerate neutrinos. A.4 Detailed neutrinophilic scalar thermodynamics
In this appendix we present a detailed comparison between the Full solution to the Liouvilleequation and the Fast approach presented in this work.First, in Figure 11 we show the CPU time used for each computation as a function ofΓ eff . The comparison in terms of ρ φ , ρ ν , n φ and n ν a function of T ν is displayed in Figure 13.We can appreciate and excellent overall agreement.Finally, in Figure 12 we show the comparison between the two approaches at the levelof the neutrino and φ distribution functions. We clearly see that there is good agreementbetween the two although there are some differences. One appreciates that the Fast solutionsoverestimate the energy/number density of φ particles at low momenta y (cid:46) − y (cid:38)
4. However, this mismatch is compensated when integrating overmomenta since the total energy and number densities agree very well in the two cases. Withrespect to the neutrino distribution function, we notice that the relative difference at the levelof dρ ν /dy ν is below 0.3% independently of momenta and of temperature. In particular, wenotice that for T (cid:28) m φ , the neutrino distributions agree with a remarkable 0 .
1% accuracy.Note that these very small differences have no impact for CMB observations [91].
Figure 11 . Left:
CPU time spent by the integrator in solving the Liouville equation for the neutrino- φ system (4.14) (red) and by solving the system of differential equations (4.5) (black). Right: fractionaltime spent by the integrator when solving the Liouville equation for f ν and f φ as a function of T ν . – 37 – y φ ≡ a p φ / MeV10 − − − − − − d ρ φ / T γ ∝ p φ E φ f φ T ν / m φ = 5 T ν / m φ = 2 T ν / m φ = 1 / T ν / m φ = 1 / Γ eff = 10 FullFast y ν ≡ a p ν / MeV − − − − − × d ( ρ ν − ρ F r ee ν ) / ρ F r ee ν T ν / m φ = 5 T ν / m φ = 2 T ν / m φ = 1 / T ν / m φ = 1 / Γ eff = 10 FullFast y φ ≡ a p φ / MeV10 − − − − − − d ρ φ / T γ ∝ p φ E φ f φ T ν / m φ = 5 T ν / m φ = 2 T ν / m φ = 1 / T ν / m φ = 1 / Γ eff = 1 FullFast y ν ≡ a p ν / MeV − − − − − × d ( ρ ν − ρ F r ee ν ) / ρ F r ee ν T ν / m φ = 5 T ν / m φ = 2 T ν / m φ = 1 / T ν / m φ = 1 / Γ eff = 1 FullFast y φ ≡ a p φ / MeV10 − − − − − − d ρ φ / T γ ∝ p φ E φ f φ T ν / m φ = 5 T ν / m φ = 2 T ν / m φ = 1 / T ν / m φ = 1 / Γ eff = 0.1 FullFast y ν ≡ a p ν / MeV − − − − − × d ( ρ ν − ρ F r ee ν ) / ρ F r ee ν T ν / m φ = 5 T ν / m φ = 2 T ν / m φ = 1 / T ν / m φ = 1 / Γ eff = 0.1 FullFast
Figure 12 . Snapshots of f ν and f φ for Γ eff = 0 . , ,
10 as a function of T ν /m φ . Solid lines correspondto the solution of the exact Liouville equation (4.14) while dashed lines correspond to the solutionof the simplified and fast approach to it developed in this study (2.9). Left panel: dρ φ ( y φ ) /T asa function of comoving momenta y . This quantity is proportional to E φ p φ f φ . Right panel: relativedifference between the distribution function of neutrinos as compared with purely freely streamingones with T γ /T ν = 1 . – 38 – − T ν / m φ − − − | ρ F u ll ν − ρ F a s t ν | / ρ F a s t ν ( % ) Γ eff − T ν / m φ − − | n F u ll ν − n F a s t ν | / n F a s t ν ( % ) Γ eff − T ν / m φ − − | ρ F u ll φ − ρ F a s t φ | / ρ F a s t φ ( % ) Γ eff − T ν / m φ − − | n F u ll φ − n F a s t φ | / n F a s t φ ( % ) Γ eff Figure 13 . Relative difference of ρ ν , n ν , ρ φ and n φ as a function of T ν between the full solution to theLiouville equation (4.14) and the simplified and fast approach to it developed in this study (2.9). Theregions of the lines that are dashed represent the cosmologically relevant ones, in which ρ φ /ρ ν > − . A.5 Evolution equations in the Maxwell-Boltzmann approximation
In this appendix we particularize the generic system of Equations (2.9) for the case in whichthe species follows Maxwell-Boltzmann statistics. We explicitly write down the results inthe limits in which T (cid:29) m and T (cid:28) m and compare the temperature evolution of a non-interacting particle with negligible chemical potential with the Bose-Einstein and Fermi-Diraccases. Finally, as an illustrative example, we re-derive the well-known baryon temperatureevolution as relevant for recombination. Simplified Evolution Equations
We can substantially simplify Equations (2.9) for a species that follows MB statistics, sincein that scenario: dn ( T, µ ) dT = ρ − µnT , dn ( T, µ ) dµ = nT , (A.24) dρ ( T, µ ) dT = σ − µρT , dρ ( T, µ ) dµ = ρT , (A.25)– 39 –here we have defined σ ≡ g π (cid:90) ∞ m dE E (cid:112) E − m e − ( E − µ ) /T . (A.26)In this case, we can easily work out the temperature and chemical potential evolution (2.9)to find: dTdt = T (cid:20) − Hnp − δnδt ρ + δρδt n (cid:21)(cid:30) (cid:2) nσ − ρ (cid:3) , (A.27) dµdt = T (cid:20) H ( n ( µp + σ ) − ρ ( p + ρ )) + δnδt ( µρ − σ ) + δρδt ( ρ − µn ) (cid:21) (cid:30) (cid:2) ρ − nσ (cid:3) . (A.28)In the Maxwell-Boltzmann approximation, the quantities n , ρ , p and σ can be written interms of modified Bessel functions as outlined in (A.38). Since we are working in the MBapproximation, Equations (A.27) and (A.28) can be written in a very compact manner forthe case of ultrarelativistic species ( T (cid:29) m ) or non-relativistic ones ( T (cid:28) m ). The T (cid:29) m and T (cid:28) m cases For ultrarelativistic species ( T (cid:29) m ), the evolution equations read: dTdt = − H T + (cid:20) ρ δρδt − n δnδt (cid:21) T , (A.29a) dµdt = − H µ + (cid:20) ρ δρδt − n δnδt (cid:21) µ + (cid:20) n δnδt − ρ δρδt (cid:21) T , (A.29b)where ρ = g π T e µ/T .While, for non-relativistic particles, ( T (cid:28) m ): dTdt = − H T + 23 n (cid:20) δρδt − m δnδt (cid:21) , (A.30a) dµdt = − H ( µ − m ) + 23 n (cid:20) δρδt − m δnδt (cid:21) µ − mT , (A.30b)since nσ − ρ (cid:39) T np , σ (cid:39) mρ , p (cid:39) T n and ρ (cid:39) mn , and n (cid:39) g (cid:0) T m π (cid:1) / e ( µ − m ) /T inthe T (cid:28) m limit. Clearly the set of equations (A.29) and (A.30) reproduce the correctscale factor behaviour in the absence of interactions since d/dt = aHd/da . See for instanceEquations 3.80, 3.81 and 3.82 of [48].In Figure 14 we show dT /dt / ( HT ) for a decoupled species ( δρ/δt = 0, and δn/δt = 0)with negligible chemical potential, µ = 0. The Fermi-Dirac and Bose-Einstein cases areobtained by numerically evaluating (2.9a), while the Maxwell-Boltzmann case is a resultof (A.27). From Figure 14 we can clearly appreciate that for T (cid:29) m , dT /dt = − HT and for T (cid:28) m , dT /dt = − HT . We therefore recover the correct behaviour [48] and also we noticethat the Fermi-Dirac and Bose-Einstein cases are relatively similar.– 40 – − − − T / m − d T / d t / ( H T ) MBFDBE
Figure 14 . − dT /dt/ ( HT ) for a massive decoupled relic as a function of T /m . We have chosen µ = 0.The different lines correspond to a particle that follows Fermi-Dirac (red), Bose-Einstein (blue) orMaxwell-Boltzmann (black) statistics. An application: Baryon temperature evolution
The evolution equations (A.30) allow us to calculate, for example, how the backgroundbaryon temperature evolves with time. This is clearly not a new result [92], but serves asan application of the formulae derived above. Since baryons are tightly coupled to electrons,and electrons are tightly coupled to photons via highly efficient elastic scatterings: δn/δt = 0and hence, we can safely neglect chemical potentials. From Equation (A.30a), we obtain thatthe temperature evolution for the electron-baryon fluid is: dT b dt = − H T b + 23 1 n b δρ bγ δt . (A.31)Then, we are simply left with calculating δρ bγ /δt . Since baryons are strongly coupled toelectrons via Rutherford scattering, we can assume that baryons and electrons share the sametemperature, T e = T b , and therefore δρ bγ /δt = δρ eγ /δt = − δρ γe /δt . We simply calculate δρ γe /δt by considering the collision term for electron-photon Compton scattering [93, 94]: C γ [ f γ ] = x e n e σ T m e p ∂∂p (cid:20) p (cid:18) T b ∂∂p f γ + f γ (1 + f γ ) (cid:19)(cid:21) , (A.32)where x e is the free electron fraction and σ T is the Thompson scattering cross section. Byconsidering f γ = (cid:2) − e p/T γ (cid:3) − and integrating over photon momenta, 2 / (2 π ) p dp , onefinds [93]: δρ bγ δt = 4 n e ρ γ σ T x e m e ( T γ − T b ) . (A.33)By plugging this expression into Equation (A.31) and by considering that n b = ρ b /µ M (where µ M is the mean molecular baryonic-electron weight), we recover the well-known baryon tem-perature evolution equation [92, 95]: dT b dt = − H T b + 83 µ M m e ρ γ ρ b x e n e σ T ( T γ − T b ) . (A.34)– 41 – .6 Thermodynamic formulae We explicit the relevant thermodynamic formulae for the number, energy and pressure densi-ties and their derivatives for particles, with g internal degrees of freedom, that follow Fermi-Dirac, Bose-Einstein and Maxwell-Boltzmann statistics. We particularize also for the case of m = 0. We remind the reader that the distribution functions are given by: f FD ( E ) = 11 + e ( E − µ ) /T , f BE ( E ) = 1 − e ( E − µ ) /T , f MB ( E ) = e − ( E − µ ) /T . (A.35)Note that the entropy density for a perfect fluid is s = ρ + p − µnT [10]. In addition, we note thatthe formulae for n , ρ , p and their derivatives can also be written as an infinite sum of Besselfunctions. This can be useful to find these quantities without integration, see e.g. [96]. Fermi-Dirac n FD = g π (cid:90) ∞ m dE E √ E − m e ( E − µ ) /T + 1 , (A.36a) ρ FD = g π (cid:90) ∞ m dE E √ E − m e ( E − µ ) /T + 1 , (A.36b) p FD = g π (cid:90) ∞ m dE ( E − m ) / e ( E − µ ) /T + 1 , (A.36c) ∂n FD ∂T = g π (cid:90) ∞ m dE E (cid:112) E − m ( E − µ )4 T cosh − (cid:18) E − µ T (cid:19) , (A.36d) ∂ρ FD ∂T = g π (cid:90) ∞ m dE E (cid:112) E − m ( E − µ )4 T cosh − (cid:18) E − µ T (cid:19) , (A.36e) ∂n FD ∂µ = g π (cid:90) ∞ m dE E (cid:112) E − m (cid:20) T cosh (cid:18) E − µT (cid:19) + 2 T (cid:21) − , (A.36f) ∂ρ FD ∂µ = g π (cid:90) ∞ m dE E (cid:112) E − m (cid:20) T cosh (cid:18) E − µT (cid:19) + 2 T (cid:21) − . (A.36g) Bose-Einstein n BE = g π (cid:90) ∞ m dE E √ E − m e ( E − µ ) /T − , (A.37a) ρ BE = g π (cid:90) ∞ m dE E √ E − m e ( E − µ ) /T − , (A.37b) p BE = g π (cid:90) ∞ m dE ( E − m ) / e ( E − µ ) /T − , (A.37c) ∂n BE ∂T = g π (cid:90) ∞ m dE E (cid:112) E − m ( E − µ )4 T sinh − (cid:18) E − µ T (cid:19) , (A.37d) ∂ρ BE ∂T = g π (cid:90) ∞ m dE E (cid:112) E − m ( E − µ )4 T sinh − (cid:18) E − µ T (cid:19) , (A.37e) ∂n BE ∂µ = g π (cid:90) ∞ m dE E (cid:112) E − m T sinh − (cid:18) E − µ T (cid:19) , (A.37f) ∂ρ BE ∂µ = g π (cid:90) ∞ m dE E (cid:112) E − m T sinh − (cid:18) E − µ T (cid:19) . (A.37g)– 42 – axwell-Boltzmann n MB = g m T e µ/T π K (cid:16) mT (cid:17) , (A.38a) ρ MB = g m T e µ/T π (cid:104) mK (cid:16) mT (cid:17) + 3 T K (cid:16) mT (cid:17)(cid:105) , (A.38b) p MB = g m T e µ/T π T K (cid:16) mT (cid:17) , (A.38c) σ MB ≡ g π (cid:90) ∞ m dE E (cid:112) E − m e − ( E − µ ) /T (A.38d)= g m T e µ/T π (cid:104) m K (cid:16) mT (cid:17) + 3 T mK (cid:16) mT (cid:17)(cid:105) , (A.38e) dn MB dT = ρ MB − µn MB T , dn MB dµ = n MB T , (A.38f) dρ MB dT = σ MB − µρ MB T , dρ MB dµ = ρ MB T . (A.38g)where K j are modified Bessel functions of the j kind. Massless case
Massless Thermodynamics , m = 0 , x ≡ e µ/T Quantity Fermi-Dirac Bose-Einstein Maxwell-Boltzmann n − g T π Li ( − x ) g T π Li ( x ) g T π xρ − g T π Li ( − x ) g T π Li ( x ) g T π xp ρ/ ρ/ ρ/ ∂n/∂T g T ( µ Li ( − x ) − T Li ( − x )) π g T (3 T Li ( x ) − µ Li ( x )) π g T (3 T − µ ) π x∂ρ/∂T g T ( µ Li ( − x ) − T Li ( − x )) π g T (4 T Li ( x ) − µ Li ( x )) π g T (4 T − µ ) π x∂n/∂µ − g T π Li ( − x ) g T π Li ( x ) g T π x∂ρ/∂µ − g T π Li ( − x ) g T π Li ( x ) g T π x (A.39)where Li j are Polylogarithms of order j . Massless case with µ = 0Massless Thermodynamics , m = 0 , µ = 0Quantity Fermi-Dirac Bose-Einstein Maxwell-Boltzmann n g ζ (3) π T g ζ (3) π T g π T ρ g π T g π T g π T p ρ/ ρ/ ρ/ ∂n/∂T n/T n/T n/T∂ρ/∂T ρ/T ρ/T ρ/T (A.40)where ζ (3) (cid:39) . .7 Number and energy density transfer rates in the MB approximation In this appendix we work out the number and energy transfer rates for decays, annihilationand scattering processes in the Maxwell-Boltzmann approximation. By numerically integrat-ing the relevant collision terms we also provide spin-statistics corrections to such rates.
Decays
The collision term for a 1 ↔ a → a particle reads [51]: C a ( p a ) = − Γ a m a m ∗ m a E a p a (cid:90) E + E − dE A ( E a , E , E a − E ) , (A.41)where Γ a is the rest frame decay width, m ∗ = m a − m + m ) + ( m − m ) /m a , E ± = m ∗ m H (cid:34) E a (cid:115) m m ∗ ± p a (cid:35) , (A.42)and A = f a [1 ± f ] [1 ± f ] − f f [1 ± f a ] , (A.43)where the + sign corresponds to bosons and the − sign corresponds to fermions.In the particular case in which the decay products are massless, the kinematical regionsimplifies to E ± = ( E a ± p a ) with m ∗ = m a . In addition, if f = f and if we neglect thestatistical terms in the collision integral and in the distribution functions we simply find C a ( p a ) = − Γ a m a E a (cid:20) e µa − EaTa − e µ − EaT (cid:21) . (A.44)If we consider that the a particle is a boson and that the 1 and 2 particles are identicalmassless fermions, we find C a ( p a ) = − Γ a m a E a Tp a (cid:16) e EaTa + µT − e EaT + µaTa (cid:17)(cid:16) e EaT − e µT (cid:17) (cid:16) e EaTa − e µaTa (cid:17) log (cid:16) e Ea − pa T + e µT (cid:17) (cid:16) e Ea +2 µ + pa T + e E a /T (cid:17)(cid:16) e Ea + pa T + e µT (cid:17) (cid:16) e Ea +2 µ − pa T + e E a /T (cid:17) (A.45)The integral over momenta, assuming Maxwell-Boltzmann statistics (A.44), yields: δn a δt = Γ a g a m a π (cid:20) T e µT K (cid:16) m a T (cid:17) − T a e µaTa K (cid:18) m a T a (cid:19)(cid:21) , (A.46) δρ a δt = Γ a g a m a π (cid:20) T e µT K (cid:16) m a T (cid:17) − T a e µaTa K (cid:18) m a T a (cid:19)(cid:21) . (A.47)From the MB equations above, and taking T a = T and µ a = µ = 0, we recover the well-knowncases: (cid:104) Γ (cid:105) a → ≡ n a δn a δt (cid:12)(cid:12)(cid:12)(cid:12) a → = Γ a K ( mT ) K ( mT ) , (A.48) (cid:104) Γ (cid:105) → a ≡ n δn a δt (cid:12)(cid:12)(cid:12)(cid:12) → a = Γ a g a g (cid:16) mT (cid:17) K (cid:16) mT (cid:17) , (A.49)– 44 – − m / T S p i n - S t a t i s t i cs C o rr e c t i on R a t e s a r e M a x i m a l δρ/δ t δ n /δ t B ↔ FF F ↔ F B Figure 15 . Statistical factor correction to the Maxwell-Boltzmann number and energy density ratesfor decays B → F F and F → F (cid:48) B . We have calculated this corrections assuming negligible chemicalpotentials and ∆ T /T = 0 .
01. We highlight the region in which the rates are maximal, T ∼ m/ which in relevant limits read (see e.g. Equations 6.8 in [48]): (cid:104) Γ (cid:105) a → (cid:12)(cid:12) T (cid:29) m = Γ a m T , (cid:104) Γ (cid:105) a → (cid:12)(cid:12) T (cid:28) m = Γ a , (A.50) (cid:104) Γ (cid:105) → a (cid:12)(cid:12) T (cid:29) m = Γ a g a g mT , (cid:104) Γ (cid:105) → a (cid:12)(cid:12) T (cid:28) m = Γ a g a g (cid:114) π (cid:16) mT (cid:17) / e − mT , (A.51)where (cid:104)(cid:105) represents thermal averaging.The number and energy density exchange rates including quantum statistics cannotbe written analytically for decay processes. In Figure 15, by numerically integrating therelevant collision terms, we show the spin-statistics corrections to the MB transfer rates. Wecan appreciate that when the rates are maximal (which is the most relevant scenario), thecorrections are within 30%. Annihilations
We calculate the number and energy density transfer rates for 2 ↔ δρδt = − (cid:88) spins (cid:90) d ˜ p d ˜ p d ˜ p d ˜ p (2 π ) δ ( p + p − p − p ) |M| E ( f f − f f ) , (A.52)while the number density transfer reads [98]: δnδt = − (cid:88) spins (cid:90) d ˜ p d ˜ p d ˜ p d ˜ p (2 π ) δ ( p + p − p − p ) |M| ( f f − f f ) , (A.53)where d ˜ p i = d p i / (2 E i (2 π ) ) and M is the amplitude for the 1 + 2 → ρ/δt for ↔ the case: T = T = T and T = T = T (cid:48) For the particular process of interest, T = T = T and T = T = T (cid:48) . Hence, within theMB approximation f f = e − ( E + E ) /T and f f = e − ( E + E ) /T (cid:48) = e − ( E + E ) /T (cid:48) as a resultof energy conservation. Since f f − f f does not depend upon p and p , the integral overthese variables can be performed, leading to [98, 99]: (cid:88) spins (cid:90) d ˜ p d ˜ p (2 π ) δ ( p + p − p − p ) | M | = 4 g g p √ s σ ( s ) . (A.54)Where g and g are the number of internal degrees of freedom of particles 1 and 2 respectively, s is the centre of mass energy squared, p = [ s − ( m + m ) ] / [ s − ( m − m ) ] / / (2 √ s ),and σ ( s ) is the usual cross section of the process (namely, that summed over final spin statesand averaged over initial spins). This leads to δρδt = − g g (cid:90) d ˜ p d ˜ p E (cid:104) e − E E T − e − E E T (cid:48) (cid:105) p √ s σ ( s ) . (A.55)By rewriting the integral in terms of the convenient variables E + ≡ E + E , E − ≡ E − E and s , the phase space density is: d ˜ p d ˜ p = 1(2 π ) dE + dE − ds , (A.56)with the appropriate limits of integration: s ≥ ( m + m ) , (A.57) E + ≥ √ s , (A.58) (cid:12)(cid:12)(cid:12)(cid:12) E − − E + m − m s (cid:12)(cid:12)(cid:12)(cid:12) ≤ p (cid:115) E − ss . (A.59)The energy transfer rate then reads: δρδt = − g g (cid:90) π ) dE + dE − ds E + + E − (cid:20) e − E + T − e − E + T (cid:48) (cid:21) p √ s σ ( s ) . (A.60)The first integral to be performed is over E − , yielding: (cid:90) E + + E − dE − = 2 E + p (cid:115) E s − (cid:18) m − m s (cid:19) . (A.61)Then, we perform the integral over E + to find: (cid:90) e − E + T dE + (cid:18)(cid:90) E + + E − dE − (cid:19) = 2 p √ s T (cid:0) s + m − m (cid:1) K (cid:18) √ sT (cid:19) , (A.62)where K is a modified Bessel function of the second kind.This allows us to write the energy transfer rate in terms of a single integral over s : δρδt = g g π (cid:90) ∞ s min ds p (cid:0) s + m − m (cid:1) σ ( s ) (cid:20) T (cid:48) K (cid:18) √ sT (cid:48) (cid:19) − T K (cid:18) √ sT (cid:19)(cid:21) , (A.63)– 46 –ates [ δρ/δt ] stats (cid:14) [ δρ/δt ] MB [ δn/δt ] stats (cid:14) [ δn/δt ] MB Initial State Final State |M | |M | |M | |M | |M | |M | F + F F + F 0.88 0.79 0.64 0.85 0.74 0.57F + F B + B 1.02 1.06 1.24 1.03 1.08 1.35B + B B + B 1.20 1.54 3.40 1.28 1.83 6.44
Table 6 . Statistical factors corrections to the annihilation rates as compared to the Maxwell-Boltzmann case for a four-point interaction and for m i = 0. M is a matrix element, |M n | ∝ p n ,where p n represents any combination of 4-momenta scalar products with dimension of energy n . Forexample, ( p · p )( p · p ) corresponds to n = 4. These spin-statistics correction factors are calculatedassuming that ( T − T ) /T = 0 .
01 and negligible chemical potentials. where s min = min[( m + m ) , ( m + m ) ]. For the particular case in which m = m = m = m = 0, the previous expression reads δρδt = g g π (cid:90) ∞ ds s σ ( s ) (cid:20) T (cid:48) K (cid:18) √ sT (cid:48) (cid:19) − T K (cid:18) √ sT (cid:19)(cid:21) . (A.64)The effect of quantum statistics in the collision depends upon the mass of the speciesinvolved in them. The larger m/T is the more they resemble the MB formulas presentedhere. Table 6 outlines the effect on the rates by assuming that all the relevant species aremassless and have a very small temperature difference of ∆ T ≡ ( T − T ) /T = 0 . An example: Neutrino-electron interactions
We can consider the case of the Standard Model neutrino-electron interactions. In that case,the neutrino-electron annihilation cross section is simply given by σ ν α ¯ ν α → e + e − ( s ) = 23 G F π ( g L + g R ) s , (A.65)and by using Equation (A.63) we find: δρδt (cid:12)(cid:12)(cid:12)(cid:12) ν α , ν α ¯ ν α → e + e − = 64 G F (cid:0) g L + g R (cid:1) π ( T γ − T ν ) , (A.66)where we have used the fact that g ν α = g ¯ ν α = 1. Equation (A.66) exactly matches EquationA.8 of [1] and in order to find the joint energy transfer for neutrinos and antineutrinos (asconsidered in Section 3) one should simply multiply this expression by 2. δn/δt for ↔ the case: T = T = T and T = T = T (cid:48) Again, by working in terms of E + , E − and s we find: δnδt = g g π (cid:90) ∞ s min ds p √ s σ ( s ) (cid:20) T (cid:48) K (cid:18) √ sT (cid:48) (cid:19) − T K (cid:18) √ sT (cid:19)(cid:21) . (A.67)Since δnδt = (cid:104) σv (cid:105) ( n − n ), this expression clearly matches the usual annihilation cross sectionformula: (cid:104) σv (cid:105) = 18 T m m K ( m /T ) K ( m /T ) (cid:90) ∞ s min ds s / K (cid:18) √ sT (cid:19) λ (cid:20) , m s , m s (cid:21) σ ( s ) , (A.68)where λ ( x, y, z ) = x + y + z − xy − xz − yz is the Kallen function.Note that in the limit in which m = m = m and T (cid:28) m , δnδt (cid:39) m δρδt for annihilationinteractions. – 47 – ρ /δt S |M| f F F stat f F B stat f BB stat68 π M ∗ T T ( T − T ) ( p · p )( p · p ) /M ∗ π M ∗ T T ( T − T ) ( p · p )( p · p ) /M ∗ π M ∗ T T ( T − T ) ( p · p )( p · p ) /M ∗ π M ∗ T T ( T − T ) ( p · p ) /M ∗ π M ∗ T T ( T − T ) ( p · p ) /M ∗ π M ∗ T T ( T − T ) ( p · p ) /M ∗ π M ∗ T T ( T − T ) ( p · p ) /M ∗ π M ∗ T T ( T − T ) ( p · p ) /M ∗ π M ∗ T T ( T − T ) ( p · p ) /M ∗ λ π T T ( T − T ) λ Table 7 . Scattering energy exchange rates for scattering processes 12 →
12 where m = m = 0 and g = 1. λ is a generic coupling constant and M (cid:63) is a mass scale. M is the matrix element of theprocess and S is the symmetry factor, see [30] for more details. The third, fourth and fifth columnscorrespond to the spin-statistics corrections to these rates for the scattering of various particles thathave negligible chemical potentials and ( T − T ) /T = 0 . Scatterings
The number density transfer rate for scattering interactions vanishes by definition since scat-tering processes conserve particle number. The energy transfer rate does not necessarilyvanish, although if annihilation or decay interactions are active then scattering interactionsrepresent a subdominant contribution to the energy transfer since the typical energy ex-changed in an scattering interaction is the momentum transferred ∼ q while the energytransferred in an annihilation or decay is ∼ E .The problem with scattering rates is that the collision term for such interactions istypically rather complicated, see e.g. [19, 74]. This is a mere consequence of the fact that forscattering interactions the collision term is proportional to e − E /T e − E /T − e − E /T e − E /T and cannot be further simplified as in the case of annihilation interactions, particularly formassive particles. Therefore, to actually calculate the exact collision term or energy transferrate, one should resort to the rather involved methods of [32] in order to reduce the phasespace to 2 dimensions (see also [55, 56]). We also point the reader to [100] where the phasespace is reduced to 1 dimension when all species are massless. We, however, warn the readerthat typos exist in that reference, see footnote 3 of [1]. Given the difficulty of the calculationand that no generic expression exists, here, by using the methods of [32] we outline the energytransfer rate for scattering processes involving massless species interacting via a four-pointinteraction. They are reported in Table 7 together with spin-statistics corrections to them.As an example of such formulae, one can consider the scattering interactions betweenneutrinos of different flavor. The matrix element for the process ν α ν β → ν α ν β where α (cid:54) = β is S |M| = 2 G F ( p p )( p p ) [30]. By using the first row of Table 7 we can relate 1 /M (cid:63) → G F . For the case of ν e and summing over µ and τ flavors, we find: δρ ν e δt (cid:12)(cid:12)(cid:12)(cid:12) ν e ν µ ↔ ν e ν µ + ν e ν τ ↔ ν e ν τ = 48 G F π T ν e T ν µ (cid:2) T ν µ − T ν e (cid:3) ,,