# QED corrections to the thermal neutrino interaction rate

G Jackson, M Laine

► To cite this version:

G Jackson, M Laine. QED corrections to the thermal neutrino interaction rate. Journal of High Energy Physics, 2024, 05, pp.089. 10.1007/JHEP05(2024)089. hal-04355096

HAL Id: hal-04355096

https://hal.science/hal-04355096v1

Submitted on 4 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RECEIVED: December 20, 2023

REVISED: March 8, 2024

ACCEPTED: March 19, 2024

PUBLISHED: May 8, 2024

# QED corrections to the thermal neutrino interaction rate

G. Jackson ID a,b and M. Laine ID c

E-mail: jackson@subatech.in2p3.fr, laine@itp.unibe.ch

**ABSTRACT:** Motivated by precision computations of neutrino decoupling at MeV temperatures, we show how QED corrections to the thermal neutrino interaction rate can be related to the electron-positron spectral function as well as an effective $\bar{\nu}\nu\gamma$ vertex. The spectral function is needed both in a timelike and in a spacelike domain, and for both of its physical polarization states (transverse and longitudinal with respect to spatial momentum). Incorporating an NLO evaluation of this spectral function, an estimate of the $\bar{\nu}\nu\gamma$ vertex, and HTL resummation of scatterings mediated by soft Bose-enhanced t-channel photons, we compute the interaction rate as a function of the neutrino momentum and flavour. Effects on the –(0...2)% level are found, noticeably smaller than a previous estimate of a related quantity.

**KEYWORDS:** Early Universe Particle Physics, Finite Temperature or Finite Density, Neutrino Interactions

ARXIV ePRINT: 2312.07015

<footnote><sup>a</sup>Institute for Nuclear Theory, University of Washington,
Box 351550, Seattle, WA 98195-1550, U.S.A.</footnote>

<footnote><sup>b</sup>SUBATECH, Nantes Université, IMT Atlantique, IN2P3/CNRS,
4 rue Alfred Kastler, La Chantrerie BP 20722, 44307 Nantes, France</footnote>

<footnote><sup>c</sup>AEC, Institute for Theoretical Physics, University of Bern,
Sidlerstrasse 5, CH-3012 Bern, Switzerland</footnote>
# Contents

1 Introduction 1

2 Technical setup 2

2.1 Observables 3

2.2 Notation for the practical computation 3

2.3 Interactions according to a Fermi effective theory 4

2.4 Contractions for neutrino self-energy 7

2.5 Carrying out a Matsubara sum and angular integrals 8

2.6 Projecting out the interaction rate 9

2.7 HTL resummation of soft t-channel photons 11

3 NLO computation and result 12

4 Conclusions 14

A Evaluation at leading order 17

A.1 Additional simplifications of the leading-order (LO) expression 17

A.2 Numerical crosscheck of the full integral representation 19

A.3 Real part of the vector correlator 20

B Leading-logarithmic part of the NLO correction 20

C Numerical results for coefficient functions 21

## 1 Introduction

The energy density carried by relativistic degrees of freedom at the time of primordial nucleosynthesis or photon decoupling is parametrized by an effective number of neutrino species, $N_{\text{eff}}$. It can be inferred experimentally, and computed theoretically, either in the Standard Model, or in a given extension thereof. In the Standard Model, the naive estimate is 3, but the precisely computed value differs from this. The current best estimate, obtained after many decades of work, reads $N_{\text{eff}}^{(\text{SM})} \approx 3.043 \pm 0.001$ [1–4]. Were the observed value to eventually differ from $N_{\text{eff}}^{(\text{SM})}$, this could be an indication of physics beyond the Standard Model.

The latest ingredient included in the Standard Model prediction of $N_{\text{eff}}$ originates from QED corrections to the rate at which neutrinos interact with the electron-positron plasma at the time of their decoupling, $T \sim (1...3)\text{MeV}$ [4]. This result was based on previous work in a different (astrophysical) environment [5]. Concretely, the observable evaluated is

$$Q = \left. \frac{\delta\rho}{\delta t} \right|_{e^+e^- \to \bar{\nu}\nu(\gamma)} = \langle (\epsilon_{e^+} + \epsilon_{e^-}) \sum |\mathcal{M}|^2_{e^+e^- \to \bar{\nu}\nu(\gamma)} \rangle_{\text{phase space}}, \quad (1.1)$$

where $\rho$ denotes the energy density. The last form stands for the phase-space average of the electron-positron energy loss, weighted by the matrix element squared of this process.
Given that the computation of refs. [4, 5] is quite complicated, and that it entails certain approximations, it seems well motivated to carry out an independent analysis of the rate at which neutrinos interact with the QED plasma. However, there is also the conceptual issue that observables should in principle be defined without reference to Boltzmann equations, which have a limited range of applicability. Therefore, we consider a quantity which is related to eq. (1.1) but *not* equivalent, with the benefit that it can be unambiguously defined beyond leading order. Concretely, this is achieved by considering the imaginary part of the retarded neutrino self-energy, frequently called the thermal neutrino interaction rate. Physically, this quantity can be argued to affect the time evolution of a neutrino density matrix (cf. section 2.1).

Our computation entails some practical differences compared with refs. [4, 5]. First of all, $e^{+}e^{-} \to \bar{\nu}\nu(\gamma)$ represents only a subclass of the processes that we consider (cf. figure 2 on page 13 for the full set); we add reactions like $e\nu \to e\bar{\nu}$ including its 1-loop virtual corrections, or $e\gamma \to e\bar{\nu}\nu$, or a logarithm $\sim \ln(T/m_e)$ that can be associated with a loop-induced $\bar{\nu}\nu\gamma$ vertex. Second, we include Pauli blocking of the final-state neutrinos, and Bose enhancement of both real and virtual photons, and Hard Thermal Loop (HTL) resummation required for properly treating soft virtual photons, which leads to another logarithm $\sim \ln(1/\alpha_{\text{em}})$. Third, eq. (1.1) represents a momentum average, whereas we determine the differential interaction rate, as a function of momentum ($k$) with respect to the heat bath. But there are simplifications on our side as well, notably we assume that the neutrino and electron ensembles carry the same temperature, which is physically the case only at the start of the decoupling process; and we omit the electron mass, which is again best justified at the highest temperatures.

In view of these differences, both in the observable itself and on the technical side, it may not be surprising that our results differ from those in ref. [4] (cf. section 4). Notably, restricting to the so-called s-channel processes, corresponding to those evaluated in ref. [4], we find the opposite sign of the relative NLO correction; a significantly smaller overall magnitude; and a larger dependence on the neutrino flavour.

Our presentation is organized as follows. In section 2, we describe our starting point, the essential methods used, and the main steps of our computation. A reader not interested in such details is advised to skip directly to section 3, where our NLO results are given and elaborated upon. Conclusions are offered in section 4. Appendix A displays computational details related to the leading-order (LO) evaluation, appendix B an HTL-resummed computation of the leading-logarithmic part of the NLO result, and appendix C supplementary plots.

## 2 Technical setup

The computation of this paper is technical in nature, and in this section we specify the steps that it entails. Our tool is thermal field theory in the imaginary-time formalism. Parts of the results may also be obtained from Boltzmann equations, however the latter are not well suited to incorporating temperature-dependent virtual effects, even though at the NLO level virtual effects can be as important as real corrections (and cancel parts of them).
## 2.1 Observables

To describe the non-equilibrium time evolution of neutrinos in the early universe, it is useful to define a density matrix for a momentum mode $k \equiv |k|,$

$$ \rho_{ab} \equiv e^{i\phi_{ab}(t)} w_a^\dagger w_b, \qquad (2.1) $$

where $w_a^\dagger$ and $w_b$ are creation and annihilation operators, respectively, and $a, b$ are flavour indices.<fnref n="1" /> The phase factor $\phi_{ab}$ depends on the time evolution picture chosen. After averaging over the “fast” medium degrees of freedom, the effective evolution equation for the “slow” variables is sometimes assumed to take the form

$$ \dot{\rho} \approx i[\Delta\mathcal{E}_k, \rho] + \frac{1}{2}\{\Gamma_k, \rho_{\text{eq}} - \rho\}. \qquad (2.2) $$

Here the time derivative reads $(\dots) \equiv (\partial_t - Hk\partial_k)(\dots)$, where $H$ is the Hubble rate; $\Delta\mathcal{E}_k$ is a matrix involving energy differences; and $\rho_{\text{eq}}$ is a would-be equilibrium density matrix. Equation (2.2) can be obtained from linear response theory, if the density matrix is perturbed around equilibrium in a single momentum bin $k$. In this language, our goal is to determine the matrix $\Gamma_k$.<fnref n="2" /> We denote its diagonal components (in the interaction basis) by $\Gamma_{k;\nu_a}$.

Actually, the form sketched in eq. (2.2) is oversimplified. Through weak interactions, the neutrinos interact with electrons and positrons through pair production, pair annihilation, and elastic scattering. The temperature of the electron-positron plasma is denoted by $T$. However, the neutrinos also experience interactions with each other. As they are falling out of equilibrium, the effective temperatures of the neutrino ensembles start to differ from $T$. In this situation, the scatterings off neutrinos drive the neutrinos not towards a universal equilibrium distribution, but rather towards the scatterers' ensemble. Having this in mind, we disentangle the full interaction rate into partial ones,

$$ \Gamma_{k;\nu_a} \approx \sum_i \Gamma_{k;\nu_a}^{(i)}, \qquad (2.3) $$

where $i$ labels the different ensembles. We employ this representation in section 3, however assume the presence of a universal temperature in all the ensembles.

## 2.2 Notation for the practical computation

In order to avoid excessive imaginary units and minus signs, we employ Euclidean conventions for the practical computation, with the path integral weight appearing as

$$ e^{iS_M} \rightarrow e^{-S_E}. \qquad (2.4) $$

Minkowskian ($t$) and Euclidean ($\tau$) time coordinates are formally related via the Wick rotation $\tau \leftrightarrow it$, whereas for the corresponding momenta we have $k_n \leftrightarrow -ik^0$, where the Matsubara frequencies read $k_n = 2n\pi T$ for bosons, and $k_n = (2n + 1)\pi T$ for fermions, with $n \in \mathbb{Z}$.

<footnote n="1">¹Often the flavour indices are placed in the opposite order, to have a text-book appearance of the overall sign in the first term of eq. (2.2) [6].</footnote>

<footnote n="2">²Another important effect are medium modifications to $\Delta\mathcal{E}_k$ [7]. They are of $\mathcal{O}(G_F)$, where $G_F$ denotes the Fermi coupling, whereas $\Gamma_k$ is of $\mathcal{O}(G_F^2)$.</footnote>
Whenever in Minkowskian spacetime, the signature (+---) is employed. Minkowskian four-momenta are denoted by $\mathcal{K} = (k^0, \mathbf{k})$, Euclidean ones by $K = (k_n, \mathbf{k})$. We use an implicit notation where the symbol used indicates whether analytic continuation has been carried out, for example

$$ \Sigma_{\mathcal{K}} \equiv \Sigma_K|_{k_n \to -i[k^0 + i0^+]}, \quad \Sigma_k \equiv \Sigma_{\mathcal{K}}|_{k^0 \to k}. \qquad (2.5) $$

Euclidean Dirac matrices, identified by the fact that all indices are down, are defined as

$$ \gamma_0 \equiv \gamma^0, \quad \gamma_i \equiv -i\gamma^i. \qquad (2.6) $$

It follows that

$$ A\gamma_{\mu}B C\gamma_{\mu}D = A\gamma_{\mu}B C\gamma^{\mu}D, \quad i\not\!K = \not\!K. \qquad (2.7) $$

The Fermi coupling reads at tree level, or leading order (LO),

$$ G_F \stackrel{\text{LO}}{=} \frac{1}{\sqrt{2}v^2} \stackrel{\text{LO}}{=} \frac{g_2^2}{4\sqrt{2}m_W^2}, \quad v \approx 246 \text{ GeV}, \qquad (2.8) $$

and the Weinberg angle is parametrized by

$$ s_W^2 \equiv \sin^2 \theta_W \stackrel{\text{LO}}{=} \frac{g_1^2}{g_1^2 + g_2^2}, \qquad (2.9) $$

where $g_1$ and $g_2$ are the gauge couplings of $U_Y(1)$ and $SU_L(2)$, respectively. As we are carrying out a loop computation, the precise definitions of $G_F$ and $s_W^2$ need to be revisited later on. The numerical values used for their $\overline{\text{MS}}$ renormalized versions at the scale $\bar{\mu} = m_e$ are specified in the next section.

## 2.3 Interactions according to a Fermi effective theory

At temperatures $T \sim$ MeV, all hadrons, as well as the $\mu$ and $\tau$ leptons, are heavy, appearing in the thermal bath with an exponentially suppressed weight. Therefore we only need to include neutrinos, electrons, positrons and photons as dynamical degrees of freedom. Working in the interaction basis, we may thus take the Lagrangian

$$ L_E \supset \frac{G_F}{4\sqrt{2}} \left\{ \bar{\nu}_a \gamma_\mu (1-\gamma_5) \nu_a \bar{\nu}_b \gamma_\mu (1-\gamma_5) \nu_b + 4 \bar{\nu}_e \gamma_\mu (1-\gamma_5) \ell_e \bar{\ell}_e \gamma_\mu (1-\gamma_5) \nu_e \right. \\ \left. - 2 \bar{\nu}_a \gamma_\mu (1-\gamma_5) \nu_a \bar{\ell}_e \gamma_\mu (1-4s_W^2 - \gamma_5) \ell_e \right\} + \delta L_E, \qquad (2.10) $$

$$ \delta L_E \supset \frac{ieG_F}{4\sqrt{2}} C_a \bar{\nu}_a \gamma_\mu (1-\gamma_5) \nu_a \partial_\nu F_{\nu\mu}, \qquad (2.11) $$

as a starting point. Here $a, b$ are flavour indices, and a sum over repeated indices is implied. The operators in eq. (2.10) are generated at tree-level, and we have also adopted tree-level values for the couplings (see below). The operator in eq. (2.11), where $F_{\nu\mu}$ is the QED field strength and $e$ is the electromagnetic coupling, is generated at 1-loop level, through the diagrams shown in figure 1 (cf., e.g., refs. [8, 9]). The coefficient $C_a$ includes a
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Standard Model processes generating the effective $\bar{\nu}_a\nu_a\gamma$ vertex in eq. (2.11). A cross denotes a counterterm, IR indicates that the amplitude is evaluated in a kinematic domain in which the effective theory (eft) is applicable, and a thick solid line in the eft contribution stands for an electron or positron. The subtraction removes low-energy contributions, which are regenerated in a more complete form through figure 2 (including thermal corrections).</figcaption></figure>

computable divergence and anomalous dimension, but its absolute value cannot be determined perturbatively, as it gets a contribution from hadronic effects through the $Z\gamma$ bubble. We return to a discussion of its effects at the end of this section (cf. eq. (2.17)).

The representation in eq. (2.10) can be brought into a more helpful form by making use of a Fierz identity. This is particularly transparent, if we work in the Weyl representation of the Euclidean Dirac matrices,

$$ \gamma_0 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \gamma_m = \begin{pmatrix} 0 & -i\sigma_m \\ i\sigma_m & 0 \end{pmatrix}, \quad 1 - \gamma_5 = \begin{pmatrix} 2 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (2.12) $$

Inserting $\sum_m (\sigma_m)_{ij}(\sigma_m)_{kl} = 2\delta_{il}\delta_{jk} - \delta_{ij}\delta_{kl}$, this gives

$$ [\gamma_{\mu}(1-\gamma_5)]_{ij}[\gamma_{\mu}(1-\gamma_5)]_{kl} = 8 \begin{pmatrix} 0 & 0 \\ \delta_{ij}\delta_{kl} & -\delta_{il}\delta_{jk} \end{pmatrix}. \qquad (2.13) $$

This structure is antisymmetric in $j \leftrightarrow l$. This is compensated for by another minus sign from the anticommutation of fermionic fields, whereby eq. (2.10) can finally be turned into

$$ L_E \supset \frac{G_F}{4\sqrt{2}} \left\{ \bar{\nu}_a \gamma_\mu (1-\gamma_5) \nu_a \bar{\nu}_b \gamma_\mu (1-\gamma_5) \nu_b \right. \\ \left. \qquad + 2 \bar{\nu}_a \gamma_\mu (1-\gamma_5) \nu_a \bar{\ell}_e \gamma_\mu [2\delta_{a,e} - 1 + 4s_W^2 + (1-2\delta_{a,e})\gamma_5] \ell_e \right\} + \delta L_E. \qquad (2.14) $$

This is the form that we employ for our practical computations.

Let us now return to the operator $\delta L_E$, defined in eq. (2.11). It can be written in the same form as the operators in eq. (2.14) by making use of equations of motion (eom), leading to

$$ \delta L_E|_{\text{eom}} \supset \frac{G_F}{4\sqrt{2}} \left\{ e^2 C_a \bar{\nu}_a \gamma_\mu (1-\gamma_5) \nu_a \bar{\ell}_e \gamma_\mu \ell_e \right\}. \qquad (2.15) $$

Thereby eq. (2.15) can be viewed as an $O(e^2)$ correction to the operators in eq. (2.14). Indeed, ref. [9] determined the NLO coefficients of the 4-fermion operators after this redefinition, by fixing the low-energy hadronic contribution to the $Z\gamma$ bubble from experimentally measured quantities [10]. However, the goal of the present paper is to compute the percentual change from QED corrections. Therefore, it is more transparent to keep the coefficients of the 4-fermion operators at their electroweak values, and track the QED corrections, such as those originating from eq. (2.11), separately. To achieve this, we need to reconstruct the value of $C_a$ from the 4-fermion coefficients estimated in ref. [9].
In order to achieve this, let us re-express the neutrino-electron interaction from eqs. (2.14) and (2.15) in the same form as in ref. [9],

$$ L_E + \delta L_E \supset 2\sqrt{2}G_F \bar{\nu}_a \gamma_\mu a_L \nu_a \bar{\ell}_e \gamma_\mu \left\{ \right. \\ \qquad + \left[ \delta_{a,e} \left( s_W^2 + \frac{e^2 C_a}{8} \right) + (1 - \delta_{a,e}) \left( s_W^2 + \frac{e^2 C_a}{8} \right) \right] a_R \\ \qquad + \left[ \delta_{a,e} \left( s_W^2 + \frac{1}{2} + \frac{e^2 C_a}{8} \right) + (1 - \delta_{a,e}) \left( s_W^2 - \frac{1}{2} + \frac{e^2 C_a}{8} \right) \right] a_L \left. \right\} \ell_e, \quad (2.16) $$

where $a_R \equiv (1 + \gamma_5)/2$ and $a_L \equiv (1 - \gamma_5)/2$ are the right and left-handed projectors. The coefficients $g_{eR}, \dots, g_{(\mu, \tau)L}$ correspond to a notation employed in the literature (cf., e.g., ref. [11]), where often no distinction is made between $a = \mu$ and $a = \tau$.

After renormalization, the finite parts of the coefficients $g_{aR}$ and $g_{aL}$ become running parameters. We fix them through their values at the $\overline{\text{MS}}$ scheme renormalization scale $\bar{\mu} = m_e$. The logarithmic running, induced by $C_a$ (cf. eq. (2.17)), matches a corresponding thermal logarithm, and ultimately yields the $\ln(m_e/T)$ visible in eq. (3.4).

To be concrete, the coefficients $2\sqrt{2}G_F g_{aR,L}$ from eq. (2.16) agree with the coefficients $c_{R,L}^{\nu_a \ell'}$ listed in table 4 of ref. [9]. Relevant for us is the case that only the charged flavour $\ell' = \ell_e$ is light. In fact, the comparison is overconstrained, as we have parametrized six independent coefficients $c_{R,L}^{\nu_a \ell_e}$, $a \in \{e, \mu, \tau\}$, through the three parameters $C_a$, and thereby omitted some (non-QED) electroweak corrections. Nevertheless, all values from ref. [9] can be conservatively encompassed by setting $G_F|_{\bar{\mu}=m_e} \simeq 1.1664 \times 10^{-5}/\text{GeV}^2$, $s_W^2|_{\bar{\mu}=m_e} \simeq 0.2386$, and

$$ \left. \frac{e^2 C_a}{8} \right|_{\text{bare}} = (2\delta_{a,e} - 1 + 4s_W^2) \frac{e^2}{3} \frac{\mu^{-2\epsilon}}{(4\pi)^2} \left( \frac{1}{\epsilon} + \ln \frac{\bar{\mu}^2}{m_e^2} \right) + (-0.01 \dots 0.01). \quad (2.17) $$

The first term contains the correct divergence (i.e. counterterm) and renormalization scale dependence; the last parentheses represent the uncertainty from our simplified procedure. We remark that this uncertainty is of the order of the smallest physical coupling entering our computation, cf. eq. (2.19).

An important point to appreciate is that QED physics is more transparent in the original vector-axial rather than in the new left-right basis. After taking care of renormalization in the left-right basis, we can transfer back to the vector-axial basis, through

$$ g_{aR}a_R + g_{aL}a_L = \underbrace{\frac{g_{aR} + g_{aL}}{2}}_{\equiv g_{aV}} + \underbrace{\frac{g_{aR} - g_{aL}}{2}}_{\equiv g_{aA}} \gamma_5. \quad (2.18) $$

One reason for the simplification is that QED is a vectorlike theory, and therefore the QED running, induced by $C_a$, affects the vector coupling, but cancels in the axial one, cf. eq. (2.15). Another reason is that the vector coupling is very small for $a \neq e$,

$$ g_{(\mu, \tau)\text{V}} = s_W^2 - \frac{1}{4} + \frac{e^2 C_{(\mu, \tau)}}{8} \quad \underset{\bar{\mu} = m_e}{\overset{\overline{\text{MS}}}{\sim}} -0.01 + (-0.01 \dots 0.01). \quad (2.19) $$
A virtual electron-positron loop couples to a photon only through $g_{aV}$. Therefore the flavours $a \neq e$ are almost insensitive to scatterings mediated by photon exchange. As we will see, this induces a significant difference to the thermal interaction rates felt by $\nu_e$ and $\nu_{a \neq e}$.

## 2.4 Contractions for neutrino self-energy

The next task is to determine the medium-modified (i.e. full) neutrino propagator. After resumming 1-loop self-energy insertions, and showing explicitly the chiral projectors in accordance with the appearance of only left-handed neutrinos in eq. (2.14), the inverse propagator is expressed as

$$ \Delta_{K;\nu_a}^{(\text{full})} = a_R (i\not K + \Sigma_{K;\nu_a}) a_L. \qquad (2.20) $$

If we furthermore define the electron currents

$$ J_{\mu} \equiv \bar{l}_{e} \gamma_{\mu} l_{e}, \quad J_{\mu}^{5} \equiv \bar{l}_{e} \gamma_{\mu} \gamma_{5} l_{e}, \qquad (2.21) $$

and their fermion-line connected (“c”) correlation functions

$$ V_{P;\mu\nu} \equiv \int_X e^{iP \cdot X} \langle J_\mu(X) J_\nu(0) \rangle_c, \qquad (2.22) $$

$$ A_{P;\mu\nu} \equiv \int_X e^{iP \cdot X} \langle J_{\mu}^{5}(X) J_{\nu}^{5}(0) \rangle_c, \qquad (2.23) $$

$$ M_{P;\mu\nu} \equiv \int_X e^{iP \cdot X} \langle J_\mu(X) J_\nu^5(0) + J_\mu^5(X) J_\nu(0) \rangle_c, \qquad (2.24) $$

then the $O(G_F^2)$ contribution to the self-energy becomes

$$ \Sigma_{K;\nu_a} \supset \frac{G_F^2}{8} \gamma_\mu (1-\gamma_5) \Delta_{K+P;\nu_a}^{-1} \mathcal{S}_{P;\mu\nu}^a \gamma_\nu (1-\gamma_5), \qquad (2.25) $$

$$ \begin{aligned} \mathcal{S}_{P;\mu\nu}^a = & \sum_b \text{Tr} \left[ \Delta_{P+Q;\nu_b}^{-1} \gamma_\mu (1-\gamma_5) \Delta_{Q;\nu_b}^{-1} \gamma_\nu (1-\gamma_5) \right] \\ & - \gamma_\nu (1-\gamma_5) \Delta_{P+Q;\nu_a}^{-1} \gamma_\mu (1-\gamma_5) \Delta_{Q;\nu_a}^{-1} \\ & - \left[ (2\delta_{a,e} - 1 + 4s_W^2)^2 V_{P;\mu\nu} + (1-2\delta_{a,e})(2\delta_{a,e} - 1 + 4s_W^2) M_{P;\mu\nu} + A_{P;\mu\nu} \right] \\ & + e^2 (2\delta_{a,e} - 1 + 4s_W^2)^2 V_{P;\mu\rho}^{\text{LO}} \Delta_{P;\rho\sigma}^{-1} V_{P;\sigma\nu}^{\text{LO}} \\ & + e^2 (2\delta_{a,e} - 1 + 4s_W^2) C_a V_{P;\mu\rho}^{\text{LO}} \Delta_{P;\rho\sigma}^{-1} (P_\sigma P_\nu - P^2 \delta_{\sigma\nu}). \end{aligned} \qquad (2.26) $$

Here $\Delta_{K;\nu_a}^{-1} \equiv (i\not K)^{-1}$ is the tree-level neutrino propagator, $\Delta_{P;\rho\sigma}^{-1}$ is the tree-level photon propagator, and sum-integrals over repeated four-momenta $(P, Q)$ are implied. The last two lines of eq. (2.26) include fermion-line “disconnected” contributions, originating from the last term in figure 2(b), and “LO” stands for a leading-order evaluation of eq. (2.22).

Now, particles can to a good approximation be treated as massless if their mass is less than $\sim \pi T$. If we consider temperatures $T \gtrsim 1$ MeV, the electron mass does satisfy this requirement. Therefore, the electron mass can be set to zero, simplifying the analysis.

Let us rewrite eq. (2.26) in the massless limit. The neutrino propagators anticommute with $\gamma_5$, and in the massless limit the electron propagators do the same. The mixed correlator
from eq. (2.24) vanishes. Taking also into account that the open $(1 - \gamma_5)$'s of the second term of eq. (2.26) can be transported into the projectors in eq. (2.25), we thereby obtain

$$ \lim_{m_e \to 0} \Sigma_{K;\nu_a} \supset \frac{G_F^2}{4} \gamma_\mu \Delta_{K+P;\nu_a}^{-1} \tilde{\mathcal{S}}_{P;\mu\nu}^a \gamma_\nu (1 - \gamma_5), \qquad (2.27) $$

$$ \begin{aligned} \tilde{\mathcal{S}}_{P;\mu\nu}^a ={}& (2 \sum_b) \operatorname{Tr} [\Delta_{P+Q;\nu_b}^{-1} \gamma_\mu \Delta_{Q;\nu_b}^{-1} \gamma_\nu] - 4 \gamma_\nu \Delta_{P+Q;\nu_a}^{-1} \gamma_\mu \Delta_{Q;\nu_a}^{-1} \\ & - [ (2\delta_{a,e} - 1 + 4s_W^2)^2 + 1 ] V_{P;\mu\nu} \\ & + e^2 (2\delta_{a,e} - 1 + 4s_W^2)^2 V_{P;\mu\rho}^{\mathrm{LO}} \Delta_{P;\rho\sigma}^{-1} V_{P;\sigma\nu}^{\mathrm{LO}} \\ & + e^2 (2\delta_{a,e} - 1 + 4s_W^2) C_a V_{P;\mu\rho}^{\mathrm{LO}} \Delta_{P;\rho\sigma}^{-1} (P_\sigma P_\nu - P^2 \delta_{\sigma\nu}). \end{aligned} \qquad (2.28) $$

## 2.5 Carrying out a Matsubara sum and angular integrals

Our next goal is to carry out the outermost Matsubara sum-integral in eq. (2.27), over $P$. In order to simplify the notation, let us denote the second part of the integrand by $\Pi_P$. We write it in a spectral representation,

$$ \Pi_P = \int_{-\infty}^{\infty} \frac{dp^0}{\pi} \frac{\operatorname{Im} \Pi_P}{p^0 - ip_n}, \quad \operatorname{Im} \Pi_P \equiv \operatorname{Im} \Pi_P|_{p_n \to -i[p^0 + i0^+]}, \qquad (2.29) $$

where $p_n$ is bosonic. Then we are faced with the structure

$$ \begin{aligned} \Sigma_K &= \int_P \frac{\alpha(k_n + p_n) + \beta(\mathbf{k} + \mathbf{p})}{(k_n + p_n)^2 + (\mathbf{k} + \mathbf{p})^2} \Pi_P \\ &= \int_{\mathbf{p}} \int_{-\infty}^{\infty} \frac{dp^0}{\pi} T \sum_{p_n, \{r_n\}} \delta_{0,k_n+p_n-r_n} \frac{\alpha r_n + \beta(\mathbf{k} + \mathbf{p})}{r_n^2 + \epsilon_{kp}^2} \frac{\operatorname{Im} \Pi_P}{p^0 - ip_n}, \end{aligned} \qquad (2.30) $$

where $k_n$ and $r_n$ are fermionic, and $\epsilon_{kp}^2 \equiv (\mathbf{k} + \mathbf{p})^2$. We may now write $\delta_{0,k_n+p_n-r_n} = T \int_0^{1/T} d\tau e^{i(k_n+p_n-r_n)\tau}$, and carry out the sums over $p_n$ and $r_n$. Furthermore the integral over $\tau$ is doable, yielding

$$ \Sigma_K = \int_{\mathbf{p}} \int_{-\infty}^{\infty} \frac{dp^0}{\pi} \frac{1}{2\epsilon_{kp}} \left\{ \begin{aligned} & [i\alpha\epsilon_{kp} - \beta(\mathbf{k} + \mathbf{p})] \frac{n_F(\epsilon_{kp}) + n_B(p^0)}{ik_n + p^0 - \epsilon_{kp}} \\ & + [i\alpha\epsilon_{kp} + \beta(\mathbf{k} + \mathbf{p})] \frac{1 - n_F(\epsilon_{kp}) + n_B(p^0)}{ik_n + p^0 + \epsilon_{kp}} \end{aligned} \right\} \operatorname{Im} \Pi_P, \qquad (2.31) $$

where $n_B$ and $n_F$ denote the Bose and Fermi distributions, respectively.

As a next step, we continue to Minkowskian frequency, $ik_n \to k^0 + i0^+$, and take the imaginary part. We also substitute $p^0 \to -p^0$, $\mathbf{p} \to -\mathbf{p}$, and note that bosonic spectral functions are odd under this substitution. Finally, $\epsilon_{kp}$ is redefined as $\epsilon_{kp} \equiv |\mathbf{k} - \mathbf{p}|$, and we recall that $n_B(-x) = -1 - n_B(x)$, $n_F(-x) = 1 - n_F(x)$. This yields

$$ \begin{aligned} \operatorname{Im} \Sigma_K ={}& \int_{\mathbf{p}} \int_{-\infty}^{\infty} \frac{dp^0}{2\epsilon_{kp}} [i\alpha(p^0 - k^0) + \beta(\mathbf{k} - \mathbf{p})] [1 - n_F(k^0 - p^0) + n_B(p^0)] \\ & \times [\delta(k^0 - p^0 - \epsilon_{kp}) - \delta(k^0 - p^0 + \epsilon_{kp})] \operatorname{Im} \Pi_P. \end{aligned} \qquad (2.32) $$
Subsequently, we may integrate over the angles of $\mathbf{p}$, removing the Dirac-$\delta$'s. At this point we also go to the light cone, $k^0 = k$, relevant for on-shell neutrinos. The first Dirac-$\delta$ can be seen to be realized for $|p^0| < p$, i.e. it corresponds to a $t$-channel process. The second Dirac-$\delta$ is realized for $p^0 > p$, and corresponds to an $s$-channel process. For both channels, the Dirac-$\delta$'s imply

$$ \mathbf{k} \cdot \mathbf{p} = \frac{k^2 + p^2 - \epsilon_{kp}^2}{2} \stackrel{k^0=k}{=} k p^0 - \frac{\mathcal{P}^2}{2}. \qquad (2.33) $$

Introducing

$$ p_{\pm} \equiv \frac{p^0 \pm p}{2}, \qquad (2.34) $$

eq. (2.32) can be turned into

$$ \mathrm{Im} \, \Sigma_k = \frac{1}{4\pi^2 k} \left( \underbrace{\int_{-\infty}^{0} dp_{-} \int_{0}^{k} dp_{+}}_{\textit{t-channel}} - \underbrace{\int_{0}^{k} dp_{-} \int_{k}^{\infty} dp_{+}}_{\textit{s-channel}} \right) \times p \left[ i\alpha(p^0 - k) + \beta(\mathbf{k} - \mathbf{p}) \right] \left[ 1 - n_F(k - p^0) + n_B(p^0) \right] \mathrm{Im} \, \Pi_{\mathcal{P}}. \qquad (2.35) $$

In this section we have dealt with the outer sum-integral in eq. (2.27). However, the same procedure can also be used for the inner one, except that then the four-momentum $\mathcal{P}$ is not light-like. In some cases, it also turns out to be helpful to redefine variables, whereby a clean separation into an outer and inner sum-integration no longer applies (cf. appendix A.1).

## 2.6 Projecting out the interaction rate

The object in eq. (2.35) is a matrix in Dirac space, cf. eq. (2.27). Assuming that all the sum-integrals have been carried out, going over to Minkowskian signature (whereby $i\mathcal{K} \to \mathcal{K}$), and redundantly displaying chiral projectors in accordance with eq. (2.20), the result has the form [12]

$$ \Sigma_{k;\nu_a} = a_R (a\mathcal{K} + b\not\!\!\nu) a_L, \qquad (2.36) $$

where $u \equiv (1, 0)$ denotes the medium four-velocity.

The interaction rate influences the time evolution of on-shell active neutrinos (cf. eq. (2.2)). In terms of eq. (2.36), it can be found in the imaginary part of $b$,

$$ b = b_r + \frac{i}{2} \Gamma_{k;\nu_a}, \quad b_r \in \mathbb{R}. \qquad (2.37) $$

The real part $b_r$, which gives a medium modification to the effective Hamiltonian $\Delta\mathcal{E}_k$ appearing in eq. (2.2), is of $\mathcal{O}(G_F)$ [7]. The imaginary part, which is of $\mathcal{O}(G_F^2)$ and interests us here, can be projected out as

$$ \Gamma_{k;\nu_a} = \frac{1}{k} \mathrm{Im} \, \mathrm{Tr} [\mathcal{K} \Sigma_{k;\nu_a}]. \qquad (2.38) $$

As a next step, we put together the ingredients from eqs. (2.27), (2.28), (2.35) and (2.38). First we consider the terms from eq. (2.28) that have a factorized Dirac trace, as is indeed
the case in all but one term. We then find

$$ \Gamma_{k;\nu_a} \supset \mathrm{Im} \left\{ \frac{G_F^2}{4k} \frac{\mathrm{Tr} [(i\not K)\gamma_\mu(-i)(\not K + \not P)\gamma_\nu(1-\gamma_5)]}{(K+P)^2} \tilde{S}_{P;\mu\nu}^a \right\} \\ \rightarrow \mathrm{Im} \left\{ \frac{G_F^2}{k} \frac{2K_\mu K_\nu - \delta_{\mu\nu} K \cdot (K+P)}{(K+P)^2} \tilde{S}_{P;\mu\nu}^a \right\}, \qquad (2.39) $$

where we noted that the part involving $\gamma_5$ drops out, if $\tilde{S}_{P;\mu\nu}^a$ is symmetric in $\mu \leftrightarrow \nu$, and that the parts proportional to $P_\mu$ and $P_\nu$ drop out, if $\tilde{S}_{P;\mu\nu}^a$ is transverse.

Subsequently, we make use of eq. (2.35), with $\tilde{S}_{P;\mu\nu}^a$ playing the role of $\Pi_P$. Inserting eq. (2.33), this leads to

$$ \Gamma_{k;\nu_a} \supset \frac{G_F^2}{8\pi^2 k^2} \left( \underbrace{\int_{-\infty}^0 \mathrm{d}p_- \int_0^k \mathrm{d}p_+}_{t\text{-channel}} - \underbrace{\int_0^k \mathrm{d}p_- \int_k^\infty \mathrm{d}p_+}_{s\text{-channel}} \right) \times p [4K_\mu K_\nu - \delta_{\mu\nu} \mathcal{P}^2] [1 - n_F(k-p^0) + n_B(p^0)] \mathrm{Im} \tilde{S}_{P;\mu\nu}^a. \qquad (2.40) $$

To complete the transition to Minkowskian spacetime, we recall the rules from eq. (2.7), resulting in

$$ \Gamma_{k;\nu_a} \supset \frac{G_F^2}{8\pi^2 k^2} \left( \underbrace{\int_{-\infty}^0 \mathrm{d}p_- \int_0^k \mathrm{d}p_+}_{t\text{-channel}} - \underbrace{\int_0^k \mathrm{d}p_- \int_k^\infty \mathrm{d}p_+}_{s\text{-channel}} \right) \times p [-4K_\mu K_\nu - \eta_{\mu\nu} \mathcal{P}^2] [1 - n_F(k-p^0) + n_B(p^0)] \mathrm{Im} \tilde{S}_{P}^{a;\mu\nu}. \qquad (2.41) $$

As a final ingredient, we take into account that the spectral functions originating from the factorized terms in eq. (2.28) are transverse. Anticipating section 3, let us focus on

$$ \mathrm{Im} \tilde{S}_{P}^{a;\mu\nu} \supset -[(2\delta_{a,e} - 1 + 4s_W^2)^2 + 1] \mathrm{Im} V_P^{\mu\nu}. \qquad (2.42) $$

The vector channel spectral function can be decomposed as

$$ \mathrm{Im} V_P^{\mu\nu} = P_T^{\mu\nu} \rho_T + P_L^{\mu\nu} \rho_L, \qquad (2.43) $$

where the projectors have been defined as

$$ P_T^{\mu\nu} \equiv -\eta_i^\mu \eta_j^\nu \left( \delta_{ij} - \frac{p_i p_j}{p^2} \right), \quad P_L^{\mu\nu} \equiv \eta^{\mu\nu} - \frac{\mathcal{P}^\mu \mathcal{P}^\nu}{\mathcal{P}^2} - P_T^{\mu\nu}. \qquad (2.44) $$

Inserting this representation, and making use of eq. (2.33) in order to resolve the scalar products $\mathbf{k} \cdot \mathbf{p}$ originating from $\mathcal{K}_\mu \mathcal{K}_\nu \mathcal{P}_T^{\mu\nu}$, the contraction in eq. (2.41) takes the form

$$ \begin{aligned} &[-4\mathcal{K}_\mu \mathcal{K}_\nu - \eta_{\mu\nu} \mathcal{P}^2] \mathrm{Im} V_P^{\mu\nu} \\ &= -\frac{\mathcal{P}^2}{2} \left\{ (2\rho_T + \rho_L) \left[ 1 + \left( \frac{2k-p^0}{p} \right)^2 \right] + \rho_L \left[ 1 - 3 \left( \frac{2k-p^0}{p} \right)^2 \right] \right\}. \end{aligned} \qquad (2.45) $$

This representation will be employed in section 3, and can also be crosschecked numerically, as demonstrated in section A.2.
## 2.7 HTL resummation of soft t-channel photons

The computation described above employs tree-level propagators. While the result is UV finite, one of the integrals turns out to be logarithmically IR divergent. In the language of eq. (2.41), the IR divergence emerges from the t-channel contribution, when the momentum flowing through the photon propagator, $\Delta_{P;\rho\sigma}^{-1}$ in eq. (2.28), is “soft”, $|p^0| \le p \ll \pi T$. The technical reason is that the vector correlators turn into “Hard Thermal Loops” (HTLs) in this domain, $V_P^{LO} \sim T^2$. Then they do not regulate that IR regime any more, as would be the case in vacuum, where they behave as $\sim P^2$. In addition, soft photons are Bose-enhanced, with $n_B(p^0) \approx T/p^0$.

HTLs appear not only in the vertices $V_P^{LO}$, but they also influence the self-energies of the photons. When all effects scaling as $\sim T^2$ are included, we talk about HTL resummation [13, 14]. HTL resummation introduces screening effects (Debye screening, Landau damping) to the photon propagator, and these regulate the regime of soft momentum exchange.

In our case, when only one contribution is affected, HTL resummation is fairly simple to implement. The relevant structures originate from the IR limits of the real and imaginary parts of $V_P^{LO}$. The imaginary parts are obtained from eqs. (A.15) and (A.17),

$$ \rho_T^{\text{IR}} = [\rho_T^{\text{LO}}]_{p, |p^0| \ll \pi T} = + \frac{\pi p^0 \mathcal{P}^2 T^2}{12p^3} \theta(p - |p^0|), \quad (2.46) $$

$$ \rho_L^{\text{IR}} = [\rho_L^{\text{LO}}]_{p, |p^0| \ll \pi T} = - \frac{\pi p^0 \mathcal{P}^2 T^2}{6p^3} \theta(p - |p^0|). \quad (2.47) $$

The real parts can be derived from eqs. (A.22) and (A.24),

$$ \chi_T^{\text{IR}} = [\chi_T^{\text{LO}}]_{p, |p^0| \ll \pi T} = + \frac{T^2}{6p^2} \left[ (p^0)^2 - \frac{p^0 \mathcal{P}^2}{2p} \ln \left| \frac{p+p^0}{p-p^0} \right| \right], \quad (2.48) $$

$$ \chi_L^{\text{IR}} = [\chi_L^{\text{LO}}]_{p, |p^0| \ll \pi T} = - \frac{\mathcal{P}^2 T^2}{3p^2} \left[ 1 - \frac{p^0}{2p} \ln \left| \frac{p+p^0}{p-p^0} \right| \right]. \quad (2.49) $$

With the help of eqs. (2.46)–(2.49), we can construct the resummed photon propagator, denoted by $\Delta_{P;\mu\nu}^{-1*}$. After analytic continuation, its transverse part can be decomposed in analogy with eq. (2.43) (the longitudinal part depends on the gauge parameter and drops out when inserted into physical observables). For convenience we insert an overall minus sign, and denote the real and imaginary parts of the resummed propagator by $\mathcal{R}_{T,L}^*$ and $\mathcal{I}_{T,L}^*$,

$$ -\Delta_{\mathcal{P}}^{-1*;\mu\nu} = \mathbb{P}_{\mathcal{T}}^{\mu\nu} (\mathcal{R}_{\mathcal{T}}^* + i\mathcal{I}_{\mathcal{T}}^*) + \mathbb{P}_{\mathcal{L}}^{\mu\nu} (\mathcal{R}_{\mathcal{L}}^* + i\mathcal{I}_{\mathcal{L}}^*). \quad (2.50) $$

With the help of eqs. (2.46)–(2.49), we find

$$ \mathcal{R}_{T,L}^* = \text{Re} \left\{ \frac{1}{\mathcal{P}^2 - e^2(\chi_{T,L}^{\text{IR}} + i\rho_{T,L}^{\text{IR}})} \right\} = \frac{\mathcal{P}^2 - e^2\chi_{T,L}^{\text{IR}}}{(\mathcal{P}^2 - e^2\chi_{T,L}^{\text{IR}})^2 + (e^2\rho_{T,L}^{\text{IR}})^2}, \quad (2.51) $$

and correspondingly for the imaginary parts.

In order to implement the resummation in an unambiguous way, we only employ it in terms where it is necessary at $O(e^2)$, but not in terms where its effect is of higher order in $e$. Notably, in the t-channel, we only keep the contribution from $\mathcal{R}_{T,L}^*$, but not from $\mathcal{I}_{T,L}^*$, since
the latter leads to a contribution of $\mathcal{O}(e^4)$.<fnref n="3" /> In the s-channel, we set $\mathcal{R}_{T,L}^* \to \mathcal{R}_{T,L} \equiv 1/\mathcal{P}^2$, since the integral is IR-finite. Furthermore, we omit the effect from $\mathcal{I}_{T,L}^*$ in the s-channel as well. It would lead to the processes illustrated in set (h) of figure 2. These correspond to the famous plasmon contributions (cf., e.g., ref. [5]). However, they are formally of $\mathcal{O}(e^3)$ rather than $\mathcal{O}(e^2)$: apart from the vertices, the result of a $1 \leftrightarrow 2$ reaction is phase-space suppressed by the plasmon mass $\sim eT$.

After the inclusion of HTL resummation in the t-channel, the contribution from soft momentum exchange is finite. However, it is logarithmically enhanced, by $\ln(1/e^2)$. The coefficient of this logarithm can be worked out analytically, as shown in appendix B. The final result, formally the largest NLO QED correction that we find, is displayed in eq. (3.5).

### 3 NLO computation and result

We now put together the full NLO expression for the thermal interaction rate of a neutrino of flavour a and momentum k. The result is expressed in terms of partial interaction rates, in the spirit of eq. (2.3). We make use of the Weinberg angle parametrization from eq. (2.9), the simplified integrals in eqs. (A.12) and (A.13) for scatterings off neutrinos, the full representation of the $e^+e^-$ contribution from eqs. (2.41), (2.42) and (2.45), the two independent electron-positron spectral functions defined according to eq. (2.43), as well as HTL resummation from eq. (2.51). We assume all ensembles to carry the same temperature, and omit chemical potentials. Then scatterings off, pair annihilations into, or pair creations from, the same neutrino flavour, different neutrino flavours, and electrons and positrons, yield the contributions

$$ \Gamma_{k;\nu_a}^{(\nu_a)} = \{2+2\} \frac{G_F^2}{16\pi^3 k^2} [2\Xi^t + \Xi^s], \quad (3.1) $$

$$ \Gamma_{k;\nu_a}^{(\nu_b \neq \nu_a)} = \left\{ \left( 2 \sum_{b \neq a} \right) \right\} \frac{G_F^2}{16\pi^3 k^2} [2\Xi^t + \Xi^s], \quad (3.2) $$

$$ \begin{aligned} \Gamma_{k;\nu_a}^{(e^+e^-)} = & \frac{G_F^2}{8\pi^2 k^2} \left( \underbrace{\int_{-\infty}^0 dp_- \int_0^k dp_+}_{\text{t-channel}} - \underbrace{\int_0^k dp_- \int_k^\infty dp_+}_{\text{s-channel}} \right) p \mathcal{P}^2 [1 - n_F(k-p^0) + n_B(p^0)] \\ & \times \left\{ \left[ (2\delta_{a,e} - 1 + 4s_W^2)^2 + 1 \right] \left[ \rho_T^{\text{NLO}} + \rho_L^{\text{NLO}} + (\rho_T^{\text{NLO}} - \rho_L^{\text{NLO}}) \left(\frac{2k-p^0}{p}\right)^2 \right] \right. \\ & \left. + 2e^2(2\delta_{a,e} - 1 + 4s_W^2)^2 \left[ \rho_T^{\text{LO}} \chi_T^{\text{LO}} \mathcal{R}_T^* + \rho_L^{\text{LO}} \chi_L^{\text{LO}} \mathcal{R}_L^* + (\rho_T^{\text{LO}} \chi_T^{\text{LO}} \mathcal{R}_T^* - \rho_L^{\text{LO}} \chi_L^{\text{LO}} \mathcal{R}_L^*) \left(\frac{2k-p^0}{p}\right)^2 \right] \right. \\ & \left. + e^2 C_a (2\delta_{a,e} - 1 + 4s_W^2) \left[ \rho_T^{\text{LO}} + \rho_L^{\text{LO}} + (\rho_T^{\text{LO}} - \rho_L^{\text{LO}}) \left(\frac{2k-p^0}{p}\right)^2 \right] \right\}, \end{aligned} \quad (3.3) $$

respectively, where $p_\pm$ are defined according to eq. (2.34); $\mathcal{P}^2 \equiv (p^0)^2 - p^2 = 4p_+p_-$; and $n_F$ and $n_B$ denote the Fermi and Bose distribution functions. The functions $\rho_{T,L}^{\text{LO,NLO}}$ originate

<footnote n="3">³In terms of figure 2, this corresponds to the square of the t-channel amplitude in set (g), rather than its interference with the tree-level process.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) leading-order reactions between neutrinos (arrowed lines) and electrons and positrons (thick lines, with arrows omitted, understanding that both directions should be included); (b) a dispersive representation of the NLO corrections, with the vertical line denoting a cut, a wiggly line a photon, an asterisk HTL-resummation, and the filled box the loop-generated operator from eq. (2.11); (c) virtual corrections to s-channel reactions; (d) real corrections to s-channel reactions; (e) virtual corrections to t-channel reactions; (f) real corrections to t-channel reactions; (g) s and t-channel processes involving virtual photon exchange, originating from the last term in set (b); (h) s-channel processes involving a plasmon resonance, originating from the last term in set (b).</figcaption></figure>

from the imaginary part of the vector current correlator in eq. (2.22), and are discussed below; the functions $\chi_{T,L}^{LO}$ originate from its real part, and are given in appendix A.3.

Let us elaborate on the three contributions in eq. (3.3). The first line contains the NLO spectral functions $\rho_{T,L}^{NLO}$. These represent the diagrams in figure 2(c)–(f): all of these processes, and the interference terms between them, are contained in the evaluations of $\rho_{T,L}^{NLO}$ that were carried out in refs. [15, 16]. We just need to adjust the QCD coupling $g^2$ and the Casimir factor $C_F$ as $g^2C_F \to e^2 \equiv 4\pi\alpha_{em}$. We remark in passing that both spectral functions are independent of the renormalization scale at NLO.
The second line of eq. (3.3) originates from the processes in figure 2(g), once these interfere with LO contributions. In the notation of the last term of figure 2(b), the open box is replaced by a closed loop here; it is the closed loop which produces the functions $\chi_{T,L}^{LO}$. We note that the t-channel part of these processes is IR divergent in naive perturbation theory, and therefore formally the single largest NLO contribution. After HTL resummation through $\mathcal{R}_{T,L}^*$, as explained in section 2.7, it develops a logarithmic dependence on $e^2$, cf. eq. (3.5).

The third line of eq. (3.3) also originates from the processes in figure 2(g), but now the open box is replaced by the effective vertex from eq. (2.11), parametrized by the coefficient $C_a$. This cancels the UV divergences of $\chi_{T,L}^{LO}$, but also introduces a finite contribution, cf. eq. (2.17). Similarly to photon vacuum polarization, the coefficient $C_a$ is non-perturbative, as it is affected by a low-energy hadronic contribution from the $Z\gamma$ bubble, cf. section 2.3.

Apart from the cancellation of the divergence, the concrete role of $C_a$ is to replace the renormalization scale $\bar{\mu}$ in eq. (A.20) through the physical scale $m_e$. If we represent momenta as $k/T$, our results thereby obtain a logarithmic sensitivity on $m_e/T$.

For a numerical illustration, we rewrite eq. (3.3) as

$$ \frac{\Gamma_{k;\nu_a}^{(e^+e^-)}}{G_F^2 T^4 k/(8\pi^2)} = \left[ (2\delta_{a,e} - 1 + 4s_W^2)^2 + 1 \right] \left\{ \underbrace{A}_{\text{from } \rho_{T,L}^{LO}} + \underbrace{e^2 B}_{\text{from } \rho_{T,L}^{NLO}} \right\} \\ + e^2 (2\delta_{a,e} - 1 + 4s_W^2)^2 \left\{ \frac{1}{6\pi^2} \underbrace{\left[C - \left(2\ln\frac{m_e}{T} + \frac{5}{3}\right)A\right]}_{\text{from } C_a \text{ and } \chi_{T,L|_{T}}^{LO}|_{\text{vac}}} + \underbrace{D}_{\text{from } \chi_{T,L|_{T}}^{LO}\mathcal{R}_{T,L}^*} \right\} \\ + e^2 (2\delta_{a,e} - 1 + 4s_W^2) \underbrace{\left\{8(-0.01 \dots 0.01)A\right\}}_{\text{from uncertainty of } C_a} \tag{3.4} $$

The couplings are evaluated at the scale $\bar{\mu} = m_e$ as specified above eq. (2.17). The representation is reliable as long as the logarithm $\ln(m_e/T)$ is of order unity. The coefficients can be decomposed into t and s-channel contributions, $A = A^{(t)} + A^{(s)}$. Furthermore, $D^{(t)}$ is logarithmically sensitive to the coupling, and we express it as (cf. appendix B)

$$ D^{(t)} = D_{1}^{(t)} \ln\left(\frac{1}{e^2}\right) + D_{2}^{(t)}, \quad D_{1}^{(t)} = -\frac{\pi T}{9k}. \tag{3.5} $$

The resulting values are illustrated numerically in appendix C. Including the flavour factors, and the uncertainty from our estimate for $C_a$, the final results are shown in figure 3.

## 4 Conclusions

The general context of our investigation is the observable known as $N_{\text{eff}}$, parametrizing the energy density of the universe at the time of primordial nucleosynthesis or photon decoupling through an effective number of massless neutrino species. In order to scrutinize the current theoretical uncertainty of the Standard Model prediction of $N_{\text{eff}}$, ref. [4] recently carried out an approximate NLO computation of the energy transfer rate defined in eq. (1.1). They found a $\sim -5\%$ QED correction to the energy transfer rate in the temperature range $T \in (1...3)$ MeV.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Numerical results for the percentual NLO QED contribution to the $e^{+}e^{-}$ part of the neutrino interaction rate, cf. eq. (3.3), in massless QED, with $\alpha_{em} = 1/137$. The uncertainty, indicated with a band, reflects the low-energy hadronic contribution to the coefficient $C_a$, as specified in eq. (2.17). Top row: dependence on $k/T$ for fixed $T$/MeV, for $a = e$ (left) and $a \neq e$ (right). Bottom row: dependence on $T$/MeV for fixed $k/T$, for $a = e$ (left) and $a \neq e$ (right).</figcaption></figure>

When this correction was inserted into momentum-averaged kinetic equations describing neutrino decoupling [11], the fourth significant digit of the established Standard Model prediction of $N_{\text{eff}}$ [1–3] was asserted to decrease by one unit.

The purpose of the present paper has been to compute NLO QED corrections to the thermal interaction rate felt by active neutrinos in the same temperature range. The thermal interaction rate is defined by eq. (2.38) and influences neutrino physics via a momentum-dependent kinetic equation, of which eq. (2.2) is a simplified prototype. Therefore, it does not coincide with the energy transfer rate in eq. (1.1). However, it gets contributions from
the same diagrams (on our side, this is the subset called s-channel or inelastic processes, cf. figure 2), and should reflect similar physics.

In order to streamline the task, we have made a few simplifications concerning the overall statistical ensemble. In particular, we assumed the presence of a universal temperature, a vanishing electron mass, and vanishing lepton chemical potentials. However, we undertook no simplifications on the technical side, including in particular the full Bose and Fermi distributions, all processes depicted in figure 2, logarithmic running induced by the effective vertex in eq. (2.11), and HTL resummation of soft t-channel photon exchange, which is otherwise IR divergent due to Bose enhancement. Moreover we have computed a differential rate, displaying its momentum dependence.

Our numerical result is shown in figure 3. We find a $\sim -(0...2)%$ contribution, dominated by the t-channel diagrams in figure 2(e,f,g), represented by the coefficients $B^{(t)}$, $C^{(t)}$ and $D^{(t)}$, illustrated numerically in appendix C. If we restrict to the s-channel processes, contributing to the observable in eq. (1.1), the magnitude of the correction remains similar, but its sign is positive. The overall magnitude of the s-channel correction is noticeably smaller than that found in ref. [4], and its sign is opposite. It is well-known from other contexts that incomplete evaluations of NLO corrections may overestimate a result, because cancellations can go amiss (cf., e.g., ref. [17]), however we are not sure if this could be an explanation here. It might also be that the different observable of eq. (1.1), and/or effects related to the electron mass, could be responsible for the larger correction seen in ref. [4].

Another difference between our results and those in ref. [4] is that we observe a visible dependence of the relative NLO correction on the neutrino flavour $a$, whereas ref. [4] stated the difference between the flavours to be less than 1%. We remark that flavour dependence originates from the second and third lines of eqs. (3.3) and (3.4), corresponding to the processes in figure 2(g), since the dependence on $a$ differs from that of the leading-order term, on the first line. To be explicit, the t-channel contribution contains the diagram

<figure><img src="image_6.png" /></figure>

where “V” indicates that the QED vertices are purely vectorlike. As discussed around eq. (2.19), the flavours $a \neq e$ are almost blind to the vectorlike coupling, and therefore to this process. The same holds for the s-channel reflections of this process, depicted as the first two diagrams in figure 2(g). At the same time, the flavour $a = e$ gets a visible contribution from the vectorlike coupling, in both the t and s channels, for instance through the logarithm $\sim \ln(T/m_e)$ in eq. (3.4). Therefore, different flavours do feel different QED corrections.

Finally, we list possible directions of future research. On the mundane side, the effects of a finite electron mass could be incorporated at $O(e^2)$, by making use of the formalism developed in ref. [18], however this could turn out to be quite tedious numerically.

A challenge of a more conceptual nature would be to generalize the framework so that non-equilibrium ensembles can be inserted for the different neutrino flavours. The QED plasma, inducing the NLO corrections that we have computed here, would however maintain the same form as in the present investigation.
Last but not least, a significant role for the physics under consideration is played by the coefficient $C_a$, originating from the operator in eq. (2.11). It would be important to have a good estimate of the magnitude of its finite part, but this is non-trivial, due to low-energy hadronic loops that affect $Z\gamma$ mixing. Our study made use of the range in eq. (2.17), extracted from refs. [9, 10], however the error, visible as bands in figure 3, could conceivably be reduced through further theoretical and phenomenological work.

## Acknowledgments

G.J. thanks Wouter Dekens for a pointer to ref. [9]. M.L. thanks Miguel Escudero for discussions about the motivation for this work [4], for drawing our attention to ref. [5], and for comments on the manuscript, and Dietrich Bödeker for many helpful discussions. G.J. was funded by the U.S. Department of Energy (DOE), under grant No. DE-FG02-00ER41132, and by the Agence Nationale de la Recherche (France), under grant ANR-22-CE31-0018 (AUTOTHERM).

**Note added.** The very recent ref. [19] considered a subset of the processes in our study, but included the electron mass. Its conclusions are similar to ours, confirming in particular that $m_e$ can generally be set to zero at decoupling temperatures.

## A Evaluation at leading order

In this appendix we show how the expression for the thermal neutrino interaction rate can be simplified at leading order, by making use of symmetrizations of the sum-integration variables (cf. section A.1). We demonstrate how this simplification permits for a numerical crosscheck of the full integral representation in eq. (2.45) (cf. section A.2). Finally, we specify the real part of the correlator from eq. (2.22), as this is needed in eq. (3.3) (cf. section A.3).

### A.1 Additional simplifications of the leading-order (LO) expression

Let us return to eqs. (2.27), (2.28) and (2.38), and consider the most basic scenario, namely assuming the absence of any chemical potentials,<fnref n="4" /> and setting the temperature of each neutrino flavour equal. In this limit, the integrand develops additional symmetries, and consequently substitutions of sum-integration variables allow us to put the result in a simple form.

As a starting point, eqs. (2.27), (2.28) and (2.38) yield (omitting contributions of $O(e^2)$)

$$ \frac{\Gamma_{k;\nu_a}^{\text{LO}}}{G_F^2/(4k)} = \left\{ (2\sum_b) + \left[ (2\delta_{a,e} - 1 + 4s_W^2)^2 + 1 \right] \right\} \times \text{Im} \left\{ \frac{\text{Tr}[(i\not{K})\gamma_\mu(-i)(\not{K} + \not{P})\gamma_\nu(1-\gamma_5)]\text{Tr}[(\not{-i})(\not{P} + \not{Q})\gamma_\mu(-i)(\not{Q})\gamma_\nu(1-\gamma_5)]}{(K+P)^2(P+Q)^2Q^2} \right\} + (-4) \text{Im} \left\{ \frac{\text{Tr}[(i\not{K})\gamma_\mu(-i)(\not{K} + \not{P})\gamma_\nu(-i)(\not{P} + \not{Q})\gamma_\mu(-i)(\not{Q})\gamma_\nu(1-\gamma_5)]}{(K+P)^2(P+Q)^2Q^2} \right\}. \quad (A.1) $$

<footnote n="4">⁴We have set the chemical potentials to vanish all along, but stress the point here, because it would otherwise play an important role, removing the symmetry between particles and antiparticles.</footnote>
The first simplification originates from the substitution $Q \to -P - Q$. This leaves the denominator invariant, but symmetrizes the appearance of Feynman slashes in the numerator. After some steps, this can be seen to guarantee that the $\gamma_5$-matrices drop out.

Subsequently, taking the Dirac traces, we are left with various scalar products of four-momenta. Some of these can be immediately completed into squares,

$$K \cdot (K + P) = \frac{1}{2} \left[ \overbrace{K^2}^{\to 0} + \overbrace{(K + P)^2}^{\text{no cut}} - P^2 \right], \qquad (A.2)$$

$$Q \cdot (P + Q) = \frac{1}{2} \left[ \overbrace{Q^2}^{\text{no cut}} + \overbrace{(P + Q)^2}^{\text{no cut}} - P^2 \right]. \qquad (A.3)$$

There are further scalar products which are not immediately completed into squares, however they can also be brought into the desired form by permutations of momenta,

$$\mathbb{P}_1 \equiv \{P \to -K + Q, Q \to -P - Q\}, \quad \mathbb{P}_2 \equiv \{P \to -K + Q, Q \to K + P\}. \qquad (A.4)$$

Both leave the denominator invariant, whereas the “non-canonical” scalar products in the numerator get transferred as

$$\begin{aligned} K \cdot (P+Q) Q \cdot (K+P) &\xrightarrow{\mathbb{P}_1} K \cdot (K+P) (P+Q) \cdot Q, \\ K \cdot Q (K+P) \cdot (P+Q) &\xrightarrow{\mathbb{P}_2} K \cdot (K+P) Q \cdot (P+Q). \end{aligned} \qquad (A.5)$$

Therefore we can again make use of eqs. (A.2) and (A.3). Implementing these steps and collecting the terms together, we find that eq. (A.1) can be reduced to

$$\frac{\Gamma_{k;\nu_a}^{\text{LO}}}{G_F^2/(4k)} = \left\{ (2\sum_b) + \left[ (2\delta_{a,e} - 1 + 4s_W^2)^2 + 1 \right] + 2 \right\} \text{Im} \left\{ \frac{(-16)P^4}{(K+P)^2(P+Q)^2 Q^2} \right\}. \qquad (A.6)$$

It remains to carry out the Matsubara sums. The outer sum-integral can be extracted from the $\beta$-coefficient in eq. (2.35), yielding

$$\begin{aligned} \mathrm{Im} \left\{ \frac{P^4}{(K+P)^2 (P+Q)^2 Q^2} \right\} = \frac{1}{4\pi^2 k} \left( \underbrace{\int_{-\infty}^{0} dp_- \int_{0}^{k} dp_+}_{\text{t-channel}} - \underbrace{\int_{0}^{k} dp_- \int_{k}^{\infty} dp_+}_{\text{s-channel}} \right) \\ \times p \mathcal{P}^4 [1 - n_F(k-p^0) + n_B(p^0)] \mathrm{Im} \left\{ \frac{1}{(P+Q)^2 Q^2} \right\}. \end{aligned} \qquad (A.7)$$

The inner sum-integral can be worked out in a similar fashion, giving at first

$$\begin{aligned} \mathrm{Im} \left\{ \frac{1}{(P+Q)^2 Q^2} \right\} = \frac{1}{16\pi p} \left\{ \underbrace{\theta(p^0 - p)}_{\text{t-channel}} \int_{p_-}^{p_+} dq [1 - n_F(q) - n_F(p^0 - q)] \right. \\ \left. - 2 \underbrace{\theta(p - |p^0|)}_{\text{t-channel}} \int_{p_+}^{\infty} dq [n_F(q - p^0) - n_F(q)] \right\}. \end{aligned} \qquad (A.8)$$
Here we have made use of $n_F(-x) = 1 - n_F(x)$ in order to represent the Fermi distributions with positive arguments. Subsequently the remaining integrals can be carried out in terms of polylogarithms,

$$l_{1f}(\omega) = \ln(1 + e^{-\omega/T}), \quad l_{2f}(\omega) = \mathrm{Li}_2(-e^{-\omega/T}), \quad l_{3f}(\omega) = \mathrm{Li}_3(-e^{-\omega/T}), \quad (A.9)$$

with the latter two defined in anticipation of eq. (A.17). This converts eq. (A.8) into

$$\mathrm{Im} \left\{ \frac{1}{(P+Q)^2 Q^2} \right\} = \underbrace{\theta(p^0-p)}_{\text{s-channel}} \frac{p+2T[l_{1f}(p_+) - l_{1f}(p_-)]}{16\pi p} - 2 \underbrace{\theta(p-|p^0|)}_{\text{t-channel}} \frac{T[l_{1f}(-p_-) - l_{1f}(p_+)]}{16\pi p} . \quad (A.10)$$

Finally we can insert eq. (A.10) into eq. (A.7), obtaining a representation for $\Gamma_{k;\nu_a}^{LO}$ similar to that found in the literature (cf., e.g., ref. [20] and references therein),

$$\Gamma_{k;\nu_a}^{LO} = \left\{ (2\Sigma_b) + \left[ (2\delta_{a,e} - 1 + 4s_W^2)^2 + 1 \right] + 2 \right\} \frac{G_F^2}{16\pi^3 k^2} [2\Xi^t + \Xi^s], \quad (A.11)$$

$$\Xi^t = \int_{-\infty}^{0} dp_{-} \int_{0}^{k} dp_{+} \mathcal{P}^4 [1 - n_F(k - p^0) + n_B(p^0)] \{T[l_{1f}(-p_{-}) - l_{1f}(p_{+})]\}, \quad (A.12)$$

$$\Xi^s = \int_0^k dp_- \int_k^\infty dp_+ \mathcal{P}^4 [1 - n_F(k - p^0) + n_B(p^0)] \{p + 2T[l_{1f}(p_+) - l_{1f}(p_-)]\}. \quad (A.13)$$

Here the prefactors within the curly brackets on the first line of eq. (A.11) have been kept apart, in order to indicate their origins in eq. (2.28).

## A.2 Numerical crosscheck of the full integral representation

Let us compare the simplified representation in eq. (A.11) with the more complete one following from eq. (2.45). For this we focus on the electron-positron contribution, which can be identified through the prefactor $[(2\delta_{a,e} - 1 + 4s_W^2)^2 + 1]$ in eq. (A.11).

The leading-order electron-positron spectral functions, needed in eq. (2.45), read (the notation is from eq. (A.9))

$$\rho_V^{LO} \equiv -(2\rho_T + \rho_L)^{LO} \quad (A.14)$$

$$\stackrel{p_+ > 0}{=} \frac{\mathcal{P}^2}{4\pi p} \left\{ p\theta(p_-) + 2T[l_{1f}(p_+) - l_{1f}(|p_-|)] \right\}, \quad (A.15)$$

$$\rho_{00}^{LO} \equiv -\frac{p^2}{\mathcal{P}^2} (\rho_L)^{LO} \quad (A.16)$$

$$\stackrel{p_+ > 0}{=} \frac{1}{12\pi p} \left\{ p^3 \theta(p_-) + 12pT^2 \left[ l_{2f}(p_+) + \mathrm{sign}(p_-) l_{2f}(|p_-|) \right] + 24T^3 \left[ l_{3f}(p_+) - l_{3f}(|p_-|) \right] \right\}. \quad (A.17)$$
The limits at $p, |p^0| \ll \pi T$ can be found in eqs. (2.46) and (2.47). The values at $p_+ < 0$ could be obtained by making use of the antisymmetry in $p_0 \to -p_0$, however they are not needed in practice. Inserting these, eqs. (2.41), (2.42) and (2.45) become

$$ \begin{aligned} \Gamma_{k;\nu_a}^{\text{LO}} \supset & -[(2\delta_{a,e} - 1 + 4s_W^2)^2 + 1] \frac{G_F^2}{8\pi^2 k^2} \left( \int_{-\infty}^{0} dp_- \int_{0}^{k} dp_+ - \int_{0}^{k} dp_- \int_{k}^{\infty} dp_+ \right) \\ & \times \frac{p \mathcal{P}^2}{2} \left[ 1 - n_F(k-p^0) + n_B(p^0) \right] \\ & \times \left\{ \rho_V^{\text{LO}} \left[ 1 + \left( \frac{2k-p^0}{p} \right)^2 \right] + \frac{\rho_{00}^{\text{LO}} \mathcal{P}^2}{p^2} \left[ 1 - 3 \left( \frac{2k-p^0}{p} \right)^2 \right] \right\}. \end{aligned} \quad (A.18) $$

Even if eqs. (A.11) and (A.18) look different, with the latter containing a much more complicated integrand, both can be evaluated numerically, confirming their equivalence.

### A.3 Real part of the vector correlator

Let us decompose the real part of the vector correlator from eq. (2.22) like in eq. (2.43),

$$ \text{Re} V_{\mathcal{P}}^{\mu\nu} = \mathbb{P}_{\text{T}}^{\mu\nu} \chi_{\text{T}} + \mathbb{P}_{\text{L}}^{\mu\nu} \chi_{\text{L}}. \quad (A.19) $$

The results have divergent vacuum parts and finite thermal parts. Writing $D = 4 - 2\epsilon$ and denoting the scale parameter of the $\overline{\text{MS}}$ scheme by $\bar{\mu}^2 \equiv 4\pi\mu^2 e^{-\gamma_E}$, the vacuum parts read

$$ \chi_{\text{T}}^{\text{LO}}|_{\text{vac}} = \chi_{\text{L}}^{\text{LO}}|_{\text{vac}} = -\frac{4\mathcal{P}^2}{3} \frac{\mu^{-2\epsilon}}{(4\pi)^2} \left( \frac{1}{\epsilon} + \ln \frac{\bar{\mu}^2}{\mathcal{P}^2} + \frac{5}{3} + \mathcal{O}(\epsilon) \right). \quad (A.20) $$

As for the thermal parts, they are best represented through the same linear combinations as in eqs. (A.15) and (A.17). For these we find

$$ \chi_V^{\text{LO}}|_T \equiv -(2\chi_T + \chi_L)^{\text{LO}}|_T \quad (A.21) $$

$$ = 4 \int_q \frac{n_F(q)}{q} \left\{ -2 - \frac{\mathcal{P}^2}{4pq} \ln \left| \frac{1 - [(p+2q)/p^0]^2}{1 - [(p-2q)/p^0]^2} \right| \right\}, \quad (A.22) $$

$$ \chi_{00}^{\text{LO}}|_T \equiv -\frac{p^2}{\mathcal{P}^2} (\chi_L)^{\text{LO}}|_T \quad (A.23) $$

$$ = 4 \int_q \frac{n_F(q)}{q} \left\{ 1 + \frac{\mathcal{P}^2 + 4q^2}{8pq} \ln \left| \frac{1 - [(p+2q)/p^0]^2}{1 - [(p-2q)/p^0]^2} \right| - \frac{p^0}{2p} \ln \left| \frac{1 - [2q/(p-p^0)]^2}{1 - [2q/(p+p^0)]^2} \right| \right\}, \quad (A.24) $$

where $\int_q \equiv \int_0^\infty dq \, q^2/(2\pi^2)$. The limits at $p, |p^0| \ll \pi T$ can be found in eqs. (2.48) and (2.49).

### B Leading-logarithmic part of the NLO correction

We show here how the coefficient of the logarithm in eq. (3.5), $D_1^{(t)}$, can be derived analytically, starting from the HTL-resummed part of eq. (3.3). Physically, this corresponds to an interference term between the amplitude in eq. (4.1) and a tree-level t-channel process.
Parametrizing the interaction rate according to eq. (3.4), and considering the contribution from the domain of small spatial momenta, $0 < p < p_*$, where $eT \ll p_* \ll \{k, \pi T\}$ is an intermediate cutoff, we can write

$$ D^{(t)} \supset \frac{1}{k^3 T^4} \overbrace{\frac{1}{2} \int_0^{p_*} dp \int_{-p}^{+p} dp^0}^{\supset \int_{-\infty}^0 dp_- \int_0^k dp_+} 2p \mathcal{P}^2 \overbrace{\frac{T}{p^0}}^{\supset n_B(p^0)} (\rho_T^{\text{IR}} \chi_T^{\text{IR}} \mathcal{R}_T^* - \rho_L^{\text{IR}} \chi_L^{\text{IR}} \mathcal{R}_L^*) \overbrace{\frac{4k^2}{p^2}}^{\supset \left(\frac{2k-p^0}{p}\right)^2} \\ \stackrel{p^0=xp}{=} \frac{4}{kT^3} \int_0^{p_*} dp \, p \int_{-1}^{+1} \frac{dx}{x} \frac{(x^2-1)}{x} (\rho_T^{\text{IR}} \chi_T^{\text{IR}} \mathcal{R}_T^* - \rho_L^{\text{IR}} \chi_L^{\text{IR}} \mathcal{R}_L^*)_{p^0=xp}. \quad (B.1) $$

Here, for $n \in \{\text{T,L}\}$, eq. (2.51) can be expressed as

$$ \mathcal{R}_n^* = \text{Re} \left\{ \frac{1}{p^2(x^2-1) - e^2(\chi_n^{\text{IR}} + i\rho_n^{\text{IR}})} \right\} = \text{Re} \left\{ \frac{1}{p^2 - \Pi_n} \right\} \frac{1}{x^2-1}, \quad (B.2) $$

where we have defined

$$ \Pi_n \equiv \frac{e^2(\chi_n^{\text{IR}} + i\rho_n^{\text{IR}})}{x^2 - 1}. \quad (B.3) $$

We note that after the substitution $p^0 = xp$, the self-energies $\chi_n^{\text{IR}}$ and $\rho_n^{\text{IR}}$ from eqs. (2.46)-(2.49), and consequently $\Pi_n$ from eq. (B.3), are functions only of $x$. Therefore, we can carry out the integral over $p$,

$$ \int_0^{p_*} \frac{dp}{p^2 - \Pi_n} = \frac{1}{2} \ln \left( \frac{p_*^2 - \Pi_n}{-\Pi_n} \right) = \frac{1}{2} \left[ \ln \left( \frac{p_*^2}{e^2 T^2} \right) + \ln \left( \frac{e^2 T^2}{-\Pi_n} \right) + \mathcal{O} \left( \frac{-\Pi_n}{p_*^2} \right) \right]. \quad (B.4) $$

The first term yields

$$ D^{(t)} \supset \frac{2}{kT^3} \ln\left(\frac{p_*^2}{e^2 T^2}\right) \int_{-1}^{+1} \frac{dx}{x} (\rho_T^{\text{IR}} \chi_T^{\text{IR}} - \rho_L^{\text{IR}} \chi_L^{\text{IR}})_{p^0=xp} \\ = \frac{\pi T}{36k} \ln\left(\frac{p_*^2}{e^2 T^2}\right) \int_{-1}^{+1} dx (x^2-1) \left[4 - 3x^2 + \frac{3x(x^2-1)}{2} \ln\frac{1+x}{1-x}\right] \\ = -\frac{\pi T}{9k} \ln\left(\frac{p_*^2}{e^2 T^2}\right). \quad (B.5) $$

The coefficient of $\ln(1/e^2)$ yields $D_1^{(t)}$ from eq. (3.5).

## C Numerical results for coefficient functions

Our main results are captured by the coefficients A, B, C and D, defined in eq. (3.4). Examples of their numerical values are displayed in table 1, permitting for a precise evaluation of the full result at selected momenta. The overall pattern may be easier to appreciate from a plot, whence we show one in figure 4.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.
<table><caption>**Table 1.** The coefficients from eq. (3.4) and (3.5), as a function of the momentum $k/T$. The same results are represented as plotted curves in figure 4.</caption><thead><tr><th>k/T</th><th>A<sup>(s)</sup></th><th>A<sup>(t)</sup></th><th>B<sup>(s)</sup></th><th>B<sup>(t)</sup></th><th>C<sup>(s)</sup></th><th>C<sup>(t)</sup></th><th>D<sup>(s)</sup></th><th>D<sub>1</sub><sup>(t)</sup></th><th>D<sub>2</sub><sup>(t)</sup></th></tr></thead><tbody><tr><td>0.5</td><td>1.80721</td><td>7.31870</td><td>0.142</td><td>-0.192</td><td>3.120</td><td>4.132</td><td>0.098</td><td>-0.698</td><td>-0.654</td></tr><tr><td>1.0</td><td>1.66528</td><td>6.70427</td><td>0.061</td><td>-0.308</td><td>4.013</td><td>8.224</td><td>0.037</td><td>-0.349</td><td>-0.495</td></tr><tr><td>2.0</td><td>1.53556</td><td>6.14451</td><td>0.042</td><td>-0.372</td><td>4.741</td><td>11.56</td><td>0.011</td><td>-0.175</td><td>-0.319</td></tr><tr><td>3.0</td><td>1.50561</td><td>6.01450</td><td>0.037</td><td>-0.419</td><td>5.244</td><td>13.66</td><td>0.005</td><td>-0.116</td><td>-0.238</td></tr><tr><td>4.0</td><td>1.50956</td><td>6.03093</td><td>0.035</td><td>-0.459</td><td>5.681</td><td>15.38</td><td>0.003</td><td>-0.087</td><td>-0.192</td></tr><tr><td>5.0</td><td>1.52175</td><td>6.08366</td><td>0.034</td><td>-0.493</td><td>6.059</td><td>16.86</td><td>0.002</td><td>-0.070</td><td>-0.160</td></tr><tr><td>6.0</td><td>1.53387</td><td>6.13650</td><td>0.033</td><td>-0.522</td><td>6.381</td><td>18.11</td><td>0.002</td><td>-0.058</td><td>-0.138</td></tr><tr><td>7.0</td><td>1.54392</td><td>6.18043</td><td>0.032</td><td>-0.545</td><td>6.656</td><td>19.19</td><td>0.001</td><td>-0.050</td><td>-0.121</td></tr><tr><td>8.0</td><td>1.55185</td><td>6.21507</td><td>0.032</td><td>-0.565</td><td>6.894</td><td>20.12</td><td>0.001</td><td>-0.044</td><td>-0.108</td></tr><tr><td>9.0</td><td>1.55808</td><td>6.24214</td><td>0.032</td><td>-0.583</td><td>7.103</td><td>20.94</td><td>0.001</td><td>-0.039</td><td>-0.098</td></tr><tr><td>10.0</td><td>1.56303</td><td>6.26349</td><td>0.031</td><td>-0.598</td><td>7.288</td><td>21.67</td><td>0.001</td><td>-0.035</td><td>-0.089</td></tr><tr><td>15.0</td><td>1.57756</td><td>6.32427</td><td>0.031</td><td>-0.653</td><td>7.988</td><td>24.42</td><td>0.001</td><td>-0.023</td><td>-0.062</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The coefficients from eqs. (3.4) and (3.5), as a function of the momentum $k/T$. The precise values at a number of selected momenta can be found in table 1, and the coefficient $D_1^{(t)}$ is given analytically in eq. (3.5). The coefficients $A$ and $C$ have been multiplied by $1/(6\pi^2) \approx 0.0169$, so that their magnitudes can be meaningfully compared with those of $B$ and $D$, cf. eq. (3.4).</figcaption></figure>
References

[1] K. Akita and M. Yamaguchi, *A precision calculation of relic neutrino decoupling*, JCAP 08
(2020) 012 [arXiv:2005.07047] [INSPIRE].

[2] J. Froustey, C. Pitrou and M.C. Volpe, *Neutrino decoupling including flavour oscillations and primordial nucleosynthesis*, JCAP 12 (2020) 015 [arXiv:2008.01074] [INSPIRE].

[3] J.J. Bennett et al., *Towards a precision calculation of $N_{\text{eff}}$ in the Standard Model. Part II. Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED*, JCAP 04 (2021) 073 [arXiv:2012.02726] [INSPIRE].

[4] M. Cielo, M. Escudero, G. Mangano and O. Pisanti, $N_{\text{eff}}$ in the Standard Model at NLO is 3.043, Phys. Rev. D 108 (2023) L121301 [arXiv:2306.05460] [INSPIRE].

[5] S. Esposito et al., *Neutrino energy loss rate in a stellar plasma*, Nucl. Phys. B 658 (2003) 217
[astro-ph/0301438] [INSPIRE].

[6] G. Sigl and G. Raffelt, *General kinetic description of relativistic mixed neutrinos*, Nucl. Phys. B 406 (1993) 423 [INSPIRE].

[7] D. Nötzold and G. Raffelt, *Neutrino dispersion at finite temperature and density*, Nucl. Phys. B 307 (1988) 924 [INSPIRE].

[8] S. Sarantakos, A. Sirlin and W.J. Marciano, *Radiative corrections to neutrino-lepton scattering in the $SU(2)_L \otimes U(1)$ theory*, Nucl. Phys. B 217 (1983) 84 [INSPIRE].

[9] R.J. Hill and O. Tomalak, *On the effective theory of neutrino-electron and neutrino-quark interactions*, Phys. Lett. B 805 (2020) 135466 [arXiv:1911.01493] [INSPIRE].

[10] O. Tomalak and R.J. Hill, *Theory of elastic neutrino-electron scattering*, Phys. Rev. D 101
(2020) 033006 [arXiv:1907.03379] [INSPIRE].

[11] M. Escudero Abenza, *Precision early universe thermodynamics made simple: $N_{\text{eff}}$ and neutrino decoupling in the Standard Model and beyond*, JCAP 05 (2020) 048 [arXiv:2001.04466] [INSPIRE].

[12] H.A. Weldon, *Effective fermion masses of order $gT$ in high-temperature gauge theories with exact chiral invariance*, Phys. Rev. D 26 (1982) 2789 [INSPIRE].

[13] E. Braaten and R.D. Pisarski, *Soft amplitudes in hot gauge theories: A general analysis*, Nucl. Phys. B 337 (1990) 569 [INSPIRE].

[14] J.C. Taylor and S.M.H. Wong, *The effective action of hard thermal loops in QCD*, Nucl. Phys. B 346 (1990) 115 [INSPIRE].

[15] G. Jackson, *Two-loop thermal spectral functions with general kinematics*, Phys. Rev. D 100
(2019) 116019 [arXiv:1910.07552] [INSPIRE].

[16] G. Jackson, *Numerical code for master integrals for thermal spectral functions*,
https://zenodo.org/record/3478143.

[17] A. Salvio, P. Lodone and A. Strumia, *Towards leptogenesis at NLO: the right-handed neutrino interaction rate*, JHEP 08 (2011) 116 [arXiv:1106.2814] [INSPIRE].

[18] G. Jackson and M. Laine, *Efficient numerical integration of thermal interaction rates*, JHEP 09
(2021) 125 [arXiv:2107.07132] [INSPIRE].

[19] M. Drewes et al., *Towards a precision calculation of $N_{\text{eff}}$ in the Standard Model III: Improved estimate of NLO corrections to the collision integral*, arXiv:2402.18481 [INSPIRE].

[20] J. Ghiglieri and M. Laine, *Neutrino dynamics below the electroweak crossover*, JCAP 07 (2016) 015 [arXiv:1605.07720] [INSPIRE].