Towards a precision calculation of N eff in the Standard Model. Part III. Improved estimate of NLO contributions to the collision integral

We compute the dominant QED correction to the neutrino-electron interaction rate in the vicinity of neutrino decoupling in the early universe, and estimate its impact on the effective number of neutrino species N eff in cosmic microwave background anisotropy observations. We find that the correction to the interaction rate is at the sub-percent level, consistent with a recent estimate by Jackson and Laine. Relative to that work we include the electron mass in our computations, but restrict our analysis to the enhanced t-channel contributions. The fractional change in N eff SM due to the rate correction is of order 10-5 or below, i.e., about a factor of 30 smaller than that recently claimed by Cielo et al., and below the nominal computational uncertainties of the current benchmark value of N eff SM = 3.0440 ± 0.0002. We therefore conclude that aforementioned number remains to be the state-of-the-art benchmark for N eff SM in the standard model of particle physics.

The effective number of neutrinos, N eff , is an important parameter in standard hot big bang cosmology [1].Defined as the energy density residing in free-streaming, ultra-relativistic particle species relative to the photon energy density in the post-neutrino decoupling early universe (i.e., at temperature T ≲ 1 MeV), the primary role of N eff is to fix the universal expansion rate at T ≲ 1 MeV up to the end of the radiation-domination epoch.Its observable consequences are many-from setting the primordial light element abundances, to influencing the correlation statistics of the cosmic microwave background (CMB) anisotropies and the large-scale matter distribution.Coupled with the fact that many beyond-the-standard-model scenarios predict N eff -like effects (e.g., light sterile neutrinos [2,3], axions [4,5], gravitational waves [6], hidden sectors [7,8], etc.), pinning down its value both theoretically and observationally has enjoyed an unwavering interest for over four decades [1,9].From the theoretical perspective, the expected value of N eff in the context of the standard model (SM) of particle physics is 3 (for three generations), plus percent-level corrections due to residual energy transfer between the quantum electrodynamic (QED) plasma and the neutrino sector during neutrino decoupling [10][11][12][13][14][15] as well as deviations of the QED plasma itself from an ideal gas [16][17][18][19][20][21].Historical estimates of these corrections have ranged from 0.011 to 0.052 [12,[22][23][24][25]. Detailed modelling in recent years [26][27][28][29], however, have drastically reduced the spread.Notably, two of us (Drewes and Wong) and collaborators reported in [29] a prediction of N SM eff = 3.0440 ± 0.0002 from a fully momentum-dependent precision transport calculation that accounted for (i) neutrino oscillations, (ii) finite-temperature corrections to the QED equation of state (EoS) to O(e 3 ), and (iii) a first estimate of finite-temperature corrections of type (a) to the weak interaction rates depicted in figure 1 (see also table 1); the error bars are due mainly to numerical resolution and experimental uncertainties in the input neutrino mixing angles.More remarkably still, this result is in perfect agreement with the independent calculations of [27,28] modelling the same physics-to five significant digits in the central value and with comparable error estimates.The precision computation of N SM eff appears therefore to have reached convergence, at least in a limited sense.
There are nonetheless reasons to be cautious.For one, while the general expectation is that the physics summarised in table 1 dominate corrections to N SM eff , as yet missing is a systematic study of possible higher-order effects that may be at least as important.Indeed, in estimating only next-to-leading-order (NLO) effects due to diagram (a) to the weak rates rather than the full range of corrections displayed in figure 1, one could argue that the computations of [27][28][29] are, at least conceptually, incomplete.Two recent works have taken a first step towards filling this gap: • Cielo et al. [30] took the finite-temperature rate corrections for e + e − → ν α να from [31], computed originally in the context of energy loss in a stellar plasma, and appealed to detailed balance to estimate the corrections to collision integrals.The claimed effect of this correction on N SM eff is quite substantial-at the ∼ 0.001 level.We have reservations about this result: Aside from mapping rate corrections from an incompatible energy

Sources of uncertainty
Numerical solution by FortEPiaNO ±0.0001 Input solar neutrino mixing angle θ 12 ±0.0001 Table 1: Leading-digit contributions from various SM corrections, in order of importance, thus far accounted that make up the final N SM eff − 3 = 0.0440 ± 0.0002 [28,29].Table adapted from [29].regime (O(1) keV ≪ m e in a stellar plasma versus O(1) MeV > m e in the early universe), the analysis of [30] also neglected (i) corrections due to diagram (d) in figure 1 (the "closed fermion loop"), (ii) corrections to elastic scattering reactions like ν α e → ν α e, as well as (iii) Pauli blocking effects of the final-state neutrinos.Of particular note is that the neglected closed fermion loop diagram (d) contains a t-channel enhancement, which should, at least naïvely, dominate the weak rate corrections.
• The more recent work of Jackson and Laine [32] considered all diagrams in figure 1 in a first-principles calculation using the imaginary-time formalism of thermal field theory which accounts for both vacuum and finite-temperature corrections, 1 along with an estimate of hadronic corrections to diagram (d) following [34].Because the imaginary-time formalism applies only to systems in thermal equilibrium, the computation of [32] effectively neglected non-equilibrium neutrino phase space distributions.Additionally, the authors assumed a negligible electron mass by setting m e = 0, which is not necessitated by the formalism.The final result, presented as a set of corrections to the leadingorder (LO) neutrino damping rate, confirms the expected t-channel enhancement in diagram (d).Nevertheless, these corrections are minute-of order 0.2 to 0.3%.While the authors of [32] did not report the corresponding change in N SM eff , it is clear that corrections of this magnitude cannot effect a shift in N eff as sizeable as that claimed in [30].
The purpose of the present work is to clarify whether or not QED corrections to the neutrino interaction rate can alter the standard-model N SM eff at the ∼ 0.001 level as claimed in [30].While this correction is small relative to the anticipated sensitivity of the next-generation CMB-S4 program to N eff , σ(N eff ) ≃ 0.02 − 0.03 [35], an accurate theoretical prediction for N SM eff at per-mille level accuracy is nonetheless desirable to justify neglecting the theoretical uncertainty in cosmological parameter inference.In this regard, our work shares a common motivation with Jackson and Laine [32].However, our work also differs from [32] in three important ways: 1. Our first-principles computation of the QED corrections uses the so-called real-time formalism and focuses on corrections of type (d) which contain a t-channel enhancement.Computation of the non-enhanced (and hence sub-dominant) diagrams (a)-(c) is postponed to a follow-up work.Like the imaginary-time formalism, the real-time formalism also automatically takes care of both vacuum and finite-temperature corrections.While in equilibrium situations the two formalisms are exactly equivalent [36], the latter has the advantage that it can easily be generalised to non-equilibrium settings [37,38].Thus, although we shall restrict the present analysis to the same equilibrium conditions as in [32] and hence provide an independent partial cross-check of their results, our calculation will also pave the way for incorporating NLO effects in neutrino decoupling codes such as FortEPianNO [26] to deliver a final verdict on N SM eff .
2. We retain a finite electron mass m e in our calculation, which represents a departure from the m e = 0 approximation made in [32].Neutrino decoupling occurs at temperatures T ∼ 1 MeV, in whose vicinity the setting of m e = 0 may not be well justified.We shall examine the validity of the approximation and its impact on N SM eff .
3. Finally, we use our NLO results to estimate the corresponding change in N SM eff .This estimate was missing in [32].
The rest of the article is organised as follows.In section 2 we describe the physical system and sketch how N SM eff is computed.Section 3 outlines the computation of QED corrections to the neutrino damping rate due to diagram (d) and presents our numerical estimates of these rate corrections.The shift in N SM eff due to these corrections is presented in section 4, and section 5 contains our conclusions.Five appendices provide details on the bosonic and fermionic propagators at finite temperatures, explicit expressions for certain thermal integrals, derivation of the 1PI-resummed photon propagator, parameterisation of the collision integrals, and functions appearing in the continuity equation.

Preliminaries
The SM effective number of neutrinos N SM eff is defined via the ratio of the energy density contained in three generations of neutrinos and antineutrinos, ρ ν , to the photon energy density ρ γ in the limit T /m e → 0, Here, T is the photon temperature, m e the electron mass, and the limit T /m e → 0 is understood to apply only to T ≫ m i , where m i is the largest neutrino mass.
To compute the precise value of N SM eff requires that we track the evolution of ρ ν and ρ γ simultaneously across the neutrino decoupling epoch, i.e., T ∼ O(10) → O(0.01) MeV.Assuming that at these temperatures the photons are held in a state of thermodynamic equilibrium together with the electrons/positrons-collectively referred to as the "QED plasma"-this is typically achieved by solving two sets of evolution equations: (i) a continuity equation that follows the total energy density of the universe, and (ii) some form of generalised Boltzmann equations-often referred to as the quantum kinetic equations (QKEs)-which describe the non-equilibrium dynamics in the neutrino sector across its decoupling.

Continuity equation
In a Friedmann-Lemaître-Robertson-Walker (FLRW) universe, the continuity equation is given by where ρ tot and P tot are, respectively, the total energy density and pressure, H ≡ (1/a)(da/dt) the Hubble expansion rate, and a is the scale factor.For the physical system at hand, ρ tot ≡ ρ QED + ρ ν and P tot ≡ P QED + P ν , where ρ QED ≡ ρ γ + ρ e subsumes the photon and the electron/positron energy densities, and similarly for P QED .The assumption of thermodynamic equilibrium in the QED sector in the time frame of interest means that the standard thermodynamic relation ρ QED = −P QED + T (∂/∂T )P QED applies.Then, the finitetemperature corrections to the QED equation of state summarised in table 1 simply refer to corrections to the QED partition function Z QED and hence P QED = (T /V ) ln Z QED that alter ρ QED +P QED = T (∂/∂T )P QED from its ideal-gas expectation.Corrections to Z QED are known to O(e3 ) for arbitrary m e and chemical potential µ [39] and to O(e 5 ) for m e = µ = 0 [40,41].Their effects on N SM eff have been estimated in references [17,20,21] to O(e 4 ).
Quantum kinetic equations State-of-the-art neutrino decoupling calculations employ the QKEs to track simultaneously the effects of in-medium neutrino flavour oscillations and particle scatterings on the one-particle reduced density matrix of the neutrino ensemble, ϱ = ϱ(p, t).Schematically, the QKEs take the form [42] where H = H(p, t) = H vac + V is the effective Hamiltonian incorporating vacuum flavour oscillations H vac and in-medium corrections from forward scattering V, 2 and is the collision integral encapsulating the non-unitary gains (Γ < = Γ < (p, t)) and losses (Γ > = Γ > (p, t)) of ϱ.In the context of a precision calculation of N SM eff , the quantities ϱ, H, and Γ ≷ are understood to be 3 × 3 hermitian matrices in flavour space, with the diagonal entries of ϱ corresponding to the occupation numbers f α (p, t) ≡ {ϱ(p, t)} αα , for α = e, µ, τ . 3We assume a CP -symmetric universe, which is well justified if any asymmetry in the lepton sector mirrors the baryon asymmetry in the observable universe [43].In practice this means it suffices to follow only one set of QKEs for the neutrinos; the antineutrinos evolve in the same way.
For the problem at hand, the gain and loss terms Γ ≷ incorporate in principle all weak scattering processes wherein at least one neutrino appears in either the initial or final state.In the temperature window of interest, however, it suffices to consider only 2 → 2 processes involving (i) two neutrinos and two electrons any way distributed in the initial and final states (labelled "νe"), and (ii) neutrino-neutrino scattering ("νν").The leading-order Γ ≷ for these processes are well known, and take the form of two-dimensional momentum integrals, The momentum integrals (2.5) are hard-coded in the purpose-built neutrino decoupling code FortEPianNO [26,29], which solves numerically the continuity equation (2.2) and the three-flavour QKEs (2.3) in their entirety in precision N SM eff computations.

Damping approximation
We would like to compute QED corrections to Γ ≷ νe , and estimate their impact on N SM eff .Ideally these corrections would take the form of corrections to the scattering kernel Π νe , to be incorporated into a neutrino decoupling code such as FortEPianNO.As a first pass, however, we may work within the damping approximation, which makes the simplifying assumption that all particle species-besides that at the momentum mode p represented by ϱ(p)-are in thermal equilibrium specified by a common temperature T equal to the photon temperature.Then, defining δϱ = ϱ(p) − ϱ eq (p), where ϱ eq (p) = diag(f is the Fermi-Dirac distribution, the collision integral (2.4) can be expanded to linear order in δϱ to give where δ αβ is a Kronecker-δ, and are the damping coefficients composed of components of the damping rate matrix Linearisation in δϱ ensures that Γ(p, T ) is diagonal in the neutrino interaction basis, and that Γ > (p, T ) = e Ep/T Γ < (p, T ) holds by detailed balance, where E p = p is the neutrino energy at the momentum mode p of interest.This is the approximation under which we shall compute Γ in section 3. Details of the derivation of (2.6) can be found in, e.g., appendix B of [29].
In the following, we shall devote section 3 to evaluating QED corrections to the damping rate {Γ νe (p, T )} αα due to the closed fermion loop (diagram (d) in figure 1).Our estimates of its impact on N SM eff will be presented in section 4.

Calculation and NLO results
We compute QED corrections to the interaction rates of the weak processes ν α e → ν α e and ν α να ↔ e + e − within the framework of Fermi theory, whose effective Lagrangian can be expressed as Here, the spinors ψ α represent neutrino fields of flavour α; J µ L/R = ψe γ µ P L/R ψ e are the leftand right-handed electron current operators, with the chiral projectors P L/R = 1 2 (1 ∓ γ5 ); the right-handed coupling reads g R = sin 2 θ W , while the left-handed coupling g α L = − 1 2 + sin 2 θ W + δ αe depends on the neutrino flavor α = e, µ, τ .Strictly speaking, the closed fermion loop in figure 1 receives in principle contributions also from quarks.This interaction is also well described by a Lagrangian of the form (3.1), with the couplings g α L and g R updated for the quarks of interest. 4We shall however omit the contributions from quarks in this work: at finite temperatures these contributions are Boltzmann-suppressed at the O(1) MeV temperatures of interest, while hadronic corrections in vacuum have been studied in [34].As QED is a vector-like theory, we also introduce the vector-axial couplings g α V,A = 1 2 (g α L ± g R ) as an alternative notation.

Evaluation of the closed fermion loop
In order to compute the damping coefficients (2.7), we must evaluate the flavour-diagonal entries of the neutrino damping rate matrix in the equilibrium approximation for ϱ, i.e., Γ α (p, T ) ≡ {Γ νe (p, T )} αα , which is given by the sum where Γ ≷ α are the α-flavoured neutrino production (<) and destruction (>) rates in the equilibrium limit.In the framework of the real-time formalism of finite-temperature QED, these rates can be extracted from the Wightman self-energies Σ ≷ , 5 where Σ ≷ are identified with the self-energies with opposite (Keldysh) contour indices:6 Σ < ≡ Σ 12 and Σ > ≡ Σ 21 .In thermal equilibrium, detailed balance is established by the fermionic Kubo-Martin-Schwinger (KMS) relation Σ > = −e p 0 /T Σ < .Then, the modedependent damping rate can also be written as where the first equality follows from the fact that Γ α can be related to the discontinuity of the retarded neutrino self-energy evaluated at the quasiparticle pole, −idiscΣ R = 2 ImΣ R = Σ > − Σ < , in accordance with the optical theorem at finite temperature.While the KMS relation makes explicit use of the fact that Σ ≷ are computed in thermal equilibrium, the realtime formalism used here can be generalised in a straightforward manner to non-equilibrium situations [37]. 7Where there is no confusion, we shall drop the flavour index α in the following.On general grounds the dominant QED correction to Γ α can be expected from diagram (d) in figure 1 in the regime where the photon propagator is on-shell.This expectation has been confirmed in [32].In the real-time formalism this process is contained in the contribution to the Wightman self-energy shown in figure 2, where, in terms of finite-temperature cutting rules, diagram (d) can be recovered from a cut through one closed electron loop and the internal neutrino line. 8We shall compute only this contribution and denote the corresponding Wightman self-energy with Σ ba NLO .For notational convenience, we further split Σ ba NLO into a sum Σ ba NLO = c,d=1,2 Σ ba,cd over the internal real-time indices c and d, with the partial self-energies Σ ba,cd given by Tr iS ad e (q)γ ρ (P Here, we have introduced the abbreviation p ≡ d 4 p /(2π) 4 for loop integrals; the traces are over the Clifford algebra; the subscripts "e" and "ν α " on the fermionic propagators refer to the associated particle; and the definitions of the propagators can be found in appendix A.
As we would like to compute Σ 12 NLO , we set the external contour indices to a = 2 and b = 1.Then, the contributions Tr / pΣ 12,21 and Tr / pΣ 12,12 , which correspond to cutting both the photon and neutrino lines, vanish by momentum conservation (see also footnote 8).Of the remaining "11" and "22" contributions, the transformation behaviour of the thermal propagators under hermitian conjugation dictates that Tr / pΣ 12,11 = Tr / pΣ i.e., the Wightman self-energy represented by figure 2 can be determined entirely through the diagonal "11" contribution.
It is instructive to recast the expressions for the self-energy (3.5)-(3.6)into the form of a standard Boltzmann collision integral commonly found in textbooks (e.g., [50]).To do so, we first identify the 4-momenta p, l, q in figure 2 with the momenta p 1 , p 2 , p 3 , p 4 of the external neutrinos and electrons of the underlying 2 → 2 process represented by diagram (d) of figure 1 via Then, writing out the propagators iS where the population factor contains the equilibrium phase space distributions of the three integrated external particles and is analogous to F ≷ in equation (2.5).The quantity T NLO plays the role of QED corrections to the LO squared matrix element, 9 and can be written in terms of the one-loop photon selfenergy Π µν ab as where we have introduced the photon momentum Since T NLO has the interpretation of a squared matrix element, we can split it into a vacuum and a thermal part according to temperature dependence, The vacuum correction T vac has no intrinsic temperature dependence in the sense that it makes no explicit reference to the temperature or to any phase space distribution.It is simply the correction to the weak matrix elements arising from the interference of the closed ) .If we were to replace TNLO with TLO in equation (3.8), we would recover the leading-order equilibrium neutrino production rates.
fermion loop diagram (d) with the LO graph in standard T = 0 quantum field theory, and can be expressed in terms of the vacuum photon self-energy as where α em is the electromagnetic fine-structure constant, and the form factor Π 2 is defined in appendix C; see equation (C.7).The simplicity of the expression follows from the fact that the integration domain is symmetric under the interchange p 3 ↔ p 4 .This symmetry, along with momentum conservation, also ensures the absence of all antisymmetric terms containing Levi-Civita symbols arising from traces of γ 5 with four or more γ-matrices.Vacuum QED corrections to the neutrino-electron interaction matrix element were previously computed in reference [33].
The thermal correction T th , on the other hand, depends explicitly on equilibrium phase space distributions, f F or f B , of the internal particles ("F " for Fermi-Dirac; "B" for Bose-Einstein), and can be thought of as a temperature-dependent correction to the squared matrix element.This T -dependence originates in the thermal part of the "11" propagator, which is proportional to δ(p 2 − m 2 )f F/B (|p 0 |) (see appendix A) and, where it is applied, effectively puts an internal line of the closed fermion loop diagram (d) in figure 1 on-shell.Purely from counting, there are altogether seven possible ways to put one, two, or all three internal lines of diagram (d) on-shell.Not all combinations contribute to T NLO : Terms proportional to f B correspond to putting the photon line on-shell, which are forbidden for diagram (d) for kinematic reasons [49,51].Similarly, putting both internal electrons on-shell leads to a purely imaginary contribution that is irrelevant to the real part of the self-energy we wish to compute.The only surviving two combinations correspond to putting either internal electron line on-shell, and are proportional to f F (|k 0 |) and f F (|k 0 + P 0 |) respectively.
Observe that T NLO contains a t-channel contribution from elastic ν α e scattering (i.e., p 0 2 < 0) that is logarithmically divergent for soft photon momenta.This divergence comes from the fact that the finite-temperature photon self-energy scales not as P 2 like in vacuum, but as T 2 in the hard-thermal loop limit which do not compensate anymore for the 1/P 2 behaviour of the photon propagator.In addition, soft photons are Bose-enhanced.To remedy the problem, we replace the tree-level photon propagator D cd µν in equation (3.10) with the fullyresummed photon propagator Dab µν . 10Furthermore, because both Dab µν and Π νσ 11 split into a longitudinal ("L") and a transverse ("T ") part, the same decomposition applies also to T th , i.e., T th = T L th + T T th , where T L,T th can be brought into the form Here, DR denotes the retarded resummed photon propagator; Π R,T ̸ =0 is the retarded thermal photon self-energy comprising the f F (|k 0 |) and f F (|k 0 + P 0 |) terms described above; we have employed the shorthand notation P i,j L/T = P µν L/T p i,µ p j,ν ; and a L,T = (3 ∓ 1)/4 further differentiates between the longitudinal and the transverse contribution.
Note that the imaginary part of the photon propagator, Im DL,T 11 , does not appear in equation (3.13) because it is formally of higher-order in α em and we are only interested in the O(α em ) corrections.We also only use the resummed photon propagator in the IR divergent tchannel contribution; where the divergence is absent, i.e., in the s-channel, we set DR → D R , where D R is the un-resummed counterpart of DR .The final expressions for Re Π 11 .Details of their derivation, along with the relevant thermal integrals, can be found in appendices B and C.

Numerical results for NLO neutrino interaction rate
We evaluate the self-energy correction (3.8) and hence the correction to the neutrino damping rate (3.4) by numerically integrating over the two remaining momenta l and q in (3.8) using the parameterisation shown in appendix D. We adopt the experimentally-determined values given by the Particle Data Group [52] for the following input parameters: • Electron mass: m e = 0.510 998 950 00(15) MeV, • Electromagnetic fine-structure constant: α −1 em (0) = 137.035999 180 (10), and • Weinberg angle: sin 2 θ W (0) MS = 0.238 63 (5).
The renormalisation scale µ R appearing in the photon self-energy of the vacuum contribution is identified with the electron mass, µ R = m e .Figure 3 shows the closed fermion loop corrections to the damping rates Γ e (p) and Γ µ,τ (p) at the mean momentum p = 3.15T .Relative to their respective LO rates, we find the corrections at temperatures T ∼ 1 → 3 MeV to fall in the range −0.2 → +0.1% and −0.0005 → +0.0002%, respectively, for ν e and ν µ,τ .We further note that: 1.At T ∼ 2 MeV, the vacuum and the thermal contributions to the correction are roughly equal in magnitude (green vs blue lines in figure 3), in contrast to the findings of [30], where finite-temperature corrections were determined to be subdominant.We note however that a direct comparison is not possible because a different set of diagrams was investigated in [30]-(a), (b), and (c) in figure 1-as opposed to the type (d) corrections considered in this work.
2. Reference [30] also found no significant flavour dependence in the rate corrections: their Γ e and Γ µ,τ corrections differ by less than 1% in the temperature regime around neutrino decoupling.Our corrections are on the other hand strongly flavour-dependent-by more than two orders of magnitude-and trace their origin to the fact that electron neutrinos experience charge current interactions while the muon-and tau-flavoured neutrinos interact only via the neutral current with the e ± thermal bath.This difference renders the corresponding vector couplings, g e V ∼ 0.49 and g µ,τ V ∼ −0.012, very roughly two orders of magnitude apart from one another.This strong flavour dependence in the rate corrections was also seen in reference [32].
3. We have computed the NLO contributions to the interaction rates in two different approximations: (i) retaining the full dependence on the electron mass (solid lines in figure 3), and (ii) in the limit m e → 0 (dashed lines).The massless calculation aims to quantify the effect of the m e = 0 approximation used in reference [32], along with the Hard Thermal Loop (HTL) approximation of the photon propagator.We observe that the error from neglecting m e is relatively minor for T ≳ 3m e , but becomes sizeable at low temperatures.In particular, in the limit T → 0 the ra-tio Γ NLO α /Γ LO α vanishes for m e > 0, but diverges for m e = 0 because of Boltzmann suppression of the LO rates.Precision computations of N SM eff track the evolution of neutrinos down to temperatures much below m e [28,29].Thus, although it is commonly understood that (electron) neutrino decoupling occurs at relativistic temperatures T ∼ (2 → 3) × m e , what we assume for m e in the NLO rate computations may yet have some impact on N eff (see section 4.2).
Figure 4 focuses on the t-channel contribution to the interaction rate, where the enhancement near the photon mass-shell occurs.We show four sets of results, based on the resummed photon propagator detailed in appendix C: (i) the complete one-loop result including a finite m e everywhere, (ii) the complete one-loop result in the limit m e → 0, (iii) using the HTL photon propagator (which does not depend on the electron mass), but with m e everywhere else, 11 and (iv) using the HTL photon propagator and setting m e = 0 everywhere.As expected, we observe that (ii) and (iv) match to a very good approximation.Indeed, since the scattering rates are dominated by the kinematic region around the t-channel singularity where photons are soft, the 1PI-resummed propagator is well-approximated by the HTL one when, in addition, we set m e = 0. On the other hand, visible differences can be discerned between (i) and (iii) at T ≲ 1 MeV, which can be explained by the fact that the HTL approximation only holds for T ≫ m e .In the lower panel of figure 4, we highlight the impact of m e by displaying the ratio of (i) to (iv).

NLO decoupling temperatures
The lower panels of the two plots in figure 3 show the LO electron neutrino and muon/tau neutrino interaction rates juxtaposed with the Hubble expansion rate.The latter is given by Pl ), where M Pl = 2.435 36 × 10 21 MeV is the reduced the Planck mass, ρ ν = 3g ν (7/8)(π 2 /30)T 4 is the energy density in three families of neutrinos with g ν = 2, and corresponding to an increase of ∼ 0.04 % for ν e and of ∼ 8 × 10 −7 % for ν µ,τ .Around the muon neutrino decoupling temperature, the vacuum and thermal contributions to the NLO rates approximately cancel, explaining the smallness of the correction to T NLO d(µ,τ ) .Given that [32] computed the NLO weak rates assuming m e = 0, it is also of interest to study how such an assumption modifies the decoupling temperature shifts.Using our

NLO effects on N SM eff
Having computed the closed fermion loop correction to the damping rate Γ α , we are now in a position to estimate its effect on N SM eff .Within the damping approximation, two avenues are available to us: • Given Γ α at a representative momentum, we have already estimated in equation (3.16) the corresponding correction to the neutrino decoupling temperature T d , defined via Γ α (T d ) = H(T d ).Then, under the assumption of instantaneous neutrino decoupling, an estimate of the change to N SM eff , δN eff ≡ N NLO eff − N LO eff , can be obtained via entropy conservation arguments.
• We may also compute δN eff by solving directly the continuity equation (2.2) and the QKEs (2.3) in the damping approximation (2.6) and neglecting neutrino oscillations.
We consider both approaches in the following.The corresponding estimates for δN eff are summarised in table 2.

Entropy conservation
The entropy conservation argument posits that entropies in a comoving volume in two decoupled sectors are separately conserved, i.e., where s να and s ≡ s γ + s e + β̸ =α s ν β + • • • denote, respectively the entropy density of a decoupled neutrino species ν α and of the QED plasma plus any neutrino species ν β that may still be in equilibrium with it, and the scale factors a 1,2 represent two different times after decoupling.We take a 1 to correspond to the time immediately after ν α decouples instantaneously from the QED plasma, i.e., a 1 = a + d , where the sector temperatures satisfy T (a 1 ) = T να (a 1 ) ≡ T d , while a 2 > a 1 is some later time to be specified below.
The assumption of instantaneous decoupling allows the neutrinos to maintain to an excellent approximation an ideal-gas equilibrium phase space distribution at all times, so that s να (a) ∝ T 3 να (a).It then follows simply from (4.1) that T να (a 2 ) = (a 1 /a 2 )T να (a 1 ) = (a 1 /a 2 )T d .In the case of the QED+ν β plasma, we parameterise its entropy density as where the entropy degree of freedom parameter gs is given in the ideal gas limit by with g γ = 2, g e = 4, and g ν β = 2.For our current purpose of estimating δN eff due to NLO contributions to the weak rates, it suffices to use the ideal-gas gs .We do note however that QED entropy density at MeV is subject in principle to a sizeable finitetemperature correction to the QED equation of state, which needs to be included in precision N SM eff calculations.See, e.g., reference [21] for details.Then, combining the above, we find an estimate of the neutrino-to-photon temperature ratio of at the later time a 2 .We use this temperature ratio in the following to provide two estimates of δN eff due to the rate corrections.
Common decoupling temperature In the first estimate, we assume all neutrino flavours to decouple effectively at the same time, an assumption that may to some extent be justified by the observed large mixing in the neutrino sector.We choose a 1 = a + d to correspond to the time immediately after ν e decoupling, and set a 2 at a time significantly after e ± annihilation, where T (a 2 ) ≪ m e .The latter leads immediately to an entropy degree of freedom of gs (a 2 ) = g γ = 2, while for the former we have Then, using the temperature ratio (4.4) and the ideal-gas relations ρ γ ∝ g γ T 4 and ρ να ∝ (7/8)g να T 4 να , we find or, equivalently in terms of N eff , where the factor 4/11 corresponds to 2/g s (T d ) evaluated in the limit T d ≫ m e .Using the LO and NLO electron neutrino decoupling temperatures in (3.15) and (3.16) respectively, we find a fractional shift of δN eff /N LO eff ≃ −1.1×10 −5 due to the rate corrections.Had we instead used the m e = 0 NLO decoupling temperatures (3.17), the correspond shift in N eff would have been δN me=0 eff /N LO eff ≃ −1.2 × 10 −5 .Thus, while setting m e = 0 leads to a ∼ 10% change in the estimate of δN eff , its ultimate impact on N SM eff appears to be small, at least within the entropy conservation picture.These results are summarised in table 2.
Flavour-dependent decoupling In the absence of neutrino oscillations, ν e and ν µ,τ decouple at different temperatures, T d(e) and T d(µ,τ ) > T d(e) .Then, to estimate δN eff requires that we consider entropy conservation across four epochs: the time immediately after ν µ,τ decoupling a 1 = a + d(µ,τ ) ; immediately before ν e decoupling a 2 = a − d(e) ; immediately after ν e decoupling a 3 = a + d(e) ; and a 4 is a time significantly after e ± annihilation.The corresponding entropy degrees of freedom are  2, we note that the flavour-dependent decoupling estimates of δN eff /N LO eff are generally about a factor of two smaller, in both the m e = 0 and m e ̸ = 0 cases.This difference is to be expected, as the QED corrections to the ν µ,τ interaction rates are negligible compared with the corrections to the ν e rates.We emphasise, however, both estimates are very crude approximations: the true shift in δN eff will probably fall somewhere in-between.

Solving the neutrino Boltzmann equations and the continuity equation
Following reference [20], we introduce the comoving quantities x ≡ m e a, y ≡ a p, and z ≡ T a, and rewrite the continuity equation (2.2) as a differential equation for the quantity z (corresponding to the photon-to-neutrino temperature ratio): where J(x/z) and Y (x/z) describe the ideal-gas behaviour of the QED plasma, while the G 1,2 (τ ) account for deviations of the QED plasma from an idea gas.Explicit forms for these expression to O(e 2 ) can be found in appendix E. Equation (4.12) also requires as input the total time derivative of the comoving neutrino density ρν ≡ ρ ν a 4 , which can generally be constructed from the neutrino density matrices via where d{ϱ(x, y)} αα /dx corresponds in general to the QKEs (2.3).Neglecting neutrino oscillations, the density matrices ϱ are diagonal in the flavour basis, in which case the QKEs (2.3) simplify to a set of classical Boltzmann equations.Then, together with the damping approximation (2.6) and in terms of the new variables, we find Equation (4.14) can be solved together with the continuity equation (4.12) for a range of momenta y covering the bulk of the neutrino population (a typical range would be y ∈ [0.01, 30]).We call this the "full-momentum" approach.Alternatively, we can simplify the equation (4.14) further by adopting the ansatz and Γ α (y) = Γ α (⟨y⟩).Identifying ⟨y⟩ with the mean momentum mode y 0 = 3.15 z(T 0 ), where T 0 is the photon temperature at initialisation (which is equal to the neutrino temperature), we can rewrite equation (4.13) in this approximation as We call this alternative the "mean-momentum" approach.
Irrespective of whether we use the full-momentum or the mean-momentum approach, the final N eff can be estimated from the solutions to ρ α and z at x → ∞ using the definition of N SM eff in equation (2.1), reproduced here in terms of the rescaled variables: Table 2 shows our estimates of δN eff /N LO eff due to QED corrections to the neutrino interaction rates using both the full-momentum and the mean-momentum approaches, with or without the electron mass in the correction.
Evidently from table 2, both the full-momentum and the mean-momentum estimates for m e ̸ = 0, δN eff /N LO eff ≃ (−7.8 → −7.9)×10 −6 , are broadly consistent with their counterpart obtained in section 4.1 from entropy arguments: in fact they sit between the common-decoupling and flavour-dependent decoupling estimates from entropy considerations.In contrast, the estimates assuming m e = 0 differ by ∼ 50% between the full-momentum (δN me=0 eff /N LO eff ≃ −2.6 × 10 −5 ) and the mean-momentum (δN me=0 eff /N LO eff ≃ −1.6 × 10 −5 ) approaches, and are furthermore ∼ 30% to a factor of four larger than those from entropy arguments.This result is consistent with expectations.As we have seen previously in figures 3 and 4, the rate correction under the m e = 0 assumption diverges as T → 0 relative the LO rate, whereas its m e ̸ = 0 counterpart tends to zero.Since neutrino decoupling in the early universe is not instantaneous but extends into the e ± -annihilation epoch at T ∼ m e , any estimate of δN eff that accounts for non-instantaneous decoupling will to some extent be sensitive to exactly what we assume for m e in the T ≲ 3m e region.Indeed, the low-T effect of m e on δN eff /N LO eff is not captured by entropy conservation arguments, which are based on one point estimate (the decoupling temperature) at T > m e , where the m e ̸ = 0 and m e = 0 rate corrections differ by less than 10%.In contrast, the full-momentum approach, which is sensitive to the largest range of temperatures, yields the strongest dependence of δN eff /N LO  eff on what we assume for the electron mass.Thus, while the overall effect of the rate corrections on N SM eff is small, we conclude that neglecting m e in their computations is strictly not a good approximation in high-precision calculations of N SM eff .

Conclusions
In this work we have computed the QED corrections to the neutrino-electron interaction rate in the vicinity of neutrino decoupling and evaluated its impact on the standard-model value of the effective number of neutrinos N SM eff .We have focused on diagram (d) in figure 1, because of the expectation of a t-channel enhancement and hence dominance over the other three diagrams.Similar corrections have also been recently investigated in Cielo et al. [30] and Jackson and Laine [32], which the former analysis [30] found to lead to a significant shift in N SM eff : from the benchmark value of 3.044 [28,29] to 3.043.Contrary [30], our first-principles calculations show that QED corrections to the neutrino interaction rates are modest.In the temperature range T ∼ 1 → 3 MeV, the corrections to electron neutrino interaction rate fall in the range −0.2 → +0.1% relative to the LO rate, while for the muon and tau neutrinos the effect is even more minute: −0.0005 → +0.0002%.These results are consistent with those reported by Jackson and Laine [32], despite differences between the formalisms (imaginary vs real-time) and some approximations (zero vs finite m e ) used.The more than two orders of magnitude difference between the relative corrections for the ν e and the ν µ,τ rates also confirms the strong flavour dependence found in [32], which was not observed in Cielo et al. [30].
Using our QED-corrected neutrino interaction rates, we proceed to estimate the corresponding shift in N SM eff under a variety of approximations and methods: via entropy conservation arguments which assume instantaneous decoupling, and by solving the Boltzmann equation in the damping approximation.Depending on the exact method/approximation used, we find the relative change in N SM eff to fall in the range δN eff /N LO eff ≃ (−0.5 → −1.1)×10 −5 .That is, relative to the current SM benchmark of N SM eff = 3.0440 ± 0.0002 [28,29], QED corrections to the neutrino interaction rates can at best shift the number in the negative direction in the fifth decimal place, and are hence completely within the quoted uncertainties.Thus, while we can confirm the sign of δN eff computed by Cielo et al. [30], even our most "optimistic" estimate is a factor of 30 smaller than their claimed correction.
It is also interesting to observe that setting m e = 0 in the rate corrections can have an O(1) impact on δN eff /N LO eff , even though the rate corrections themselves at T ∼ 1 → 3 MeV differ by less than 10%.This is because corrections assuming m e = 0 deviate from their m e ̸ = 0 counterparts at T ≲ 3m e , and diverge relative to LO rates just as the m e ̸ = 0 corrections vanish in the T → 0 limit.Since neutrino decoupling in the early universe is not instantaneous, these T ≲ 3m e effects will imprint on N SM eff despite the common understanding that neutrino decoupling occurs at T ∼ 1 MeV.Thus, while the overall contribution from QED weak rate corrections to N SM eff is ultimately small, we argue that neglecting m e their computations is strictly speaking not a good approximation in high-precision N SM eff calculations.In conclusion, our results strongly suggest that the SM benchmark value N SM eff = 3.0440± 0.0002 [28,29] is correct within the quoted uncertainties.Naturally, a full numerical solution of the QKEs by a dedicated neutrino decoupling code such as FortEPiaNO [26]-with all NLO contributions from diagrams (a) to (d) incorporated in the collision integral-would be highly desirable.However, short of a new and previously unaccounted-for effect, we believe it is unlikely that a more detailed investigation of NLO effects on the neutrino interaction rate will yield a deviation from the existing SM benchmark N SM eff that is large enough to be of any relevance for cosmological observations in the foreseeable future.while the ones for bosons are where f B (E) = 1/(e E/T − 1) is the Bose-Einstein distribution.
For the computation of the resummed photon propagator, it is useful to write down the advanced and retarded bosonic propagators, as well as the statistical propagator where the relation to the spectral function ) is another way of stating the KMS relation.In terms of the scalar propagators, the tree-level photon propagator is given by where we have defined the polarisation tensor, for a general R ξ -gauge.

B Thermal integrals
Throughout the calculations we encounter thermal integrals of the following forms To evaluate them, we choose P as the reference direction with the corresponding unit vector P, and parameterise the integration variable via where, together, P and the two unit vectors e 1 and e 2 form a complete orthogonal basis.
The contributions from e 1 and e 2 to the spatial components of K µ vanish after the azimuthal integration, giving where ω = |k 0 | = |k| 2 + m 2 e , and arise from integration over cos ϑ.Correspondingly, is the remaining, temporal component of K µ .For K µν , after application of the completeness relation δ ij = P i P j + k=1,2 e i k e j k , the spatial components become where we have introduced the two auxiliary functions The remaining components with at least one index being zero are dω ω 2P 0 ωℓ 1 (ω, P ) + P 2 ℓ 2 (ω, P ) f F (ω).

(B.8)
In all cases, the remaining single integral over ω can, for our purposes, be evaluated numerically efficiently.

C 1PI-resummed photon propagator
A crucial ingredient in our calculation is the resummed photon propagator used in our NLO thermal matrix element (3.13).In this appendix, we provide an exact derivation of the photon self-energy at the one-loop level within the real-time formalism, assuming electrons in thermal equilibrium with zero chemical potentials.From there, we extract the resummed Feynman propagator.While exact expressions for this resummed photon propagator have been long known in the simplified scenarios wherein m e = 0 [53] or in the HTL approximation (see, e.g., [54,55]), here we remove these assumptions and derive a propagator valid for all values of the electron mass m e and photon 4-momentum P .We note that, while the reference [49] also studied the photon propagator without HTL approximations (and accounting for nonzero electron chemical potentials), no explicit formulae were provided.with the corresponding projectors [36,56] P 00 T = P 0i T = P i0 T = 0, P ij T = −δ ij + P i P j , P µν L = g µν − P µ P ν P 2 − P µν T .

(C.3)
Writing the self-energy in this form makes also manifest that Π µν ab fulfills the Ward-Takahashi identity P µ Π µν ab = 0.The longitudinal projector can alternatively be expressed through the heat bath four-velocity U µ = δ µ 0 (in its rest frame) as P µν L = Ũ µ Ũ ν / Ũ 2 with Ũ µ = P 2 U µ − (U • P )P µ , so that the transverse projector becomes P µν T = g µν − P µ P ν P 2 − P µν L .In the following, we only highlight the derivation of Π 11 and Π 12 , as the remaining self-energies can be related to the first two via the bosonic KMS relation and the finite-temperature part Π µν 11,T ̸ =0 (Π µν 22,T ̸ =0 ).The renormalised vacuum part Π µν 2 (P ) = (α em /π)(P 2 g µν − P ν P µ )Π 2 (P 2 ) is textbook material and, in the MS-scheme, reads for m e → 0, with τ = 4m 2 e /P 2 (see, e.g., appendix A of [57] for the derivation).In practice, we note that the vacuum contribution to the photon self-energy is numerically irrelevant when resumming the photon propagator, as it vanishes in the soft photon exchange limit for which P 2 = 0. We can therefore safely neglect it in the resummed propagator.
In terms of the integrals K µ and K µν computed in appendix B, the corresponding real parts of the two diagonal components read After collecting all the terms according to the projectors defined in equation (C.3), we obtain for the longitudinal part with the thermal photon mass m γ = eT /3.Here, the symbol HTL ≈ refers to the result obtained in the HTL approximation [58], defined through the assumption that the quantities in the integral can be separated into hard O(πT ) and soft scales O(eT ), and assuming that hard momenta k dominate the loop integral.The external momentum P on the other hand is defined to be soft and therefore also negligible compared to the loop momentum.More precisely, to go from the fully momentum-and mass-dependent expression (C.9) (or (C.11)) to the HTL result (C.10) (or (C.12)), we have expanded in P 0 /|k| at the first step and then in P 0 /|P| in the second step, using also in the same expansion scheme.We further assume that m e ≪ T , such that the electron can be treated as effectively massless.This leaves us with an integral over ω that can be performed analytically using the relation , thereby yielding the HTL result (C.10) (or (C.12)).
The (purely imaginary) off-diagonal component Π µν 12 , on the other hand, reads where the minus sign from the fermionic loop and the one from the fact that we have a type 1 and a type 2 vertex cancel out, and For the evaluation of the integrals, we make use again of the decomposition of k shown in equation (B.2).The additional δ-function here (compared to the diagonal case) fixes then the polar angle θ to the value With that, T µν separates into a transverse and a longitudinal part according to The final result then reads (C.18) where the sum runs over positive and negative energies k 0 = ±ω.In the HTL limit we find where we have used the integral As will be explained in the subsequent section, it is convenient to compute the retarded and advanced photon self-energy Π R/A in order to perform the resummation of the photon propagator.Using the relation Im Π 11 = i 2 (Π 12 + Π 21 ), we can write the advanced and retarded self-energies as our results are in agreement with [54] after carefully sending the time-ordering parameter ϵ to zero.In the longitudinal part we differ by a term P 2 0 /|P| 2 with respect to [54], 13 but agree with [55].We note in passing that the self-energy (C.18) leads to the unphysical process of photon decay γ → e + e − at high enough temperatures, where m γ exceeds 2m e [49].In practice, this could be resolved by resumming the electron propagator.In our case, however, the relevant dynamics happen at temperatures much below that threshold, such that this resummation is not necessary.
In figure 6, we display the finite-temperature contribution to the real and the imaginary parts of the retarded transverse and longitudinal propagators, for T = 2m e and various choices of |P|/T .We compare the exact one-loop photon self-energy to (i) the equivalent selfenergy in the limit m e → 0, and (ii) the equivalent self-energy in the HTL limit.Similarly, in figure 7, we examine the impact of the temperature on the photon self-energy.We observe that the effect of the finite electron mass m e , while small for large temperatures T /m e ∼ 5, remains sizeable for temperatures around the electron neutrino decoupling temperature where T /m e ∼ 2.
Beyond the photon self-energy, other quantities of interest are the residue of the trans- 14 Longitudinal photons are also sometimes referred to as "plasmons" in the literature.It is understood that (C.34) should be replaced by the causality-respecting ϵ-prescription if the imaginary part of the photon self-energy vanishes for kinematic reasons at the given order in perturbation theory.

E Functions in the continuity equation
The continuity (4.12), written in terms of the rescaled variables x, y, z, contains the functions J(x/z), Y (x/z), and G 1,2 (x/z).We give the explicit forms of these functions here to O(e 2 ) [20], , (E.1)

2 5 19 A
Solving the neutrino Boltzmann equations and the continuity equation 18 Conclusions Bosonic and fermionic propagators at finite temperatures 20

Figure 1 :
Figure 1: The four qualitatively different QED corrections to the neutrino interaction rate corresponding to (a) modification of the electron dispersion relation, (b) virtual photon exchange, (c) thermal photon emission and absorption, and (d) corrections with a closed fermion loop.All diagrams are schematic in that all time directions are possible.
Standard-model corrections to N SM eff Leading-digit contribution m e /T d correction +0.04 O(e 2 ) FTQED correction to the QED EoS +0.01 Non-instantaneous decoupling+spectral distortion −0.006 O(e 3 ) FTQED correction to the QED EoS −0.001 Flavour oscillations +0.0005 Type (a) FTQED corrections to the weak rates ≲ 10 −4

Figure 2 :
Figure 2: The three-loop diagram containing the closed fermion loop.Explicitly labelled are the external (a and b) and summed (c and d) real-time contour indices, the momenta (k, l, p and q) and particles (e and ν α ).
self-energy (3.5)-(3.6)can be brought into the form ̸ =0 and Re DL/T 11 are given in equations (C.10), (C.12) and (C.32), which can be easily mapped to Re Π L/T R,T ̸ =0 and Re DR via Re DL/T

Figure 3 :
Figure3: Top: NLO contributions to the ν e interaction rate from the closed fermion loop with a finite electron mass (solid) or m e = 0 (dashed) at different temperatures for the mean neutrino momentum p = 3.15T .For comparison we normalise all curves to the LO rate (without QED corrections), which we always evaluate with m e ̸ = 0 to ensure a common normalisation.The total rate correction (red) is further split into the vacuum (blue) and the thermal contribution (green).At T ≫ m e , the green curve flattens out as the thermal correction contains no scale besides the temperature in this region, while the vacuum correction retains a mild dependence on the renormalisation scale µ R .The lower panel shows the LO neutrino interaction rate compared to the Hubble rate, and T d indicates the decoupling temperature, defined via Γ α (T d(e) ) = H(T d(e) ).Bottom: Same as top, but for ν µ,τ .

. 14 )αRatioFigure 4 :
Figure4: Top: The t-channel contribution (l 0 > 0) to the NLO electron neutrino interaction rate at the mean momentum p = 3.15T in different approximations, namely, using 1PIresummed (green) and HTL (magenta) photon propagators, in each case with m e ̸ = 0 (solid lines) or m e = 0 (dashed lines) in both electron loops of the self-energy diagram of figure2.As in figure3, we normalise all curves to the LO rate, which we always evaluate with m e ̸ = 0.The lower panel shows the ratio of the 1PI result for m e ̸ = 0 to the HTL result for m e = 0 as a function of the temperature.Bottom: Same as top, but for α ̸ = e.

Figure 5 :
Figure 5: One-loop contribution to the photon self-energy.

Figure 6 :
Figure 6: Comparison between the finite-temperature part of the photon self-energy at the one-loop level with (continuous lines) and without (dashed lines) the HTL approximation at 2m e and for various choices of |P|/T .

Table 2 :
Estimates of the relative correction to N SM eff due to NLO weak rate corrections, with and without the electron mass, using different methods.computations of the m e = 0 rate corrections (but with m e ̸ = 0 LO rates), we find the corresponding NLO decoupling temperatures to be i.e., a 0.05 % and 6 × 10 −6 % shift for ν e and ν µ,τ , respectively, relative to their corresponding LO decoupling temperature(3.15).
shift in N eff of δN eff /N LO eff ≃ −5.4 × 10 −6 if we use the m e ̸ = 0 corrections, or δN me=0 F (E e , T d(e) ),(4.8)andgs (a 4 ) = 2.An estimate of the ν e -to-photon energy density ratio at a 4 follows straightforwardly from the temperature ratio (4.4) and ideal-gas temperature-energy relations: For ρ νµ,τ (a 4 )/ρ γ (a 4 ), we note that ν e decoupling at a d(e) introduces a discontinuity in gs , thereby leading to a more complicated energy density ratio at a 4 ,