# Towards a precision calculation of $N_{\rm eff}$ in the Standard Model III:
Improved estimate of NLO corrections to the collision integral

Marco Drewes  Yannis Georis  Michael Klasen  Luca Paolo Wiggering  and Yvonne
Y. Y. Wong

###### Abstract

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_{\rm 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_{\rm eff}^{\rm
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_{\rm
eff}^{\rm SM}=3.0440\pm 0.0002$. We therefore conclude that aforementioned
number remains to be the state-of-the-art benchmark for $N_{\rm eff}^{\rm SM}$
in the standard model of particle physics.

CPPC-2024-01, MS-TP-24-06

##  1 Introduction

The effective number of neutrinos, $N_{\rm 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\lesssim 1$ MeV), the primary role of $N_{\rm eff}$ is to fix
the universal expansion rate at $T\lesssim 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_{\rm 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].

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.

From the theoretical perspective, the expected value of $N_{\rm 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_{\rm eff}^{\rm SM}=3.0440\pm 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 ${\cal 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 calculation of [28] modelling the same physics—to five
significant digits in the central value and with comparable error estimates.
The precision computation of $N_{\rm eff}^{\rm SM}$ appears therefore to have
reached convergence, at least in a limited sense.

Standard-model corrections to $N_{\rm eff}^{\rm SM}$ | Leading-digit contribution
---|---
$m_{e}/T_{d}$ correction | $+0.04$
$\mathcal{O}(e^{2})$ FTQED correction to the QED EoS | $+0.01$
Non-instantaneous decoupling+spectral distortion | $-0.006$
$\mathcal{O}(e^{3})$ FTQED correction to the QED EoS | $-0.001$
Flavour oscillations | $+0.0005$
Type (a) FTQED corrections to the weak rates | $\lesssim 10^{-4}$
Sources of uncertainty |
Numerical solution by FortEPiaNO | $\pm 0.0001$
Input solar neutrino mixing angle $\theta_{12}$ | $\pm 0.0001$
Table 1: Leading-digit contributions from various SM corrections, in order of
importance, thus far accounted that make up the final $N_{\rm eff}^{\rm
SM}-3=0.0440\pm 0.0002$ [29, 28]. Table adapted from [29].

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_{\rm eff}^{\rm SM}$, 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 [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^{-}\to\nu_{\alpha}\bar{\nu}_{\alpha}$ 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_{\rm eff}^{\rm SM}$ is quite substantial—at the $\sim
0.001$ level. We have reservations about this result: Aside from mapping rate
corrections from an incompatible energy regime (${\cal O}(1)$ keV $\ll m_{e}$
in a stellar plasma versus ${\cal 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 $\nu_{\alpha}e\to\nu_{\alpha}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,111We note that computations of the vacuum corrections to the
neutrino-electron interaction due to diagrams (a)–(d) already appeared
previously in reference [33]. 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 leading-order (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_{\rm eff}^{\rm SM}$, it is clear that
corrections of this magnitude cannot effect a shift in $N_{\rm eff}$ as
sizeable as that claimed in reference [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_{\rm
eff}^{\rm SM}$ at the $\sim 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_{\rm eff}$, $\sigma(N_{\rm eff})\simeq
0.02-0.03$ [35], an accurate theoretical prediction for $N_{\rm eff}^{\rm SM}$
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. 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_{\rm eff}^{\rm SM}$.

  2. 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\sim 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_{\rm eff}^{\rm SM}$.

  3. 3.

Finally, we use our NLO results to estimate the corresponding change in
$N_{\rm eff}^{\rm SM}$. 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_{\rm eff}^{\rm SM}$ 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_{\rm eff}^{\rm SM}$ due to these corrections is presented in
section 4, and section 5 contains our conclusions. Four 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, and parameterisation of the collision integrals.

##  2 Preliminaries

The SM effective number of neutrinos $N_{\rm eff}^{\rm SM}$ is defined via the
ratio of the energy density contained in three generations of neutrinos and
antineutrinos, $\rho_{\nu}$, to the photon energy density $\rho_{\gamma}$ in
the limit $T/m_{e}\to 0$,

$\left.\frac{\rho_{\nu}}{\rho_{\gamma}}\right|_{T/m_{e}\to
0}\equiv\frac{7}{8}\Big{(}\frac{4}{11}\Big{)}^{4/3}N_{\rm eff}^{\rm SM}.$
(2.1)

Here, $T$ is the photon temperature, $m_{e}$ the electron mass, and the limit
$T/m_{e}\to 0$ is understood to apply only to $T\gg m_{i}$, where $m_{i}$ is
the largest neutrino mass.

To compute the precise value of $N_{\rm eff}^{\rm SM}$ requires that we track
the evolution of $\rho_{\nu}$ and $\rho_{\gamma}$ simultaneously across the
the neutrino decoupling epoch, i.e., $T\sim{\cal O}(10)\to{\cal 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

$\frac{{\rm d}\rho_{\rm tot}}{{\rm d}t}+3H(\rho_{\rm tot}+P_{\rm tot})=0,$
(2.2)

where $\rho_{\rm tot}$ and $P_{\rm tot}$ are, respectively, the total energy
density and pressure, $H\equiv(1/a)({\rm d}a/{\rm d}t)$ the Hubble expansion
rate, and $a$ is the scale factor. For the physical system at hand, $\rho_{\rm
tot}\equiv\rho_{\rm QED}+\rho_{\nu}$ and $P_{\rm tot}\equiv P_{\rm
QED}+P_{\nu}$, where $\rho_{\rm QED}\equiv\rho_{\gamma}+\rho_{e}$ subsumes the
photon and the electron/positron energy densities, and similarly for $P_{\rm
QED}$. The assumption of thermodynamic equilibrium in the QED sector in the
time frame of interest means that the standard thermodynamic relation
$\rho_{\rm QED}=-P_{\rm QED}+T\,(\partial/\partial T)P_{\rm QED}$ applies.
Then, the finite-temperature corrections to the QED equation of state
summarised in table 1 simply refer to corrections to the QED partition
function $Z_{\rm QED}$ and hence $P_{\rm QED}=(T/V)\ln Z_{\rm QED}$ that alter
$\rho_{\rm QED}+P_{\rm QED}=T\,(\partial/\partial T)P_{\rm QED}$ from its
ideal-gas expectation. Corrections to $Z_{\rm QED}$ are known to ${\cal
O}(e^{3})$ for arbitrary $m_{e}$ and chemical potential $\mu$ [39] and to
${\cal O}(e^{5})$ for $m_{e}=\mu=0$ [40, 41]. Their effects on $N_{\rm
eff}^{\rm SM}$ have been estimated in references [17, 20, 21] to ${\cal
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, $\varrho=\varrho(p,t)$. Schematically, the QKEs take the
form [42]

$\partial_{t}\varrho-pH\partial_{p}\varrho=-{\rm i}[\mathbb{H},\varrho]+{\cal
I}[\varrho],$ (2.3)

where $\mathbb{H}=\mathbb{H}(p,t)=\mathbb{H}_{\rm vac}+\mathbb{V}$ is the
effective Hamiltonian incorporating vacuum flavour oscillations
$\mathbb{H}_{\rm vac}$ and in-medium corrections from forward scattering
$\mathbb{V}$,222Depending on context, these modifications to the in-medium
quasiparticle dispersion relations are variously known as thermal masses,
matter potentials, or refractive indices. and

${\cal
I}[\varrho]=\frac{1}{2}\Big{(}(\mathbb{1}-\varrho)\mathbb{\Gamma}^{<}-\varrho\mathbb{\Gamma}^{>}\Big{)}+{\rm
h.c.}$ (2.4)

is the collision integral encapsulating the non-unitary gains
($\mathbb{\Gamma}^{<}=\mathbb{\Gamma}^{<}(p,t)$) and losses
($\mathbb{\Gamma}^{>}=\mathbb{\Gamma}^{>}(p,t)$) of $\varrho$. In the context
of precision $N_{\rm eff}^{\rm SM}$, the quantities $\varrho$, $\mathbb{H}$,
and $\mathbb{\Gamma}^{\gtrless}$ are understood to be $3\times 3$ hermitian
matrices in flavour space, with the diagonal entries of $\varrho$
corresponding to the occupation numbers
$f_{\alpha}(p,t)\equiv\\{\varrho(p,t)\\}_{\alpha\alpha}$, for
$\alpha=e,\mu,\tau$.333Strictly speaking, an occupation number refers to a
diagonal element of $\varrho$ in the basis in which $\mathbb{H}$ is diagonal,
i.e., in the mass basis, rather than the interaction basis. This distinction
is however unimportant for the present study, and we follow the common
practice of calling, e.g., $\\{\varrho(p,t)\\}_{ee}$ the electron neutrino
occupation number, etc. We 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 $\mathbb{\Gamma}^{\gtrless}$
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\to 2$ processes
involving (i) two neutrinos and two electrons any way distributed in the
initial and final states (labelled “$\nu e$”), and (ii) neutrino-neutrino
scattering (“$\nu\nu$”). The leading-order $\mathbb{\Gamma}^{\gtrless}$ for
these processes are well known, and take the form of two-dimensional momentum
integrals,

$\displaystyle\mathbb{\Gamma}^{\gtrless}_{\nu e}(p,t)$ $\displaystyle\propto
G_{F}^{2}\int{\rm d}p_{2}\,{\rm d}p_{3}\,\mathbb{\Pi}_{\nu
e}(p,p_{2},p_{3})\,{\cal F}_{\nu
e}^{\gtrless}(\varrho(p_{2},t),f_{e}(p_{3},t),f_{e}(p_{4},t)),$ (2.5)
$\displaystyle\mathbb{\Gamma}^{\gtrless}_{\nu\nu}(p,t)$ $\displaystyle\propto
G_{F}^{2}\int{\rm d}p_{2}\,{\rm
d}p_{3}\,\mathbb{\Pi}_{\nu\nu}(p,p_{2},p_{3})\,{\cal
F}_{\nu\nu}^{\gtrless}(\varrho(p_{2},t),\varrho(p_{3},t),\varrho(p_{4},t)).$

Here, $G_{F}$ is the Fermi constant; $p_{4}$ is fixed by energy-momentum
conservation; $\mathbb{\Pi}_{\nu e}$ and $\mathbb{\Pi}_{\nu\nu}$ are
scattering kernels, which are diagonal in the neutrino interaction basis; and
${\cal F}_{\nu e}^{\gtrless}$ and ${\cal F}_{\nu\nu}^{\gtrless}$ are real-time
matrix products of $\varrho(p_{2},t)$ with $f_{e}(p_{3},t),f_{e}(p_{4},t)$ and
$\varrho(p_{3},t),\varrho(p_{4},t)$, respectively. Because $\varrho$ can
develop non-zero off-diagonal components from neutrino oscillations, ${\cal
F}^{\gtrless}$ and hence $\mathbb{\Gamma}^{\gtrless}$ are generally not
diagonal in the interaction basis.

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_{\rm eff}^{\rm SM}$ computations.

###  2.1 Damping approximation

We would like to compute QED corrections to $\mathbb{\Gamma}^{\gtrless}_{\nu
e}$, and estimate their impact on $N_{\rm eff}^{\rm SM}$. Ideally these
corrections would take the form of corrections to the scattering kernel
$\mathbb{\Pi}_{\nu 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 $\varrho(p)$—are
in thermal equilibrium specified by a common temperature $T$ equal to the
photon temperature. Then, defining $\delta\varrho=\varrho(p)-\varrho_{\rm
eq}(p)$, where $\varrho_{\rm eq}(p)={\rm diag}(f_{F}(p),f_{F}(p),f_{F}(p))$
and $f_{F}(p)=[\exp(p/T)+1]^{-1}$ is the Fermi-Dirac distribution, the
collision integral (2.4) can be expanded to linear order in $\delta\varrho$ to
give

$\\{{\cal
I}[\varrho(p)]\\}_{\alpha\beta}\simeq-\\{D(p,T)\\}_{\alpha\beta}\Big{[}\\{\varrho(p)\\}_{\alpha\beta}-\delta_{\alpha\beta}f_{F}(p)\Big{]},$
(2.6)

where $\delta_{\alpha\beta}$ is a Kronecker-$\delta$, and

$\\{D(p,T)\\}_{\alpha\beta}\equiv\frac{1}{2}\Big{[}\\{\mathbb{\Gamma}(p,T)\\}_{\alpha\alpha}+\\{\mathbb{\Gamma}(p,T)\\}_{\beta\beta}\Big{]}$
(2.7)

are the damping coefficients composed of components of the damping rate matrix

$\displaystyle\mathbb{\Gamma}(p,T)\equiv\mathbb{\Gamma}^{>}(p,T)+\mathbb{\Gamma}^{<}(p,T).$
(2.8)

Linearisation in $\delta\varrho$ ensures that $\mathbb{\Gamma}(p,T)$ is
diagonal in the neutrino interaction basis, and that
$\mathbb{\Gamma}^{>}(p,T)=e^{E_{p}/T}\mathbb{\Gamma}^{<}(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
$\mathbb{\Gamma}$ 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 $\\{\mathbb{\Gamma}_{\nu e}(p,T)\\}_{\alpha\alpha}$ due to
the closed fermion loop (diagram (d) in figure 1). Our estimates of its impact
on $N_{\rm eff}^{\rm SM}$ will be presented in section 4.

##  3 Calculation and NLO results

We compute QED corrections to the interaction rates of the weak processes
$\nu_{\alpha}e\to\nu_{\alpha}e$ and
$\nu_{\alpha}\bar{\nu}_{\alpha}\leftrightarrow e^{+}e^{-}$ within the
framework of Fermi theory, whose effective Lagrangian can be expressed as

$\displaystyle\mathcal{L}_{\mathrm{4F}}=\frac{4G_{F}}{\sqrt{2}}\left[\bar{\psi}_{\alpha}\gamma_{\mu}P_{L}\psi_{\alpha}\right]\left[g_{L}^{\alpha}J_{L}^{\mu}+g_{R}J_{R}^{\mu}\right].$
(3.1)

Here, the spinors $\psi_{\alpha}$ represent neutrino fields of flavour
$\alpha$; $J^{\mu}_{L/R}=\bar{\psi}_{e}\gamma^{\mu}P_{L/R}\psi_{e}$ are the
left- and right-handed electron current operators, with the chiral projectors
$P_{L/R}=\frac{1}{2}(1\mp\gamma^{5})$; the right-handed coupling reads
$g_{R}=\sin^{2}\theta_{W}$, while the left-handed coupling
$g^{\alpha}_{L}=-\frac{1}{2}+\sin^{2}\theta_{W}+\delta^{\alpha e}$ depends on
the neutrino flavor $\alpha=e,\mu,\tau$. 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^{\alpha}_{L}$ and $g_{R}$ updated for the quarks of
interest.444The reason is that a Fierz identity has been used to bring the
charged and neutral current contributions into the common form (3.1). We shall
however omit the contributions from quarks in this work: at finite
temperatures these contributions are Boltzmann-suppressed at the ${\cal 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^{\alpha}_{V,A}=\frac{1}{2}(g^{\alpha}_{L}\pm g_{R})$ as an
alternative notation.

###  3.1 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 $\varrho$, i.e.,
$\Gamma_{\alpha}(p,T)\equiv\\{\mathbb{\Gamma}_{\nu e}(p,T)\\}_{\alpha\alpha}$,
which is given by the sum

$\Gamma_{\alpha}=\Gamma_{\alpha}^{>}(p)+\Gamma_{\alpha}^{<}(p),$ (3.2)

where $\Gamma_{\alpha}^{\gtrless}$ are the $\alpha$-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 $\Sigma^{\gtrless}$,555There
is no term involving the plasma four-velocity $U$ because we consider the sum
over helicities. See, e.g., appendix C of [44] in the context of heavy
neutrinos.

$\Gamma_{\alpha}^{\gtrless}(p)=\mp\frac{1}{2p^{0}}\left.\mathrm{Tr}\left[\not{p}\Sigma_{\alpha}^{\gtrless}\right]\right|_{p^{0}=E_{p}},$
(3.3)

where $\Sigma^{\gtrless}$ are identified with the self-energies with opposite
(Keldysh) contour indices:666The idea to evaluate fields on a closed time path
was already implicit in Schwinger’s seminal work [45]. Here, we follow common
convention and refer to the indices marking the forward and backward contours
in the real-time formalism as Keldysh indices [46].
$\Sigma^{<}\equiv\Sigma^{12}$ and $\Sigma^{>}\equiv\Sigma^{21}$. In thermal
equilibrium, detailed balance is established by the fermionic Kubo-Martin-
Schwinger (KMS) relation $\Sigma^{>}=-e^{p^{0}/T}\Sigma^{<}$. Then, the mode-
dependent damping rate can also be written as

$\Gamma_{\alpha}(p)=\frac{1}{2p^{0}}\left.\mathrm{Tr}\left[\not{p}(\Sigma_{\alpha}^{<}-\Sigma_{\alpha}^{>})\right]\right|_{p^{0}=E_{p}}=\frac{1}{2p^{0}f_{F}(p^{0})}\left.\mathrm{Tr}\left[\not{p}\Sigma_{\alpha}^{<}\right]\right|_{p^{0}=E_{p}},$
(3.4)

where the first equality follows from the fact that $\Gamma_{\alpha}$ can be
related to the discontinuity of the retarded neutrino self-energy evaluated at
the quasiparticle pole, $-i{\rm disc}\Sigma^{R}=2{\rm
Im}\Sigma^{R}=i(\Sigma^{>}-\Sigma^{<}$), in accordance with the optical
theorem at finite temperature. While the KMS relation makes explicit use of
the fact that $\Sigma^{\gtrless}$ are computed in thermal equilibrium, the
real-time formalism used here can be generalised in a straightforward manner
to non-equilibrium situations [37].777Specifically, this entails replacing the
equilibrium distribution in the relation (A.4) between the statistical and the
spectral propagators with a dynamical function, and modifying all other
propagators accordingly. See, e.g. [47, 48] for a discussion. Where there is
no confusion, we shall drop the flavour index $\alpha$ in the following.

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 $\nu_{\alpha}$).

On general grounds the dominant QED correction to $\Gamma_{\alpha}$ 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.888 Resumming the photon
propagator (see discussion above equation (3.13) and appendix C) opens up an
additional cut through the photon and the internal neutrino lines,
corresponding to the so-called plasmon process $\gamma\to\nu\bar{\nu}$, even
if the imaginary part of the photon self-energy is neglected. This can be seen
in equation (C.32), which exhibits an on-shell $\delta$-function in the limit
of vanishing width (C.34). However, the plasmon process contributes at ${\cal
O}(\alpha_{\rm em}^{2})$ (see appendix C in [49]), which can be understood as
a phase space suppression from $m_{\gamma}\sim\alpha_{\rm em}T$. We therefore
do not consider it in this work. We shall compute only this contribution and
denote the corresponding Wightman self-energy with $\Sigma^{ba}_{\rm NLO}$.
For notational convenience, we further split $\Sigma^{ba}_{\rm NLO}$ into a
sum $\Sigma^{ba}_{\rm NLO}=\sum_{c,d=1,2}\Sigma^{ba,cd}$ over the internal
real-time indices $c$ and $d$, with the partial self-energies $\Sigma^{ba,cd}$
given by

$\mathrm{Tr}\left[\not{p}\Sigma^{ba,cd}\right]=(-1)^{a+b+c+d+1}(ie)^{2}\left(\frac{i4G_{F}}{\sqrt{2}}\right)^{2}\\\
\times\int_{l,q,k}\mathrm{Tr}\left[iS_{e}^{ad}(q)\gamma_{\rho}(P_{L}g_{L}^{\alpha}+P_{R}g_{R})iS_{e}^{da}(q+p-l)\gamma^{\mu}\right]iD_{\mu\nu}^{cd}(p-l)\\\
\times\mathrm{Tr}\left[\not{p}\gamma^{\rho}P_{L}iS_{\nu_{\alpha}}^{ba}(l)\gamma^{\sigma}P_{L}\right]\mathrm{Tr}\left[iS_{e}^{cb}(k)\gamma^{\nu}iS_{e}^{bc}(k+p-l)\gamma_{\sigma}(P_{L}g_{L}^{\alpha}+P_{R}g_{R})\right].$
(3.5)

Here, we have introduced the abbreviation
$\int_{p}\equiv\int\differential[4]{p}/(2\pi)^{4}$ for loop integrals; the
traces are over the Clifford algebra; the subscripts “$e$” and
“$\nu_{\alpha}$” 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 $\Sigma^{<}=\Sigma^{12}_{\rm NLO}$, we set the
external contour indices to $a=2$ and $b=1$. Then, the contributions
$\mathrm{Tr}\left[\not{p}\Sigma^{12,21}\right]$ and
$\mathrm{Tr}\left[\not{p}\Sigma^{12,12}\right]$, 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
$\mathrm{Tr}\left[\not{p}\Sigma^{12,11}\right]=\mathrm{Tr}\left[\not{p}\Sigma^{12,22}\right]^{\ast}$.
It then follows that

$\mathrm{Tr}\left[\not{p}\Sigma^{12}_{\rm
NLO}\right]=2\mathrm{Re}\,\mathrm{Tr}\left[\not{p}\Sigma^{12,11}\right],$
(3.6)

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\to 2$ process represented by diagram (d) of
figure 1 via

$\displaystyle p_{1}=p,\qquad p_{2}=-l,\qquad p_{3}=q+p-l,\qquad p_{4}=-q.$
(3.7)

Then, writing out the propagators $iS_{e}^{12/21}$ and $iS_{\nu}^{12/21}$
explicitly, the self-energy (3.5)–(3.6) can be brought into the form

$\mathrm{Tr}\left[\not{p}\Sigma^{12}_{\mathrm{NLO}}\right]=-(2\pi)^{3}\int_{q,l}\delta(p_{2}^{2})\delta(p_{3}^{2}-m_{e}^{2})\delta(p_{4}^{2}-m_{e}^{2})\mathcal{F}(p_{2}^{0},p_{3}^{0},p_{4}^{0})\mathcal{T}_{\rm
NLO},$ (3.8)

where the population factor

$\mathcal{F}(p_{2}^{0},p_{3}^{0},p_{4}^{0})=\left[f_{F}(|p_{2}^{0}|)-\Theta(p_{2}^{0})\right]\left[f_{F}(|p_{3}^{0}|)-\Theta(-p_{3}^{0})\right]\left[f_{F}(|p_{4}^{0}|)-\Theta(-p_{4}^{0})\right]$
(3.9)

contains the equilibrium phase space distributions of the three integrated
external particles and is analogous to $\cal{F}^{\gtrless}$ in equation (2.5).
The quantity ${\cal T}_{\rm NLO}$ plays the role of QED corrections to the LO
squared matrix element,999For completeness,
$\mathcal{T}_{\mathrm{LO}}=2^{7}G_{F}^{2}\left[(g_{L}^{\alpha})^{2}\left(p_{1}\cdot
p_{4}\right)^{2}+g_{R}^{2}\left(p_{2}\cdot
p_{4}\right)^{2}+g_{L}^{\alpha}g_{R}m_{e}^{2}\left(p_{1}\cdot
p_{2}\right)\right]$. If we were to replace ${\cal T}_{\rm NLO}$ with ${\cal
T}_{\rm LO}$ in equation (3.8), we would recover the leading-order equilibrium
neutrino production rates. and can be written in terms of the one-loop photon
self-energy $\Pi_{ab}^{\mu\nu}$ as

$\mathcal{T}_{\rm
NLO}=-16G_{F}^{2}g_{V}^{\alpha}\mathrm{Re}\left\\{\mathrm{Tr}\left[(\not{p}_{3}+m_{e})\gamma^{\mu}\
(m_{e}-\not{p}_{4})\gamma^{\rho}\left(g_{L}^{\alpha}P_{L}+g_{R}P_{R}\right)\right]\right.\\\
\left.\times\mathrm{Tr}\left[\not{p}_{1}\gamma_{\rho}P_{L}\not{p}_{2}\gamma_{\sigma}P_{L}\right]D^{11}_{\mu\nu}(P)\Pi_{11}^{\nu\sigma}(P)\right\\},$
(3.10)

where we have introduced the photon momentum $P=p_{1}+p_{2}$.

Since ${\cal T_{\rm NLO}}$ has the interpretation of a squared matrix element,
we can split it into a vacuum and a thermal part according to temperature
dependence,

${\cal T}_{\rm NLO}={\cal T}_{\rm vac}+{\cal T}_{\rm th}.$ (3.11)

The vacuum correction ${\cal T}_{\rm 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 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

$\displaystyle\mathcal{T}_{\rm vac}=-2^{10}G_{F}^{2}\,\alpha_{\rm
em}\,[(g^{\alpha}_{V})^{2}/(4\pi)]\left[m^{2}_{e}(p_{1}\cdot
p_{2})+2(p_{1}\cdot p_{4})^{2}\right]\mathrm{Re}\,\Pi_{2}(P^{2})$ (3.12)

where $\alpha_{\rm em}$ is the electromagnetic fine-structure constant, and
the form factor $\Pi_{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}\leftrightarrow p_{4}$. This
symmetry, along with momentum conservation, also ensures the absence of all
antisymmetric terms containing Levi-Civita symbols arising from traces of
$\gamma^{5}$ with four or more $\gamma$-matrices. Vacuum QED corrections to
the neutrino-electron interaction matrix element were previously computed in
reference [33].

The thermal correction ${\cal T}_{\rm 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 $\delta(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 ${\cal T}_{\rm NLO}$: Terms proportional to $f_{B}$ correspond to putting
the photon line on-shell, which are forbidden for diagram (d) for kinematic
reasons [51, 49]. 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 ${\cal T}_{\rm NLO}$ contains a $t$-channel contribution from
elastic $\nu_{\alpha}e$ scattering (i.e., $p_{2}^{0}<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}_{\mu\nu}$ in equation (3.10) with the
fully-resummed photon propagator $\bar{D}^{ab}_{\mu\nu}$.101010In principle,
it is also necessary to replace the tree-level photon propagator in equation
(3.12) with the fully resummed one. However, since the vacuum correction is
not dominated by small photon momenta, the choice of photon propagator has in
this case a negligible impact on the numerical outcome. Furthermore, because
both $\bar{D}^{ab}_{\mu\nu}$ and $\Pi_{11}^{\nu\sigma}$ split into a
longitudinal (“$L$”) and a transverse (“$T$”) part, the same decomposition
applies also to ${\cal T}_{\rm th}$, i.e., $\mathcal{T}_{\rm
th}=\mathcal{T}^{L}_{\rm th}+\mathcal{T}^{T}_{\rm th}$, where
$\mathcal{T}^{L,T}_{\rm th}$ can be brought into the form

$\mathcal{T}^{L/T}_{\rm
th}=-2^{8}G_{F}^{2}(g_{V}^{\alpha})^{2}\mathrm{Re}\,\bar{D}^{L/T}_{R}(P)\mathrm{Re}\,\Pi^{L/T}_{R,T\neq
0}(P)\left[2(p_{1}\cdot p_{4})P^{1,4}_{L/T}\right.\\\ \left.+(p_{1}\cdot
p_{2})(a_{L/T}(p_{1}\cdot p_{2})+P^{2,2}_{L,T}+P^{4,4}_{L/T})\right].$ (3.13)

Here, $\bar{D}_{R}$ denotes the retarded resummed photon propagator;
$\Pi_{R,T\neq 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^{\mu\nu}_{L/T}p_{i,\mu}p_{j,\nu}$; and $a_{L,T}=(3\mp 1)/4$
further differentiates between the longitudinal and the transverse
contribution.

Note that the imaginary part of the photon propagator,
$\mathrm{Im}\,\bar{D}^{L,T}_{11}$, does not appear in equation (3.13) because
it is formally of higher-order in $\alpha_{\rm em}$ and we are only interested
in the $\order{\alpha_{\rm em}}$ corrections. We also only use the resummed
photon propagator in the IR divergent $t$-channel contribution; where the
divergence is absent, i.e., in the $s$-channel, we set $\bar{D}_{R}\to D_{R}$,
where $D_{R}$ is the un-resummed counterpart of $\bar{D}_{R}$. The final
expressions for $\mathrm{Re}\Pi_{11,T\neq 0}^{L/T}$ and
$\mathrm{Re}\,\bar{D}^{L/T}_{11}$ are given in equations (C.10), (C.12) and
(C.32), which can be easily mapped to $\mathrm{Re}\,\Pi_{R,T\neq 0}^{L/T}$ and
$\mathrm{Re}\,\bar{D}_{R}$ via
$\mathrm{Re}\,\bar{D}^{L/T}_{R}=\mathrm{Re}\,\bar{D}^{L/T}_{11}$ and
$\mathrm{Re}\,\Pi^{L/T}_{R}=\mathrm{Re}\,\Pi^{L/T}_{11}$. Details of their
derivation, along with the relevant thermal integrals, can be found in
appendices B and C.

###  3.2 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:

  * •

Fermi’s constant:
$G_{F}=$1.166\,378\,8(6)\text{\,}\mathrm{M}\mathrm{eV}^{-2}$$,

  * •

Electron mass: $m_{e}=$0.510\,998\,950\,00(15)\text{\,}\mathrm{MeV}$$,

  * •

Electromagnetic fine-structure constant: $\alpha_{\rm
em}^{-1}(0)=$137.035\,999\,180(10)\text{\,}$$, and

  * •

Weinberg angle:
$\sin^{2}\theta_{W}(0)_{\overline{\mathrm{MS}}}=$0.238\,63(5)\text{\,}$$.

The renormalisation scale $\mu_{R}$ appearing in the photon self-energy of the
vacuum contribution is identified with the electron mass, $\mu_{R}=m_{e}$.

Figure 3: Top: NLO contributions to the $\nu_{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}\neq 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\gg 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 $\mu_{R}$. The lower panel shows the LO neutrino
interaction rate compared to the Hubble rate, and $T_{d}$ indicates the
decoupling temperature, defined via $\Gamma_{\alpha}(T_{d(e)})=H(T_{d(e)})$.
Bottom: Same as top, but for $\nu_{\mu,\tau}$.

Figure 3 shows the closed fermion loop corrections to the damping rates
$\Gamma_{e}(p)$ and $\Gamma_{\mu,\tau}(p)$ at the mean momentum $p=3.15T$.
Relative to their respective LO rates, we find the corrections at temperatures
$T\sim 1\to 3$ MeV to fall in the range $-0.2\to+0.1\%$ and
$-0.0005\to+0.0002\%$, respectively, for $\nu_{e}$ and $\nu_{\mu,\tau}$. We
further note that:

  1. 1.

At $T\sim$2\text{\,}\mathrm{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. 2.

Reference [30] also found no significant flavour dependence in the rate
corrections: their $\Gamma_{e}$ and $\Gamma_{\mu,\tau}$ 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^{\pm}$ thermal
bath. This difference renders the corresponding vector couplings,
$g_{V}^{e}\sim 0.49$ and $g_{V}^{\mu,\tau}\sim-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. 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}\rightarrow 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\gtrsim 3m_{e}$, but becomes sizeable at low temperatures. In particular, in
the limit $T\rightarrow 0$ the ratio $\Gamma_{\alpha}^{\rm
NLO}/\Gamma_{\alpha}^{\rm LO}$ vanishes for $m_{e}>0$, but diverges for
$m_{e}=0$ because of Boltzmann suppression of the LO rates. Precision
computations of $N_{\rm eff}^{\rm SM}$ 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\sim(2\to 3)\times m_{e}$, what we assume for $m_{e}$ in the
NLO rate computations may yet have some impact on $N_{\rm eff}$ (see section
4.2).

Figure 4: _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 1PI-resummed (green) and HTL (magenta) photon
propagators, in each case with $m_{e}\neq 0$ (solid lines) or $m_{e}=0$
(dashed lines) in both electron loops of the self-energy diagram of figure 2.
As in figure 3, we normalise all curves to the LO rate, which we always
evaluate with $m_{e}\neq 0$. The lower panel shows the ratio of the 1PI result
for $m_{e}\neq 0$ to the HTL result for $m_{e}=0$ as a function of the
temperature. _Bottom_ : Same as top, but for $\alpha\neq e$.

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}\rightarrow 0$, (iii) using the
HTL photon propagator (which does not depend on the electron mass), but with
$m_{e}$ everywhere else,111111Even though it is physically inconsistent to
keep $m_{e}$ in the electron propagators and not in the photon one, (iii) is
nonetheless helpful as it isolates the dependence of the rate corrections to
the 1PI-resummed vs the HTL propagator. 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\lesssim 1$ MeV, which
can be explained by the fact that the HTL approximation only holds for $T\gg
m_{e}$. In the lower panel of figure 4, we highlight the impact of $m_{e}$ by
displaying the ratio of (i) to (iv).

###  3.3 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 $H^{2}(T)=(\rho_{\rm
QED}+\rho_{\nu})/(3M^{2}_{\rm Pl})$, where $M_{\rm
Pl}=$2.435\,36\text{\times}{10}^{21}\text{\,}\mathrm{MeV}$$ is the reduced the
Planck mass, $\rho_{\nu}=3g_{\nu}(7/8)(\pi^{2}/30)T^{4}$ is the energy density
in three families of neutrinos with $g_{\nu}=2$, and

$\displaystyle\rho_{\rm
QED}=g_{\gamma}\frac{\pi^{2}}{30}T^{4}+\frac{g_{e}}{2\pi^{2}}\int_{0}^{\infty}\differential{p}p^{2}E_{e}f_{F}(E_{e})$
(3.14)

is the energy density in the QED plasma assuming an ideal gas configuration,
with $g_{\gamma}=2$, $g_{e}=4$, and $E_{e}=\sqrt{p^{2}+m_{e}^{2}}$. By solving
the equation $\Gamma_{\alpha}^{\rm LO}(T_{d(\alpha)})=H(T_{d(\alpha)})$ for
the flavour-dependent decoupling temperatures $T_{d(\alpha)}$, we determine
the LO decoupling temperatures to be

$\displaystyle T^{\mathrm{LO}}_{d(e)}\simeq$1.320\,58\text{\,}\mathrm{MeV}$,$
(3.15) $\displaystyle
T^{\mathrm{LO}}_{d(\mu,\tau)}\simeq$2.221\,722\,09\text{\,}\mathrm{MeV}$.$

The first number differs from the earlier estimate of
$T^{\mathrm{LO}}_{d(e)}\simeq 1.3453$ MeV from reference [21], which can be
traced to a different choice of the Weinberg angle. Incorporating QED
corrections to the damping rates, the decoupling temperatures shift to

$\displaystyle T^{\mathrm{NLO}}_{d(e)}\simeq$1.321\,10\text{\,}\mathrm{MeV}$,$
(3.16) $\displaystyle
T^{\mathrm{NLO}}_{d(\mu,\tau)}\simeq$2.221\,722\,11\text{\,}\mathrm{MeV}$$

corresponding to an increase of $\sim$0.04\text{\,}\mathrm{\char 37\relax}$$
for $\nu_{e}$ and of $\sim$8\text{\times}{10}^{-07}\text{\,}\mathrm{\char
37\relax}$$ for $\nu_{\mu,\tau}$. 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^{\mathrm{NLO}}_{d(\mu,\tau)}$.

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 computations of the $m_{e}=0$ rate corrections (but with
$m_{e}\neq 0$ LO rates), we find the corresponding NLO decoupling temperatures
to be

$\displaystyle
T^{\mathrm{NLO},m_{e}=0}_{d(e)}\simeq$1.321\,18\text{\,}\mathrm{MeV}$,$ (3.17)
$\displaystyle
T^{\mathrm{NLO},m_{e}=0}_{d(\mu,\tau)}\simeq$2.221\,722\,23\text{\,}\mathrm{MeV}$,$

i.e., a $0.05\text{\,}\mathrm{\char 37\relax}$ and
$6\text{\times}{10}^{-06}\text{\,}\mathrm{\char 37\relax}$ shift for $\nu_{e}$
and $\nu_{\mu,\tau}$, respectively, relative to their corresponding LO
decoupling temperature (3.15).

##  4 NLO effects on $N_{\rm eff}^{\rm SM}$

Having computed the closed fermion loop correction to the damping rate
$\Gamma_{\alpha}$, we are now in a position to estimate its effect on $N_{\rm
eff}^{\rm SM}$. Within the damping approximation, two avenues are available to
us:

  * •

Given $\Gamma_{\alpha}$ at a representative momentum, we have already
estimated in equation (3.16) the corresponding correction to the neutrino
decoupling temperature $T_{d}$, defined via $\Gamma_{\alpha}(T_{d})=H(T_{d})$.
Then, under the assumption of instantaneous neutrino decoupling, an estimate
of the change to $N_{\rm eff}^{\rm SM}$, $\delta N_{\rm eff}\equiv N_{\rm
eff}^{\rm NLO}-N_{\rm eff}^{\rm LO}$, can be obtained via entropy conservation
arguments.

  * •

We may also compute $\delta N_{\rm 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
$\delta N_{\rm eff}$ are summarised in table 2.

Method | $\delta N_{\rm eff}/N^{\rm LO}_{\rm eff}$ | $\delta N_{\rm eff}^{m_{e}=0}/N^{\rm LO}_{\rm eff}$
---|---|---
Entropy conservation (common decoupling) | $-1.1\text{\times}{10}^{-05}\text{\,}$ | $-1.2\text{\times}{10}^{-05}\text{\,}$
Entropy conservation (flavour-dependent decoupling) | $-5.4\text{\times}{10}^{-06}\text{\,}$ | $-6.1\text{\times}{10}^{-06}\text{\,}$
Boltzmann damping approximation (mean momentum) | $-7.8\text{\times}{10}^{-06}\text{\,}$ | $-1.6\text{\times}{10}^{-05}\text{\,}$
Boltzmann damping approximation (full momentum) | $-7.9\text{\times}{10}^{-06}\text{\,}$ | $-2.6\text{\times}{10}^{-05}\text{\,}$
Table 2: Estimates of the relative correction to $N^{\rm SM}_{\rm eff}$ due to
NLO weak rate corrections, with and without the electron mass, using different
methods.

###  4.1 Entropy conservation

The entropy conservation argument posits that entropies in a comoving volume
in two decoupled sectors are separately conserved, i.e.,

$\displaystyle s(a_{1})a_{1}^{3}$ $\displaystyle=s(a_{2})a_{2}^{3},$ (4.1)
$\displaystyle s_{\nu_{\alpha}}(a_{1})a_{1}^{3}$
$\displaystyle=s_{\nu_{\alpha}}(a_{2})a_{2}^{3},$

where $s_{\nu_{\alpha}}$ and $s\equiv
s_{\gamma}+s_{e}+\sum_{\beta\neq\alpha}s_{\nu_{\beta}}+\cdots$ denote,
respectively the entropy density of a decoupled neutrino species
$\nu_{\alpha}$ and of the QED plasma plus any neutrino species $\nu_{\beta}$
that may still be 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 $\nu_{\alpha}$ decouples instantaneously from
the QED plasma, i.e., $a_{1}=a_{d}^{+}$, where the sector temperatures satisfy
$T(a_{1})=T_{\nu_{\alpha}}(a_{1})\equiv 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_{\nu_{\alpha}}(a)\propto T_{\nu_{\alpha}}^{3}(a)$. It
then follows simply from (4.1) that
$T_{\nu_{\alpha}}(a_{2})=(a_{1}/a_{2})T_{\nu_{\alpha}}(a_{1})=(a_{1}/a_{2})T_{d}$.
In the case of the QED+$\nu_{\beta}$ plasma, we parameterise its entropy
density as

$s(a)=\frac{2\pi^{2}}{45}\tilde{g}_{s}(a)T^{3}(a),$ (4.2)

where the entropy degree of freedom parameter $\tilde{g}_{s}$ is given in the
ideal gas limit by

$\tilde{g}_{s}(a)\equiv
g_{\gamma}+\frac{45}{4\pi^{4}}\frac{g_{e}}{T^{4}(a)}\int_{0}^{\infty}{\rm
d}p\,p^{2}\left(E_{e}+\frac{p^{2}}{3E_{e}}\right)f_{F}(E_{e},T(a))+\frac{7}{8}\sum_{\beta\neq\alpha}g_{\nu_{\beta}},$
(4.3)

with $g_{\gamma}=2$, $g_{e}=4$, and $g_{\nu_{\beta}}=2$. For our current
purpose of estimating $\delta N_{\rm eff}$ due to NLO contributions to the
weak rates, it suffices to use the ideal-gas $\tilde{g}_{s}$. We do note
however that QED entropy density at $T_{d}\sim{\cal O}(1)$ MeV is subject in
principle to a sizeable finite-temperature correction to the QED equation of
state, which needs to be included in precision $N_{\rm eff}^{\rm SM}$
calculations. See, e.g., reference [21] for details.

Then, combining the above, we find an estimate of the neutrino-to-photon
temperature ratio of

$\frac{T_{\nu_{\alpha}}(a_{2})}{T(a_{2})}=\left[\frac{\tilde{g}_{s}(a_{2})}{\tilde{g}_{s}(a_{d}^{+})}\right]^{1/3}$
(4.4)

at the later time $a_{2}$. We use this temperature ratio in the following to
provide two estimates of $\delta N_{\rm 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 $\nu_{e}$ decoupling, and set $a_{2}$
at a time significantly after $e^{\pm}$ annihilation, where $T(a_{2})\ll
m_{e}$. The latter leads immediately to an entropy degree of freedom of
$\tilde{g}_{s}(a_{2})=g_{\gamma}=2$, while for the former we have

$\tilde{g}_{s}(a_{d}^{+})=2+\frac{45}{\pi^{4}T_{d(e)}^{4}}\int_{0}^{\infty}{\rm
d}p\,p^{2}\left(E_{e}+\frac{p^{2}}{3E_{e}}\right)f_{F}(E_{e},T_{d(e)}).$ (4.5)

Then, using the temperature ratio (4.4) and the ideal-gas relations
$\rho_{\gamma}\propto g_{\gamma}T^{4}$ and
$\rho_{\nu_{\alpha}}\propto(7/8)g_{\nu_{\alpha}}T_{\nu_{\alpha}}^{4}$, we find

$\left.\frac{\rho_{\nu}}{\rho_{\gamma}}\right|_{T/m_{e}\to
0}=\sum_{\alpha}\frac{\rho_{\nu_{\alpha}}(a_{2})}{\rho_{\gamma}(a_{2})}=3\times\frac{7}{8}\left[\frac{2}{\tilde{g}_{s}(a_{d}^{+})}\right]^{4/3},$
(4.6)

or, equivalently in terms of $N_{\rm eff}$,

$N_{\rm
eff}=3\times\left[\frac{11}{4}\frac{2}{\tilde{g}_{s}(a_{d}^{+})}\right]^{4/3},$
(4.7)

where the factor $4/11$ corresponds to $2/\tilde{g}_{s}(T_{d})$ evaluated in
the limit $T_{d}\gg m_{e}$.

Using the LO and NLO electron neutrino decoupling temperatures in (3.15) and
(3.16) respectively, we find a fractional shift of $\delta N_{\rm eff}/N_{\rm
eff}^{\rm LO}\simeq$-1.1\text{\times}{10}^{-05}\text{\,}$$ due to the rate
corrections. Had we instead used the $m_{e}=0$ NLO decoupling temperatures
(3.17), the correspond shift in $N_{\rm eff}$ would have been $\delta
N^{m_{e}=0}_{\rm eff}\simeq$-3.8\text{\times}{10}^{-05}\text{\,}$$, or $\delta
N_{\rm eff}^{m_{e}=0}/N_{\rm eff}^{\rm
LO}\simeq$-1.2\text{\times}{10}^{-05}\text{\,}$$. Thus, while setting
$m_{e}=0$ leads to a $\sim 10\%$ change in the estimate of $\delta N_{\rm
eff}$, its ultimate impact on $N_{\rm eff}^{\rm SM}$ 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, $\nu_{e}$ and $\nu_{\mu,\tau}$
decouple at different temperatures, $T_{d(e)}$ and $T_{d(\mu,\tau)}>T_{d(e)}$.
Then, to estimate $\delta N_{\rm eff}$ requires that we consider entropy
conservation across four epochs: the time immediately after $\nu_{\mu,\tau}$
decoupling $a_{1}=a_{d(\mu,\tau)}^{+}$; immediately before $\nu_{e}$
decoupling $a_{2}=a_{d(e)}^{-}$; immediately after $\nu_{e}$ decoupling
$a_{3}=a_{d(e)}^{+}$; and $a_{4}$ is a time significantly after $e^{\pm}$
annihilation. The corresponding entropy degrees of freedom are

$\displaystyle\tilde{g}_{s}(a_{d(\mu,\tau)}^{+})$
$\displaystyle=2+\frac{45}{\pi^{4}T_{d(\mu,\tau)}^{4}}\int_{0}^{\infty}{\rm
d}p\,p^{2}\left(E_{e}+\frac{p^{2}}{3E_{e}}\right)f_{F}(E_{e},T_{d(\mu,\tau)})+\frac{7}{8}g_{\nu_{e}},$
(4.8) $\displaystyle\tilde{g}_{s}(a_{d(e)}^{-})$
$\displaystyle=2+\frac{45}{\pi^{4}T_{d(e)}^{4}}\int_{0}^{\infty}{\rm
d}p\,p^{2}\left(E_{e}+\frac{p^{2}}{3E_{e}}\right)f_{F}(E_{e},T_{d(e)})+\frac{7}{8}g_{\nu_{e}},$
$\displaystyle\tilde{g}_{s}(a_{d(e)}^{+})$
$\displaystyle=2+\frac{45}{\pi^{4}T_{d(e)}^{4}}\int_{0}^{\infty}{\rm
d}p\,p^{2}\left(E_{e}+\frac{p^{2}}{3E_{e}}\right)f_{F}(E_{e},T_{d(e)}),$

and $\tilde{g}_{s}(a_{4})=2$. An estimate of the $\nu_{e}$-to-photon energy
density ratio at $a_{4}$ follows straightforwardly from the temperature ratio
(4.4) and ideal-gas temperature-energy relations:

$\frac{\rho_{\nu_{e}}(a_{4})}{\rho_{\gamma}(a_{4})}=\frac{7}{8}\left[\frac{\tilde{g}_{s}(a_{4})}{\tilde{g}_{s}(a_{d(e)}^{+})}\right]^{4/3}=\frac{7}{8}\left[\frac{2}{\tilde{g}_{s}(a_{d(e)}^{+})}\right]^{4/3}.$
(4.9)

For $\rho_{\nu_{\mu,\tau}}(a_{4})/\rho_{\gamma}(a_{4})$, we note that
$\nu_{e}$ decoupling at $a_{d(e)}$ introduces a discontinuity in
$\tilde{g}_{s}$, thereby leading to a more complicated energy density ratio at
$a_{4}$,

$\frac{\rho_{\nu_{\mu,\tau}}(a_{4})}{\rho_{\gamma}(a_{4})}=\frac{7}{8}\left[\frac{2}{\tilde{g}_{s}(a_{d(e)}^{+})}\frac{\tilde{g}_{s}(a_{d(e)}^{-})}{\tilde{g}_{s}(a_{d(\mu,\tau)}^{+})}\right]^{4/3}.$
(4.10)

Then, combining equations (4.9) and (4.10) to form
$\rho_{\nu}=\sum_{\alpha}\rho_{\nu_{\alpha}}$, we find

$N_{\rm
eff}=\left(\frac{11}{4}\frac{2}{\tilde{g}_{s}(a_{d(e)}^{+})}\right)^{4/3}\left[1+2\times\left(\frac{\tilde{g}_{s}(a_{d(e)}^{-})}{\tilde{g}_{s}(a_{d(\mu,\tau)}^{+})}\right)^{4/3}\right]\,$
(4.11)

for our estimate.

Evaluating (4.11) for the LO decoupling temperatures (3.15), we find $N^{\rm
LO}_{\rm eff}\simeq 3.027867$. Including QED corrections to the interaction
rates leads to a shift of $\delta N_{\rm
eff}\simeq$-1.6\text{\times}{10}^{-05}\text{\,}$$ (or $\delta N_{\rm
eff}/N_{\rm eff}^{\rm LO}\simeq$-5.4\text{\times}{10}^{-06}\text{\,}$$) if we
use the $m_{e}\neq 0$ corrections, or $\delta N_{\rm
eff}\simeq$-1.9\text{\times}{10}^{-05}\text{\,}$$ ($\delta N^{m_{e}=0}_{\rm
eff}/N_{\rm eff}^{\rm LO}\simeq$-6.1\text{\times}{10}^{-06}\text{\,}$$ using
the $m_{e}=0$ corrections. Compared with the common-decoupling estimates
displayed in table 2, we note that the flavour-dependent decoupling estimates
of $\delta N_{\rm eff}/N_{\rm eff}^{\rm LO}$ are generally about a factor of
two smaller, in both the $m_{e}=0$ and $m_{e}\neq 0$ cases. This difference is
to be expected, as the QED corrections to the $\nu_{\mu,\tau}$ interaction
rates are negligible compared with the corrections to the $\nu_{e}$ rates. We
emphasise, however, both estimates are very crude approximations: the true
shift in $\delta N_{\rm eff}$ will probably fall somewhere in-between.

###  4.2 Solving the neutrino Boltzmann equations and the continuity equation

Following reference [20], we introduce the comoving quantities $x\equiv
m_{e}\,a$, $y\equiv a\,p$, and $z\equiv T\,a$, and rewrite the continuity
equation (2.2) as a differential equation for the quantity $z$ (corresponding
to the photon temperature):

$\displaystyle\derivative{z}{x}$
$\displaystyle=\frac{\frac{x}{z}J(x/z)-\frac{1}{2z^{3}}\derivative{\bar{\rho}_{\nu}}{x}+G_{1}(x/z)}{\frac{x^{2}}{z^{2}}J(x/z)+Y(x/z)+\frac{2\pi^{2}}{15}+G_{2}(x/z)}\,,$
(4.12)

where $J(x/z)$ and $Y(x/z)$ describe the ideal-gas behaviour of the QED
plasma, while the $G_{1,2}(\tau)$ account for deviations of the QED plasma
from an idea gas. Explicit forms for these expression to ${\cal 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 $\bar{\rho}_{\nu}$, which can generally be
constructed from the neutrino density matrices via

$\frac{{\rm d}\bar{\rho}_{\nu}}{{\rm
d}x}=\frac{1}{\pi^{2}}\int\mathrm{d}y~{}y^{3}\sum_{\alpha}\derivative{\\{\varrho(x,y)\\}_{\alpha\alpha}}{x},$
(4.13)

where ${\rm d}{\\{\varrho(x,y)\\}_{\alpha\alpha}}/{\rm d}x$ corresponds in
general to the QKEs (2.3). Neglecting neutrino oscillations, the density
matrices $\varrho$ 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

$\displaystyle\derivative{\\{\varrho(x,y)\\}_{\alpha\alpha}}{x}$
$\displaystyle\simeq-\frac{\Gamma_{\alpha}(x,y)}{xH(x)}\Big{[}\\{\varrho(x,y)\\}_{\alpha\alpha}-f_{F}((y/z)T)\Big{]}.$
(4.14)

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\in[0.01,30]$). We call this the “full-momentum” approach.
Alternatively, we can simplify the equation (4.14) further by adopting the
ansatz

$\\{\varrho(y)\\}_{\alpha\alpha}=\\{\varrho(\langle
y\rangle)\\}_{\alpha\alpha}\frac{f_{F}(yT)}{f_{F}(\langle y\rangle T)},$
(4.15)

and $\Gamma_{\alpha}(y)=\Gamma_{\alpha}(\langle y\rangle)$. Identifying
$\langle y\rangle$ with the mean momentum mode $y_{0}=3.15z(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

$\frac{{\rm d}\bar{\rho}_{\nu}}{{\rm
d}x}=-\frac{7\pi^{2}}{120xH(x)}\left[\sum_{\alpha}\Gamma_{\alpha}(\langle
y\rangle)\left(\frac{\\{\varrho(x,\langle
y\rangle)\\}_{\alpha\alpha}}{f_{F}(\langle y\rangle
T)}-z^{3}(x)\right)\right].$ (4.16)

We call this alternative the “mean-momentum” approach.

Irrespective of whether we use the full-momentum or the mean-momentum
approach. The final $N_{\rm eff}$ can be estimated from the solutions to
$\rho_{\alpha}$ and $z$ at $x\to\infty$ using the definition of $N^{\rm
SM}_{\rm eff}$ in equation (2.1), reproduced here in terms of the rescaled
variables:

$N_{\rm
eff}=\left.\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\frac{30}{2\pi^{2}}\,\left[\frac{z(x=0)}{z(x)}\right]^{4}\bar{\rho}_{\nu}(x)\right|_{x\to\infty}.$
(4.17)

Table 2 shows our estimates of $\delta N_{\rm eff}/N_{\rm eff}^{\rm LO}$ 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}\neq 0$, $\delta N_{\rm eff}/N_{\rm eff}^{\rm
LO}\simeq(-7.8\to-7.9)\times 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 $\sim 50\%$ between the full-momentum ($\delta N_{\rm eff}^{m_{e}=0}/N_{\rm
eff}^{\rm LO}\simeq-2.6\times 10^{-5}$) and the mean-momentum ($\delta N_{\rm
eff}^{m_{e}=0}/N_{\rm eff}^{\rm LO}\simeq-1.6\times 10^{-5}$) approaches, and
are furthermore $\sim 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\to 0$ relative the LO rate, whereas its $m_{e}\neq
0$ counterpart tends to zero. Since neutrino decoupling in the early universe
is not instantaneous but extends into the $e^{\pm}$-annihilation epoch at
$T\sim m_{e}$, any estimate of $\delta N_{\rm eff}$ that accounts for non-
instantaneous decoupling will to some extent be sensitive to exactly what we
assume for $m_{e}$ in the $T\lesssim 3m_{e}$ region. Indeed, the low-$T$
effect of $m_{e}$ on $\delta N_{\rm eff}/N_{\rm eff}^{\rm LO}$ 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}\neq 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 $\delta N_{\rm eff}/N_{\rm eff}^{\rm LO}$ on what we
assume for the electron mass. Thus, while the overall effect of the rate
corrections on $N_{\rm eff}^{\rm SM}$ is small, we conclude that neglecting
$m_{e}$ in their computations is strictly not a good approximation in high-
precision calculations of $N_{\rm eff}^{\rm SM}$.

##  5 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_{\rm eff}^{\rm SM}$. 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_{\rm eff}^{\rm SM}$: 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\sim
1\to 3$ MeV, the corrections to electron neutrino interaction rate fall in the
range $-0.2\to+0.1\%$ relative to the LO rate, while for the muon and tau
neutrinos the effect is even more minute: $-0.0005\to+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 $\nu_{e}$ and
the $\nu_{\mu,\tau}$ 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_{\rm eff}^{\rm SM}$ 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_{\rm eff}^{\rm SM}$ to fall in the range $\delta N_{\rm
eff}/N_{\rm eff}^{\rm LO}\simeq(-0.5\to-1.1)\times 10^{-5}$. That is, relative
to the current SM benchmark of $N_{\rm eff}^{\rm SM}=3.0440\pm 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 $\delta N_{\rm 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 ${\cal O}(1)$ impact on $\delta N_{\rm eff}/N_{\rm
eff}^{\rm LO}$, even though the rate corrections themselves at $T\sim 1\to 3$
MeV differ by less than $10\%$. This is because corrections assuming $m_{e}=0$
deviate from their $m_{e}\neq 0$ counterparts at $T\lesssim 3m_{e}$, and
diverge relative to LO rates just as the $m_{e}\neq 0$ corrections vanish in
the $T\to 0$ limit. Since neutrino decoupling in the early universe is not
instantaneous, these $T\lesssim 3m_{e}$ effects will imprint on $N_{\rm
eff}^{\rm SM}$ despite the common understanding that neutrino decoupling
occurs at $T\sim 1$ MeV. Thus, while the overall contribution from QED weak
rate corrections to $N_{\rm eff}^{\rm SM}$ is ultimately small, we argue that
neglecting $m_{e}$ their computations is strictly speaking not a good
approximation in high-precision $N_{\rm eff}^{\rm SM}$ calculations.

In conclusion, our results strongly suggest that the SM benchmark value
$N_{\rm eff}^{\rm SM}=3.0440\pm 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_{\rm eff}^{\rm SM}$ that is large enough to be of any
relevance for cosmological observations in the foreseeable future.

## Acknowledgments

MaD and Y3W would like to thank Isabel Oldengott and Jamie McDonald for
discussions. YG acknowledges the support of the French Community of Belgium
through the FRIA grant No. 1.E.063.22F. The research of MK and LPW was
supported by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) through the Research Training Group 2149 “Strong and Weak
Interactions - from Hadrons to Dark matter”. Y3W is supported in part by the
Australian Research Council’s Future Fellowship (project FT180100031) and
Discovery Project (project DP240103130). Computational resources have been
provided by the University of Münster through the HPC cluster PALMA II funded
by the DFG (INST 211/667-1), and the supercomputing facilities of the
Université Catholique de Louvain (CISM/UCL) and the Consortium des Équipements
de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI) funded by the Fond
de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under convention
2.5020.11 and by the Walloon Region.

##  Appendix A Bosonic and fermionic propagators at finite temperatures

For vanishing chemical potentials the free thermal propagators for fermions
are given by

$\displaystyle iS_{\psi}^{11}$
$\displaystyle=\frac{i(\not{p}+m_{\psi})}{p^{2}-m_{\psi}^{2}+i\epsilon}-2\pi(\not{p}+m_{\psi})\delta(p^{2}-m_{\psi}^{2})f_{F}(|p^{0}|)=\gamma^{0}(iS_{\psi}^{22})^{\dagger}\gamma^{0},$
(A.1) $\displaystyle iS_{\psi}^{12}$
$\displaystyle=-2\pi(\not{p}+m_{\psi})\delta(p^{2}-m_{\psi}^{2})\left[f_{F}(|p^{0}|)-\Theta(-p^{0})\right],$
$\displaystyle iS_{\psi}^{21}$
$\displaystyle=-2\pi(\not{p}+m_{\psi})\delta(p^{2}-m_{\psi}^{2})\left[f_{F}(|p^{0}|)-\Theta(p^{0})\right],$

while the ones for bosons are

$\displaystyle i\Delta^{11}$
$\displaystyle=\frac{i}{p^{2}-m^{2}+i\epsilon}+2\pi\delta(p^{2}-m^{2})f_{B}(|p^{0}|)=(i\Delta^{22})^{\ast},$
(A.2) $\displaystyle i\Delta^{12}$
$\displaystyle=2\pi\delta(p^{2}-m^{2})\left[f_{B}(|p^{0}|)+\Theta(-p^{0})\right],$
$\displaystyle i\Delta^{21}$
$\displaystyle=2\pi\delta(p^{2}-m^{2})\left[f_{B}(|p^{0}|)+\Theta(p^{0})\right],$

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,

$\displaystyle i\Delta^{R/A}=\frac{i}{p^{2}-m^{2}\pm
i\mathrm{sgn}(p^{0})\epsilon},$ (A.3)

as well as the statistical propagator

$i\Delta^{+}\equiv\Delta^{12}+\Delta^{21}=-2\pi
i\left[1+2f_{B}(|p^{0}|)\right]\delta(p^{2}-m^{2})=i\left[1+2f_{B}(p^{0})\right]\Delta^{-},$
(A.4)

where the relation to the spectral function
$\Delta^{-}=i(\Delta^{21}-\Delta^{12})$ is another way of stating the KMS
relation. In terms of the scalar propagators, the tree-level photon propagator
is given by

$iD^{ab}_{\mu\nu}(p)=id_{\mu\nu}\Delta^{ab}(p)|_{m=0},$ (A.5)

where we have defined the polarisation tensor,

$d^{\mu\nu}(p)=-g^{\mu\nu}+(1-\xi)\frac{p^{\mu}p^{\nu}}{p^{2}}=-P_{L}^{\mu\nu}-P_{T}^{\mu\nu}-\xi\frac{p^{\mu}p^{\nu}}{p^{2}},$
(A.6)

for a general $R_{\xi}$-gauge.

##  Appendix B Thermal integrals

Throughout the calculations we encounter thermal integrals of the following
forms

$\displaystyle
K^{\mu}(P)=2e^{2}\int\frac{\differential[4]{k}}{(2\pi)^{3}}k^{\mu}\frac{\delta(k^{2}-m_{e}^{2})f_{F}(|k^{0}|)}{(P+k)^{2}-m_{e}^{2}+i\epsilon},$
(B.1) $\displaystyle
K^{\mu\nu}(P)=2e^{2}\int\frac{\differential[4]{k}}{(2\pi)^{3}}k^{\mu}k^{\nu}\frac{\delta(k^{2}-m_{e}^{2})f_{F}(|k^{0}|)}{(P+k)^{2}-m_{e}^{2}+i\epsilon}.$

To evaluate them, we choose $\mathbf{P}$ as the reference direction with the
corresponding unit vector $\hat{\mathbf{P}}$, and parameterise the integration
variable via

$\mathbf{k}=|\mathbf{k}|(\cos\vartheta\hat{\mathbf{P}}+\cos\varphi\sin\vartheta\mathbf{e}_{1}+\sin\varphi\sin\vartheta\mathbf{e}_{2}),$
(B.2)

where, together, $\hat{\mathbf{P}}$ and the two unit vectors $\mathbf{e}_{1}$
and $\mathbf{e}_{2}$ form a complete orthogonal basis.

The contributions from $\mathbf{e}_{1}$ and $\mathbf{e}_{2}$ to the spatial
components of $K^{\mu}$ vanish after the azimuthal integration, giving

$\mathbf{K}=\frac{\mathbf{P}}{|\mathbf{P}|^{3}}\frac{\alpha_{\rm
em}}{4\pi}\int_{m_{e}}^{\infty}\differential{\omega}\left[P^{2}\ell_{1}(\omega,P)+2\omega
P_{0}\ell_{2}(\omega,P)-8|\mathbf{k}||\mathbf{P}|\right]f_{F}(\omega),$ (B.3)

where $\omega=|k^{0}|=\sqrt{|\mathbf{k}|^{2}+m_{e}^{2}}$, and

$\displaystyle\ell_{1}(\omega,P)=\ln\left|\frac{(P^{2}+2|\mathbf{k}||\mathbf{P}|)^{2}-4P_{0}^{2}\omega^{2}}{(P^{2}-2|\mathbf{k}||\mathbf{P}|)^{2}-4P_{0}^{2}\omega^{2}}\right|,$
(B.4)
$\displaystyle\ell_{2}(\omega,P)=\ln\left|\frac{P^{4}-4(P_{0}\omega+|\mathbf{P}||\mathbf{k}|)^{2}}{P^{4}-4(P_{0}\omega-|\mathbf{P}||\mathbf{k}|)^{2}}\right|$

arise from integration over $\cos\vartheta$. Correspondingly,

$K^{0}=\frac{\alpha_{\rm
em}}{2\pi|\mathbf{P}|}\int_{m_{e}}^{\infty}\differential{\omega}\omega\,\ell_{2}(\omega,P)f_{F}(\omega)$
(B.5)

is the remaining, temporal component of $K^{\mu}$.

For $K^{\mu\nu}$, after application of the completeness relation
$\delta^{ij}=\hat{\mathbf{P}}^{i}\hat{\mathbf{P}}^{j}+\sum_{k=1,2}e_{k}^{i}e_{k}^{j}$,
the spatial components become

$K^{ij}=\hat{\mathbf{P}}^{i}\hat{\mathbf{P}}^{j}G_{1}+(\delta^{ij}-\hat{\mathbf{P}}^{i}\hat{\mathbf{P}}^{j})(G_{2}-G_{1}/2),$
(B.6)

where we have introduced the two auxiliary functions

$\displaystyle G_{1}=\frac{\alpha_{\rm
em}}{8\pi|\mathbf{P}|^{3}}\int_{m_{e}}^{\infty}\differential{\omega}\left[(P^{4}+4P_{0}^{2}\omega^{2})\ell_{1}(\omega,P)+4P^{2}P_{0}\,\omega\,\ell_{2}(\omega,P)-8|\mathbf{P}||\mathbf{k}|P^{2}\right]f_{F}(\omega),$
(B.7) $\displaystyle G_{2}=\frac{\alpha_{\rm
em}}{4\pi|\mathbf{P}|}\int_{m_{e}}^{\infty}\differential{\omega}|\mathbf{k}|^{2}\ell_{1}(\omega,P)f_{F}(\omega).$

The remaining components with at least one index being zero are

$\displaystyle K^{00}(P)=\frac{\alpha_{\rm
em}}{2\pi|\mathbf{P}|}\int_{m_{e}}^{\infty}\differential{\omega}\omega^{2}\,\ell_{1}(\omega,P)f_{F}(\omega),$
(B.8) $\displaystyle K^{0i}=\frac{\alpha_{\rm
em}\hat{\mathbf{P}}^{i}}{4\pi|\mathbf{P}|^{2}}\int_{m_{e}}^{\infty}\differential{\omega}\omega\left[2P_{0}\omega\ell_{1}(\omega,P)+P^{2}\ell_{2}(\omega,P)\right]f_{F}(\omega).$

In all cases, the remaining single integral over $\omega$ can, for our
purposes, be evaluated numerically efficiently.

##  Appendix 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 non-zero
electron chemical potentials), no explicit formulae were provided.

###  C.1 1-loop photon self-energy

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

At the one-loop level, the thermal photon self-energy, displayed in figure 5,
reads121212We only compute here the electron contribution to the photon self-
energy, which is by far the dominant contribution at the relevant temperatures
for $N_{\mathrm{eff}}$. It is nonetheless straightforward to extend our
results to higher temperatures where different charged particles can also
contribute.

$\Pi^{\mu\nu}_{ab}(P)=(-1)^{a+b}ie^{2}\int\frac{\differential[4]{k}}{(2\pi)^{4}}\mathrm{Tr}\left[iS^{ab}_{e}(k)\gamma^{\mu}iS_{e}^{ba}(k+P)\gamma^{\nu}\right]$
(C.1)

for general Keldysh indices $a,b\in\\{1,2\\}$. It can always be split into a
longitudinal and a transverse part,

$\Pi_{ab}^{\mu\nu}=\Pi^{L}_{ab}P_{L}^{\mu\nu}+\Pi^{T}_{ab}P_{T}^{\mu\nu},$
(C.2)

with the corresponding projectors [36, 56]

$\displaystyle
P_{T}^{00}=P_{T}^{0i}=P_{T}^{i0}=0,~{}~{}P_{T}^{ij}=-\delta^{ij}+\hat{\mathbf{P}}^{i}\hat{\mathbf{P}}^{j},$
(C.3) $\displaystyle
P_{L}^{\mu\nu}=g^{\mu\nu}-\frac{P^{\mu}P^{\nu}}{P^{2}}-P_{T}^{\mu\nu}.$

Writing the self-energy in this form makes also manifest that
$\Pi^{\mu\nu}_{ab}$ fulfills the the Ward-Takahashi identity
$P_{\mu}\Pi^{\mu\nu}_{ab}=0$. The longitudinal projector can alternatively be
expressed through the heat bath four-velocity $U^{\mu}=\delta^{\mu}_{0}$ (in
its rest frame) as
$P_{L}^{\mu\nu}=\tilde{U}^{\mu}\tilde{U}^{\nu}/{\tilde{U}^{2}}$ with
$\tilde{U}^{\mu}=P^{2}U^{\mu}-(U\cdot P)P^{\mu}$, so that the transverse
projector becomes
$P_{T}^{\mu\nu}=g^{\mu\nu}-\frac{P^{\mu}P^{\nu}}{P^{2}}-P_{L}^{\mu\nu}$.

In the following, we only highlight the derivation of $\Pi_{11}$ and
$\Pi_{12}$, as the remaining self-energies can be related to the first two via
the KMS relation

$\displaystyle\Pi_{21}(P)=e^{P_{0}/T}\Pi_{12}(P)$ (C.4)

and

$\displaystyle\Pi_{22}(P)=-\Pi_{11}(P)^{*}$ (C.5)

for a general four-vector $P$.

Starting with $\Pi^{\mu\nu}_{11}$ ($\Pi^{\mu\nu}_{22}$), its diagonal
component splits into the vacuum self-energy

$\displaystyle\Pi^{\mu\nu}_{2}=ie^{2}\int_{k}\frac{i}{k^{2}-m_{e}^{2}+i\epsilon}\frac{i}{(k+P)^{2}-m_{e}^{2}+i\epsilon}\mathrm{Tr}\left[(\not{k}+m_{e})\gamma^{\mu}(\not{k}+\not{P}+m_{e})\gamma^{\nu}\right]$
(C.6)

and the finite-temperature part $\Pi_{11,T\neq 0}^{\mu\nu}$ ($\Pi_{22,T\neq
0}^{\mu\nu}$). The renormalised vacuum part $\Pi^{\mu\nu}_{2}(P)=(\alpha_{\rm
em}/\pi)(P^{2}g^{\mu\nu}-P^{\nu}P^{\mu})\Pi_{2}(P^{2})$ is textbook material
and, in the $\overline{\mathrm{MS}}$-scheme, reads

$\displaystyle\Pi_{2}(P^{2})$
$\displaystyle=\frac{2}{3}\ln(\frac{\mu}{m_{e}})+\frac{5}{9}+\frac{\tau}{3}+\frac{1}{3}\left(1+\frac{\tau}{2}\right)\sqrt{1-\tau}\ln(\frac{\sqrt{1-\tau}-1}{\sqrt{1-\tau}+1})$
(C.7)
$\displaystyle\longrightarrow\frac{1}{3}\ln(\frac{\mu^{2}}{-P^{2}})+\frac{5}{9}\
\ \text{for}\ m_{e}\to 0,$

with $\tau=4m_{e}^{2}/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^{\mu}$ and $K^{\mu\nu}$ computed in appendix B,
the corresponding real parts of the two diagonal components read

$\displaystyle\mathrm{Re}\,\Pi_{11,T\neq 0}^{\mu\nu}$
$\displaystyle=-\mathrm{Re}\,\Pi_{22,T\neq 0}^{\mu\nu}=4(-g^{\mu\nu}K\cdot
P+2K^{\mu\nu}+K^{\mu}P^{\nu}+K^{\nu}P^{\nu}).$ (C.8)

After collecting all the terms according to the projectors defined in equation
(C.3), we obtain for the longitudinal part

$\displaystyle\mathrm{Re}\,\Pi_{11,T\neq 0}^{L}$
$\displaystyle=\frac{\alpha_{\rm
em}P^{2}}{\pi|\mathbf{P}|^{3}}\int_{m_{e}}^{\infty}\differential{\omega}\left[8k|\mathbf{P}|-\ell_{1}(\omega,P)(P^{2}+4\omega^{2})-4\ell_{2}(\omega,P)P_{0}\omega\right]f_{F}(\omega)$
(C.9)
$\displaystyle\overset{\mathrm{HTL}}{\approx}-3m_{\gamma}^{2}\left(1-\frac{P_{0}^{2}}{|\mathbf{P}|^{2}}\right)\left[1-\frac{P_{0}}{2|\mathbf{P}|}\ln\left|\frac{P_{0}+|\mathbf{P}|\
}{P_{0}-|\mathbf{P}|}\right|\right]\longrightarrow-3m_{\gamma}^{2}\ \
\text{for}\ P_{0}=0,$ (C.10)

and for the transverse one

$\displaystyle\mathrm{Re}\,\Pi_{11,T\neq 0}^{T}$
$\displaystyle=\frac{\alpha_{\rm
em}}{\pi|\mathbf{P}|^{3}}\int_{m_{e}}^{\infty}\differential{\omega}\left[2\ell_{2}(\omega,P)P_{0}P^{2}\omega-4|\mathbf{k}||\mathbf{P}|(P_{0}^{2}+|\mathbf{P}|^{2})\right.$
$\displaystyle\hskip
56.9055pt\left.+\ell_{1}(\omega,P)(|\mathbf{P}|^{2}P^{2}+2P_{0}^{2}\omega^{2}-2|\mathbf{k}|^{2}|\mathbf{P}|^{2}+\frac{1}{2}P^{4})\right]f_{F}(\omega)$
(C.11)
$\displaystyle\overset{\mathrm{HTL}}{\approx}-\frac{3}{2}m_{\gamma}^{2}\frac{P_{0}^{2}}{|\mathbf{P}|^{2}}\left[1-\left(1-\frac{|\mathbf{P}|^{2}}{P_{0}^{2}}\right)\frac{P_{0}}{2|\mathbf{P}|}\ln\left|\frac{P_{0}+|\mathbf{P}|\
}{P_{0}-|\mathbf{P}|}\right|\right]\longrightarrow 0,\ \ \text{for}\ P_{0}=0$
(C.12)

with the thermal photon mass $m_{\gamma}=eT/3$. Here, the symbol
$\overset{\mathrm{HTL}}{\approx}$ 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 $\order{\pi T}$ and soft scales
$\order{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}/|\textbf{k}|$ at
the first step and then in $P_{0}/|\textbf{P}|$ in the second step, using also

$\displaystyle\ell_{1}(\omega,P)\overset{\mathrm{HTL}}{\approx}-2\frac{|\mathbf{P}|}{|\mathbf{k}|},\
\ \
\ell_{2}(\omega,P)\overset{\mathrm{HTL}}{\approx}2\ln\left|\frac{P_{0}+|\mathbf{P}|\
}{P_{0}-|\mathbf{P}|}\right|$ (C.13)

in the same expansion scheme. We further assume that $m_{e}\ll T$, such that
the electron can be treated as effectively massless. This leaves us with an
integral over $\omega$ that can be performed analytically using the relation
$\int_{0}^{\infty}\differential{k}k\,f_{F}(k)=\pi^{2}T^{2}/12$, thereby
yielding the HTL result (C.10) (or (C.12)).

The (purely imaginary) off-diagonal component $\Pi^{\mu\nu}_{12}$, on the
other hand, reads

$\displaystyle\Pi^{\mu\nu}_{12}=i(ie)^{2}\int$
$\displaystyle\frac{\differential[4]{k}}{(2\pi)^{4}}16\pi^{2}\delta\left(k^{2}-m_{e}^{2}\right)(f_{F}(|k^{0}|)-\Theta(k^{0}))\delta\left((k+P)^{2}-m_{e}^{2}\right)$
(C.14)
$\displaystyle\times(f_{F}(|k^{0}+P^{0}|)-\Theta(-k^{0}-P^{0}))\,T^{\mu\nu},$

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

$T^{\mu\nu}=-g^{\mu\nu}k\cdot
P+2k^{\mu}k^{\nu}+k^{\nu}P^{\mu}+k^{\mu}P^{\nu}.$ (C.15)

For the evaluation of the integrals, we make use again of the decomposition of
$\mathbf{k}$ shown in equation (B.2). The additional $\delta$-function here
(compared to the diagonal case) fixes then the polar angle $\theta$ to the
value

$\cos\theta^{\ast}_{\pm}=\frac{P^{2}\pm 2\omega
P_{0}}{2|\mathbf{k}||\mathbf{P}|}.$ (C.16)

With that, $T^{\mu\nu}$ separates into a transverse and a longitudinal part
according to

$\displaystyle
T_{\pm}^{\mu\nu}=\frac{P^{2}}{2}(P_{T}^{\mu\nu}+P_{L}^{\mu\nu})-|\mathbf{k}|^{2}\sin^{2}\theta^{\ast}_{\pm}P_{T}^{\mu\nu}-\frac{P^{2}}{2|\mathbf{P}|^{2}}(P_{0}\pm
2\omega)^{2}P_{L}^{\mu\nu}.$ (C.17)

The final result then reads

$\Pi^{\mu\nu}_{12}=-i\frac{2\alpha_{\rm
em}}{|\mathbf{P}|}\int_{m_{e}}^{\infty}\differential{\omega}\sum_{\pm}\Theta(1-|\cos\theta^{\ast}_{\pm}|)(f_{F}(\omega)-\Theta(k^{0}))(f_{F}(|k^{0}+P^{0}|)-\Theta(-k^{0}-P^{0}))T_{\pm}^{\mu\nu},$
(C.18)

where the sum runs over positive and negative energies $k^{0}=\pm\omega$. In
the HTL limit we find

$\displaystyle\Pi^{L}_{12}\overset{\mathrm{HTL}}{\approx}-i\frac{\pi
e^{2}P^{2}T^{3}}{3|\mathbf{P}|^{3}}\Theta(|\mathbf{P}|-|P^{0}|),$ (C.19)
$\displaystyle\Pi^{T}_{12}\overset{\mathrm{HTL}}{\approx}i\frac{\pi
e^{2}P^{2}T^{3}}{6|\mathbf{P}|^{3}}\Theta(|\mathbf{P}|-|P^{0}|),$ (C.20)

where we have used the integral
$\int_{0}^{\infty}\differential{k}f_{F}(k)\left[1-f_{F}(k)\right]=\pi^{2}T^{3}/6$.

As will be explained in the subsequent section, it is convenient to compute
the retarded and advanced photon self-energy $\Pi_{R/A}$ in order to perform
the resummation of the photon propagator. Using the relation
$\mathrm{Im}\,\Pi_{11}=\frac{i}{2}(\Pi_{12}+\Pi_{21})$, we can write the
advanced and retarded self-energies as

$\displaystyle\Pi^{T/L}_{R/A}=\Pi^{T/L}_{11}+\Pi^{T/L}_{12/21}$
$\displaystyle=\mathrm{Re}\,\Pi^{T/L}_{11}\pm\frac{1}{2}\left(\Pi^{T/L}_{12}-\Pi^{T/L}_{21}\right)$
(C.21)
$\displaystyle=\mathrm{Re}\,\Pi^{T/L}_{11}\pm\frac{1}{2}\left(1-e^{P_{0}/T}\right)\Pi^{T/L}_{12},$

where we have used the KMS relation (C.4) at the second equality. In
particular, equation (C.21) implies that
$\mathrm{Im}\Pi^{R}=-\mathrm{Im}\Pi^{A}$. With this in mind, we can directly
compare our results for the advanced and the retarded transverse photon self-
energies in the HTL limit to the results of Carrington et al. [54]. Given that
$\Pi^{L/T}_{R,A}\overset{\mathrm{HTL}}{\approx}\mathrm{Re}\,\Pi^{T}_{11}\mp\frac{P^{0}}{2T}\Pi^{T}_{12}$
in the HTL limit, our results are in agreement with [54] after carefully
sending the time-ordering parameter $\epsilon$ to zero. In the longitudinal
part we differ by a term $P_{0}^{2}/|\mathbf{P}|^{2}$ with respect to
[54],131313This difference can be seen by comparing figure 7 to figure 4 in
[53]. but agree with [55]. We note in passing that the self-energy (C.18)
leads to the unphysical process of photon decay $\gamma\to e^{+}e^{-}$ at high
enough temperatures, where $m_{\gamma}$ 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.

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 $T=2m_{e}$ and for various choices of
$|\mathbf{P}|/T$.

Figure 7: 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 for different temperatures
$T\in\\{m_{e},2m_{e},5m_{e}\\}$ at a fixed $|\mathbf{P}|/T=0.5$.

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 $|\mathbf{P}|/T$. We compare the exact
one-loop photon self-energy to (i) the equivalent self-energy in the limit
$m_{e}\rightarrow 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}\sim 5$, remains sizeable for
temperatures around the electron neutrino decoupling temperature where
$T/m_{e}\sim 2$.

Beyond the photon self-energy, other quantities of interest are the residue of
the transverse and longitudinal photon141414Longitudinal photons are also
sometimes referred to as “plasmons” in the literature. propagators which are
quantified by

$\displaystyle Z^{T/L}=\lim_{P_{0}\rightarrow
w(|\mathbf{P}|)}\left[1+\frac{\partial\mathrm{Re}\,\Pi^{T/L}_{11}}{\partial
P_{0}^{2}}\right]^{-1},$ (C.22)

where $P_{0}=w(|\mathbf{P}|)$ defines the photon dispersion relation, i.e.,
solutions of the equation $P^{2}+\mathrm{Re}\,\Pi^{T/L}_{11}(P)=0$, which we
show for the benchmark $T=2m_{e}$ case in figure 8. The two residues take the
form

$\displaystyle Z^{L}=$ $\displaystyle\Bigg{[}1-\frac{\alpha_{\rm
em}}{\pi|\mathbf{P}|^{3}}\int_{m_{e}}^{\infty}\mathrm{d}\omega
f_{F}(\omega)\Big{[}\ell^{\prime}_{2}(\omega,P)\left(4\omega
P_{0}P^{2}\right)+\ell_{2}(\omega,P)\left(2\omega(3P_{0}-|\mathbf{P}|^{2}/P_{0})\right)$
(C.23)
$\displaystyle+\ell^{\prime}_{1}(\omega,P)\left(P^{4}+4\omega^{2}P^{2}\right)+\ell_{1}(\omega,P)\left(2P^{2}+4\omega^{2}\right)-8|\mathbf{k}||\mathbf{P}|\Big{]}_{\mid
P_{0}=w(|\mathbf{P}|)}\Bigg{]}^{-1},$

and

$\displaystyle Z^{T}=$ $\displaystyle\Bigg{[}1+\frac{\alpha_{\rm
em}}{2\pi|\mathbf{P}|^{3}}\int_{m_{e}}^{\infty}\mathrm{d}\omega
f_{F}(\omega)\Big{[}\ell^{\prime}_{2}(\omega,P)\left(4\omega
P_{0}P^{2}\right)+\ell_{2}(\omega,P)\left(2\omega(3P_{0}-|\mathbf{P}|^{2}/P_{0})\right)$
(C.24)
$\displaystyle+\ell^{\prime}_{1}(\omega,P)\left(P^{4}+4\omega^{2}P_{0}^{2}-4|\mathbf{k}|^{2}|\mathbf{P}|^{2}+2P^{2}|\mathbf{P}|^{2}\right)$
$\displaystyle+\ell_{1}(\omega,P)\left(4\omega^{2}+2P_{0}^{2}\right)-8|\mathbf{k}||\mathbf{P}|\Big{]}_{\mid
P_{0}=w(|\mathbf{P}|)}\Bigg{]}^{-1},$

where we have defined the following derivatives

$\displaystyle\ell^{\prime}_{1}(\omega,P)\equiv\frac{\partial\ell_{1}}{\partial
P_{0}^{2}}$
$\displaystyle=\frac{2(P^{2}+2|\mathbf{k}||\mathbf{P}|)-4\omega^{2}}{(P^{2}+2|\mathbf{k}||\mathbf{P}|)^{2}-4P_{0}^{2}\omega^{2}}-\frac{2(P^{2}-2|\mathbf{k}||\mathbf{P}|)-4\omega^{2}}{(P^{2}-2|\mathbf{k}||\mathbf{P}|)^{2}-4P_{0}^{2}\omega^{2}},$
(C.25)
$\displaystyle\ell^{\prime}_{2}(\omega,P)\equiv\frac{\partial\ell_{2}}{\partial
P_{0}^{2}}$
$\displaystyle=\frac{2P^{2}-8\omega(P_{0}\omega+|\mathbf{P}|\mathbf{k}|)}{P^{4}-4(P_{0}\omega+|\mathbf{P}||\mathbf{k}|)^{2}}-\frac{2P^{2}-8\omega(P_{0}\omega-|\mathbf{P}||\mathbf{k}|)}{P^{4}-4(P_{0}\omega-|\mathbf{P}||\mathbf{k}|)^{2}},$

of the logarithmic functions $\ell_{1}(\omega,P)$ and $\ell_{2}(\omega,P)$ of
equation (B.4).

Figure 8: Dispersion relation for transverse (blue lines) and longitudinal
(orange lines) photons at the one-loop level with (continuous lines) and
without (dashed lines) the HTL approximation at $T=2m_{e}$. The black dotted
line indicates the dispersion relation at zero temperature, i.e.,
$P_{0}^{2}=|\mathbf{P}|^{2}$.

###  C.2 The resummed photon propagator

The resummed photon propagator $\bar{D}^{ab}_{\mu\nu}$ is given by the Dyson-
Schwinger equation

$\bar{D}^{ab}_{\mu\nu}=D^{ab}_{\mu\nu}+\sum_{c,d}D^{ac}_{\mu\alpha}(\Pi^{\alpha\beta})^{cd}\bar{D}^{db}_{\beta\nu}.$
(C.26)

The sum over the real-time labels $c$ and $d$ can be avoided [59] by re-
expressing the equation in terms of the retarded and advanced propagators

$\displaystyle G^{R}=G_{11}-G_{12}=G_{21}-G_{22},$ (C.27) $\displaystyle
G^{A}=G_{11}-G_{21}=G_{12}-G_{22},$

defined in a similar fashion to the retarded/advanced self-energies (C.21).
This leads to the compact result

$\bar{D}^{R/A}_{\mu\nu}=D^{R/A}_{\mu\nu}+D^{R/A}_{\mu\alpha}(\Pi^{\alpha\beta})^{R/A}\bar{D}^{R/A}_{\beta\nu}.$
(C.28)

Decomposing the resummed photon propagator using the projectors from equations
(C.3),

$\bar{D}^{ab}_{\mu\nu}=d_{\mu\nu}\bar{D}^{ab}=-\bar{D}_{L}^{ab}P^{L}_{\mu\nu}-\bar{D}_{T}^{ab}P^{T}_{\mu\nu}-\xi\Delta^{ab}\frac{P_{\mu}P_{\nu}}{P^{2}},$
(C.29)

and using equation (C.28), the corresponding longitudinal and transverse
components of the resummed retarded and advanced propagator take the form

$\bar{D}_{L/T}^{R/A}=\frac{1}{(\Delta^{R/A})^{-1}+\Pi^{R/A}_{L/T}}.$ (C.30)

With $D_{11}=\frac{1}{2}(D^{A}+D^{R}+D^{+})$ and using the fact that _in
equilibrium_ , the thermal resummed statistical photon propagator
$\bar{D}_{\mu\nu}^{+}$ is related to the retarded and advanced propagators via
[54]

$\bar{D}_{\mu\nu}^{+}(P)=[1+2f_{B}(|P^{0}|)]\mathrm{sgn}(P^{0})\left[\bar{D}^{R}_{\mu\nu}-\bar{D}_{\mu\nu}^{A}\right],$
(C.31)

we finally arrive at

$\displaystyle\bar{D}^{L/T}_{11}=\frac{\Omega_{L,T}}{\Omega^{2}_{L,T}+\Gamma^{2}_{L,T}}-i[1+2f_{B}(|p^{0}|)]\mathrm{sgn}(p^{0})\frac{\Gamma_{L,T}}{\Omega^{2}_{L,T}+\Gamma^{2}_{L,T}},$
(C.32)

with

$\displaystyle\Omega^{L,T}$ $\displaystyle=p^{2}+\mathrm{Re}\,\Pi^{L/T}_{R},$
(C.33) $\displaystyle\Gamma^{L,T}$
$\displaystyle=\mathrm{Im}\,\Pi^{L,T}_{R}=-\frac{1}{2f_{B}(P^{0})}\mathrm{Im}\Pi^{L,T}_{12}.$
(C.34)

It is understood that (C.34) should be replaced by the causality-respecting
$\epsilon$-prescription if the imaginary part of the photon self-energy
vanishes for kinematic reasons at the given order in perturbation theory.

##  Appendix D Parameterisation of the collision integrals

We parameterise the momenta in the self-energy expression (3.8) as in
reference [21],

$\displaystyle\mathbf{p}=|\mathbf{p}|(0,0,1)^{T},$ (D.1)
$\displaystyle\mathbf{l}=|\mathbf{l}|(0,\sin\alpha,\cos\alpha)^{T},$
$\displaystyle\mathbf{q}=|\mathbf{q}|(\sin\theta\sin\beta,\sin\theta\cos\beta,\cos\theta)^{T},$

resulting in

$\displaystyle\mathrm{Tr}\left[\not{p}\Sigma^{12}_{\mathrm{NLO}}\right]=\int\differential[4]{q}\differential[4]{l}\delta(l^{2})\delta(q^{2}-m_{e}^{2})\delta((q+p-l)^{2}-m_{e}^{2})f(q,l,p)$
(D.2)
$\displaystyle=\left.\int\frac{\differential[3]{q}\differential[3]{l}}{4|\mathbf{l}|E_{q}}\sum_{\epsilon,\tau=\pm
1}f(q,l,p)\delta((q+p-l)^{2}-m_{e}^{2})\right|_{\overset{q^{0}=\epsilon
E_{q}}{l^{0}=\tau|l|}}$
$\displaystyle=\frac{2\pi}{4}\left.\int\differential{|\mathbf{l}|}\differential{|\mathbf{q}|}\frac{|\mathbf{q}|^{2}|\mathbf{l}|}{E_{q}}\int_{-1}^{1}\differential{\\!\cos\alpha}\sum_{\epsilon,\tau=\pm
1}\frac{\theta(b^{2}-4ac)}{\sqrt{|a|}}\int_{z_{-}}^{z^{+}}\differential{z}\frac{f(q,l,p)}{\sqrt{(z-z_{-})(z_{+}-z)}}\right|_{\overset{q^{0}=\epsilon
E_{q}}{l^{0}=\tau|\mathbf{l}|}},$

where $E_{q}=\sqrt{|\mathbf{q}|^{2}+m_{e}^{2}}$, and the function $f(q,l,p)$
acts as a placeholder for the integrand in equation (3.8). The last
$\delta$-function fixes the angle $\beta$ to

$\cos\beta_{i}=\frac{|\mathbf{p}||\mathbf{q}|\cos\theta+\epsilon\tau
E_{q}|\mathbf{l}|+\tau|\mathbf{l}||\mathbf{p}|-|\mathbf{p}||\mathbf{l}|\cos\alpha-\epsilon
E_{q}|\mathbf{p}|-|\mathbf{l}||\mathbf{q}|\cos\alpha\cos\theta}{|\mathbf{l}||\mathbf{q}|\sin\alpha\sin\theta}.$
(D.3)

The integration limits on $z=\cos\theta$,

$z_{\pm}=-\frac{b}{2a}\pm\sqrt{\left(\frac{b}{2a}\right)^{2}-\frac{c}{a}},$
(D.4)

are the roots of the quadratic polynomial $az^{2}+bz+c$, with

$\displaystyle a$ $\displaystyle=-q^{2}|\mathbf{p}-\mathbf{l}|^{2},$ (D.5)
$\displaystyle b$
$\displaystyle=-2(\cos\alpha|\mathbf{l}|-|\mathbf{p}|)|\mathbf{q}|(|\mathbf{l}||\mathbf{p}|(\cos\alpha-\tau)+\epsilon
E_{q}(|\mathbf{p}|-|\mathbf{l}|\tau)),$ $\displaystyle c$
$\displaystyle=(E_{q}^{2}(|\mathbf{p}|-|\mathbf{l}|\tau)^{2})+2\epsilon
E_{q}|\mathbf{l}||\mathbf{p}|(\cos\alpha-\tau)(-|\mathbf{p}|+|\mathbf{l}|\tau)$
$\displaystyle\hskip
56.9055pt+|\mathbf{l}|^{2}(-((\cos^{2}\alpha)|\mathbf{p}|^{2})+|\mathbf{q}|^{2}(1-\cos^{2}\alpha)+2\cos\alpha|\mathbf{p}|^{2}\tau-|\mathbf{p}|^{2})$

as the coefficients.

The $z$-integration of the vacuum and the thermal corrections can be performed
analytically following

$I^{z}=\int_{z_{-}}^{z^{+}}\differential{z}\frac{\mathcal{T}(q,l,p)}{\sqrt{(z-z_{-})(z_{+}-z)}}=\mathcal{A}\,\pi\left(\mathcal{G}^{0}-\frac{b}{2a}\mathcal{G}^{1}+\frac{3b^{2}-4ca}{8a^{2}}\mathcal{G}^{2}\right),$
(D.6)

where we have expanded the correction matrix ${\cal T}$ in powers of $z$
according to

$\mathcal{T}=\mathcal{A}\left(\mathcal{G}^{0}+\mathcal{G}^{1}z+\mathcal{G}^{2}z^{2}\right).$
(D.7)

Then, for the transverse thermal part, the coefficients read

$\displaystyle\mathcal{G}^{0}_{\mathrm{th},T}=\frac{P^{2}}{4|\mathbf{P}|^{2}}\Bigl{[}2|\mathbf{l}|^{2}P^{2}+4l_{0}|\mathbf{p}|P_{0}q_{0}+8l_{0}|\mathbf{p}|q_{0}^{2}+2l_{0}P_{0}P^{2}+2m_{e}^{2}|\mathbf{P}|^{2}+8|\mathbf{p}|^{2}P_{0}^{2}q_{0}^{2}/P^{2}$
(D.8) $\displaystyle\ \ \ \ \ \ \ \ \ \
-8|\mathbf{p}|^{2}q_{0}^{2}+4|\mathbf{p}|P_{0}^{2}q_{0}+4|\mathbf{p}|P_{0}q_{0}^{2}-2|\mathbf{p}|P^{2}q_{0}+P_{0}^{2}P^{2}+2P_{0}P^{2}q_{0}+2P^{2}q_{0}^{2}\Bigr{]},$
$\displaystyle\mathcal{G}^{1}_{\mathrm{th},T}=-\frac{|\mathbf{p}||\mathbf{q}|P^{2}}{2|\mathbf{P}|^{2}}\left[2l_{0}(P_{0}+2q_{0})+2P_{0}^{2}(1+4|\mathbf{p}|q_{0}/P^{2})-8|\mathbf{p}|q_{0}-P^{2}+2P_{0}q_{0}\right],$
$\displaystyle\mathcal{G}^{2}_{\mathrm{th},T}=2|\mathbf{p}|^{2}|\mathbf{q}|^{2},$

while for the longitudinal thermal part we have

$\displaystyle\mathcal{G}^{0}_{\mathrm{th},L}=-\frac{P^{2}}{8|\mathbf{P}|^{2}}\Bigl{[}4P_{0}q_{0}(P^{2}+2|\mathbf{p}|q_{0})+4P^{2}l^{2}+(P^{2}+4|\mathbf{p}|q_{0})P^{2}_{0}+P^{2}(P^{2}+4q^{2}_{0})$
(D.9) $\displaystyle\ \ \ \ \ \ \ \ \ \
+4l_{0}(P_{0}(P^{2}+2|\mathbf{p}|q_{0})+4|\mathbf{p}|q^{2}_{0})\Bigr{]},$
$\displaystyle\mathcal{G}^{1}_{\mathrm{th},L}=\frac{P^{2}|\mathbf{q}||\mathbf{p}|}{2|\mathbf{P}|^{2}}(l_{0}+|\mathbf{p}|)(P_{0}+2q_{0}),$
$\displaystyle\mathcal{G}^{2}_{\mathrm{th},L}=0,$

with the prefactor $\mathcal{A}^{L/T}_{\rm
th}=-2^{8}G_{F}^{2}(g_{V}^{\alpha})^{2}\mathrm{Re}\,\bar{D}^{L,T}_{11}\mathrm{Re}\,\Pi^{L,T}_{11,T\neq
0}$. For the vacuum part, the corresponding coefficients are

$\displaystyle\mathcal{G}^{0}_{\mathrm{vac}}=-m_{e}^{2}(l\cdot
p)+2|\mathbf{p}|^{2}q_{0}^{2},$ (D.10)
$\displaystyle\mathcal{G}^{1}_{\mathrm{vac}}=-4|\mathbf{q}|q_{0}|\mathbf{p}|^{2},$
$\displaystyle\mathcal{G}^{2}_{\mathrm{vac}}=2|\mathbf{p}|^{2}|\mathbf{q}|^{2},$

with $\mathcal{A}_{\rm vac}=-2^{10}G_{F}^{2}\,\alpha_{\rm
em}\,[(g^{\alpha}_{V})^{2}/(4\pi)]\,\Pi_{2}$.

##  Appendix 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 ${\cal O}(e^{2})$ [20],

$\displaystyle J(\tau)$
$\displaystyle\equiv\frac{1}{\pi^{2}}\int_{0}^{\infty}{\rm
d}\omega\,\omega^{2}\frac{\exp(\sqrt{\omega^{2}+\tau^{2}})}{[\exp(\sqrt{\omega^{2}+\tau^{2}})+1]^{2}},$
(E.1) $\displaystyle Y(\tau)$
$\displaystyle\equiv\frac{1}{\pi^{2}}\int_{0}^{\infty}{\rm
d}\omega\,\omega^{4}\frac{\exp(\sqrt{\omega^{2}+\tau^{2}})}{[\exp(\sqrt{\omega^{2}+\tau^{2}})+1]^{2}},$
$\displaystyle G_{1}(\tau)=$ $\displaystyle 2\pi\alpha_{\rm
em}\Bigg{[}\frac{1}{\tau}\left(\frac{K(\tau)}{3}+2K^{2}(\tau)-\frac{J(\tau)}{6}-K(\tau)J(\tau)\right)$
(E.2)
$\displaystyle+\left(\frac{K^{\prime}(\tau)}{6}-K(\tau)K^{\prime}(\tau)+\frac{J^{\prime}(\tau)}{6}+J^{\prime}(\tau)K(\tau)+J(\tau)K^{\prime}(\tau)\right)\Bigg{]},$
$\displaystyle G_{2}(\tau)=$ $\displaystyle-8\pi\alpha_{\rm
em}\left(\frac{K(\tau)}{6}+\frac{J(\tau)}{6}-\frac{1}{2}K^{2}(\tau)+K(\tau)J(\tau)\right)$
(E.3) $\displaystyle+2\pi\alpha_{\rm
em}\left(\frac{K^{\prime}(\tau)}{6}-K(\tau)K^{\prime}(\tau)+\frac{J^{\prime}(\tau)}{6}+J^{\prime}(\tau)K(\tau)+J(\tau)K^{\prime}(\tau)\right),$

where ′ denotes a derivative with respect to $\tau$. Expressions for
$G_{1,2}(\tau)$ at ${\cal O}(e^{3})$ can be found in reference [21].

## References

  * [1] G. Steigman, D.N. Schramm and J.E. Gunn, _Cosmological Limits to the Number of Massive Leptons_ , _Phys. Lett. B_ 66 (1977) 202.
  * [2] R. Barbieri and A. Dolgov, _Neutrino oscillations in the early universe_ , _Nucl. Phys. B_ 349 (1991) 743.
  * [3] A.M. Abdullahi et al., _The present and future status of heavy neutral leptons_ , _J. Phys. G_ 50 (2023) 020501 [2203.08039].
  * [4] L. Di Luzio, J. Martin Camalich, G. Martinelli, J.A. Oller and G. Piazza, _Axion-pion thermalization rate in unitarized NLO chiral perturbation theory_ , _Phys. Rev. D_ 108 (2023) 035025 [2211.05073].
  * [5] F. D’Eramo, F. Hajkarim and S. Yun, _Thermal QCD Axions across Thresholds_ , _JHEP_ 10 (2021) 224 [2108.05371].
  * [6] C. Caprini and D.G. Figueroa, _Cosmological Backgrounds of Gravitational Waves_ , _Class. Quant. Grav._ 35 (2018) 163001 [1801.04268].
  * [7] A. Aboubrahim, M. Klasen and P. Nath, _Analyzing the Hubble tension through hidden sector dynamics in the early universe_ , _JCAP_ 04 (2022) 042 [2202.04453].
  * [8] P. Agrawal et al., _Feebly-interacting particles: FIPs 2020 workshop report_ , _Eur. Phys. J. C_ 81 (2021) 1015 [2102.12143].
  * [9] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [10] S. Dodelson and M.S. Turner, _Nonequilibrium neutrino statistical mechanics in the expanding universe_ , _Phys. Rev. D_ 46 (1992) 3372.
  * [11] S. Hannestad and J. Madsen, _Neutrino decoupling in the early universe_ , _Phys. Rev. D_ 52 (1995) 1764 [astro-ph/9506015].
  * [12] A.D. Dolgov, S.H. Hansen and D.V. Semikoz, _Nonequilibrium corrections to the spectra of massless neutrinos in the early universe_ , _Nucl. Phys. B_ 503 (1997) 426.
  * [13] A.D. Dolgov, S.H. Hansen and D.V. Semikoz, _Nonequilibrium corrections to the spectra of massless neutrinos in the early universe: Addendum_ , _Nucl. Phys. B_ 543 (1999) 269.
  * [14] S. Esposito, G. Miele, S. Pastor, M. Peloso and O. Pisanti, _Nonequilibrium spectra of degenerate relic neutrinos_ , _Nucl. Phys. B_ 590 (2000) 539.
  * [15] J. Froustey and C. Pitrou, _Incomplete neutrino decoupling effect on big bang nucleosynthesis_ , _Phys. Rev. D_ 101 (2020) 043524 [1912.09378].
  * [16] 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 Temperature Corrections to Weak Rates_ , _Phys. Rev. D_ 26 (1982) 2694.
  * [17] A. Heckler, _Astrophysical applications of quantum corrections to the equation of state of a plasma_ , _Phys. Rev. D_ 49 (1994) 611.
  * [18] N. Fornengo, C. Kim and J. Song, _Finite temperature effects on the neutrino decoupling in the early universe_ , _Phys. Rev. D_ 56 (1997) 5123 [hep-ph/9702324].
  * [19] R.E. Lopez and M.S. Turner, _An Accurate Calculation of the Big Bang Prediction for the Abundance of Primordial Helium_ , _Phys. Rev. D_ 59 (1999) 103502 [astro-ph/9807279].
  * [20] G. Mangano, G. Miele, S. Pastor and M. Peloso, _A Precision calculation of the effective number of cosmological neutrinos_ , _Phys. Lett. B_ 534 (2002) 8 [astro-ph/0111408].
  * [21] J.J. Bennett, G. Buldgen, M. Drewes and Y.Y.Y. Wong, _Towards a precision calculation of the effective number of neutrinos $N_{\rm eff}$ in the Standard Model I: the QED equation of state_, _JCAP_ 03 (2020) 003 [1911.04504].
  * [22] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti and P.D. Serpico, _Relic neutrino decoupling including flavor oscillations_ , _Nucl. Phys. B_ 729 (2005) 221 [hep-ph/0506164].
  * [23] J. Birrell, C.-T. Yang and J. Rafelski, _Relic Neutrino Freeze-out: Dependence on Natural Constants_ , _Nucl. Phys. B_ 890 (2014) 481 [1406.1759].
  * [24] E. Grohs, G.M. Fuller, C.T. Kishimoto, M.W. Paris and A. Vlasenko, _Neutrino energy transport in weak decoupling and big bang nucleosynthesis_ , _Phys. Rev. D_ 93 (2016) 083522 [1512.02205].
  * [25] P.F. de Salas and S. Pastor, _Relic neutrino decoupling with flavour oscillations revisited_ , _JCAP_ 07 (2016) 051 [1606.06986].
  * [26] S. Gariazzo, P.F. de Salas and S. Pastor, _Thermalisation of sterile neutrinos in the early Universe in the 3+1 scheme with full mixing matrix_ , _JCAP_ 07 (2019) 014 [1905.11290].
  * [27] K. Akita and M. Yamaguchi, _A precision calculation of relic neutrino decoupling_ , _JCAP_ 08 (2020) 012 [2005.07047].
  * [28] J. Froustey, C. Pitrou and M.C. Volpe, _Neutrino decoupling including flavour oscillations and primordial nucleosynthesis_ , _JCAP_ 12 (2020) 015 [2008.01074].
  * [29] J.J. Bennett, G. Buldgen, P.F. De Salas, M. Drewes, S. Gariazzo, S. Pastor et al., _Towards a precision calculation of $N_{\rm eff}$ in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED_, _JCAP_ 04 (2021) 073 [2012.02726].
  * [30] M. Cielo, M. Escudero, G. Mangano and O. Pisanti, _Neff in the Standard Model at NLO is 3.043_ , _Phys. Rev. D_ 108 (2023) L121301 [2306.05460].
  * [31] S. Esposito, G. Mangano, G. Miele, I. Picardi and O. Pisanti, _Neutrino energy loss rate in a stellar plasma_ , _Nucl. Phys. B_ 658 (2003) 217 [astro-ph/0301438].
  * [32] G. Jackson and M. Laine, _QED corrections to the thermal neutrino interaction rate_ , 2312.07015.
  * [33] O. Tomalak and R.J. Hill, _Theory of elastic neutrino-electron scattering_ , _Phys. Rev. D_ 101 (2020) 033006 [1907.03379].
  * [34] R.J. Hill and O. Tomalak, _On the effective theory of neutrino-electron and neutrino-quark interactions_ , _Phys. Lett. B_ 805 (2020) 135466 [1911.01493].
  * [35] CMB-S4 collaboration, _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * [36] N.P. Landsman and C.G. van Weert, _Real and Imaginary Time Field Theory at Finite Temperature and Density_ , _Phys. Rept._ 145 (1987) 141.
  * [37] K.-c. Chou, Z.-b. Su, B.-l. Hao and L. Yu, _Equilibrium and Nonequilibrium Formalisms Made Unified_ , _Phys. Rept._ 118 (1985) 1.
  * [38] J. Berges, _Introduction to nonequilibrium quantum field theory_ , _AIP Conf. Proc._ 739 (2004) 3 [hep-ph/0409233].
  * [39] J.I. Kapusta and C. Gale, _Finite-temperature field theory: Principles and applications_ , Cambridge Monographs on Mathematical Physics, Cambridge University Press (2011), 10.1017/CBO9780511535130.
  * [40] C. Coriano and R.R. Parwani, _The Three loop equation of state of QED at high temperature_ , _Phys. Rev. Lett._ 73 (1994) 2398 [hep-ph/9405343].
  * [41] R.R. Parwani and C. Coriano, _Higher order corrections to the equation of state of QED at high temperature_ , _Nucl. Phys. B_ 434 (1995) 56 [hep-ph/9409269].
  * [42] G. Sigl and G. Raffelt, _General kinetic description of relativistic mixed neutrinos_ , _Nucl. Phys. B_ 406 (1993) 423.
  * [43] L. Canetti, M. Drewes and M. Shaposhnikov, _Matter and Antimatter in the Universe_ , _New J. Phys._ 14 (2012) 095012 [1204.4186].
  * [44] S. Antusch, E. Cazzato, M. Drewes, O. Fischer, B. Garbrecht, D. Gueter et al., _Probing Leptogenesis at Future Colliders_ , _JHEP_ 09 (2018) 124 [1710.03744].
  * [45] J.S. Schwinger, _Brownian motion of a quantum oscillator_ , _J. Math. Phys._ 2 (1961) 407.
  * [46] L.V. Keldysh, _Diagram technique for nonequilibrium processes_ , _Zh. Eksp. Teor. Fiz._ 47 (1964) 1515.
  * [47] B. Garbrecht and M. Garny, _Finite Width in out-of-Equilibrium Propagators and Kinetic Theory_ , _Annals Phys._ 327 (2012) 914 [1108.3688].
  * [48] M. Drewes, S. Mendizabal and C. Weniger, _The Boltzmann Equation from Quantum Field Theory_ , _Phys. Lett. B_ 718 (2013) 1119 [1202.1301].
  * [49] E. Braaten and D. Segel, _Neutrino energy loss from the plasma process at all temperatures and densities_ , _Phys. Rev. D_ 48 (1993) 1478 [hep-ph/9302213].
  * [50] E.W. Kolb and M.S. Turner, _The Early Universe_ , vol. 69 (1990), 10.1201/9780429492860.
  * [51] E. Braaten, _Neutrino emissivity of an ultrarelativistic plasma from positron and plasmino annihilation_ , _Astrophys. J._ 392 (1992) 70.
  * [52] Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * [53] A. Peshier, K. Schertler and M.H. Thoma, _One loop selfenergies at finite temperature_ , _Annals Phys._ 266 (1998) 162 [hep-ph/9708434].
  * [54] M.E. Carrington, D.-f. Hou and M.H. Thoma, _Equilibrium and nonequilibrium hard thermal loop resummation in the real time formalism_ , _Eur. Phys. J. C_ 7 (1999) 347 [hep-ph/9708363].
  * [55] M.L. Bellac, _Thermal Field Theory_ , Cambridge Monographs on Mathematical Physics, Cambridge University Press (3, 2011), 10.1017/CBO9780511721700.
  * [56] A. Rebhan, _Hard thermal loops and QCD thermodynamics_ , in _Cargese Summer School on QCD Perspectives on Hot and Dense Matter_ , pp. 327–351, 11, 2001 [hep-ph/0111341].
  * [57] M. Vanderhaeghen, J.M. Friedrich, D. Lhuillier, D. Marchand, L. Van Hoorebeke and J. Van de Wiele, _QED radiative corrections to virtual Compton scattering_ , _Phys. Rev. C_ 62 (2000) 025501 [hep-ph/0001100].
  * [58] E. Braaten and R.D. Pisarski, _Soft Amplitudes in Hot Gauge Theories: A General Analysis_ , _Nucl. Phys. B_ 337 (1990) 569.
  * [59] J. Ghiglieri, A. Kurkela, M. Strickland and A. Vuorinen, _Perturbative Thermal QCD: Formalism and Applications_ , _Phys. Rept._ 880 (2020) 1 [2002.10188].