QED corrections to the thermal neutrino interaction rate

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.


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 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 (SM) eff ≈ 3.043 ± 0.001 [1][2][3][4].Were the observed value to eventually differ from N (SM)  eff , this could be an indication of physics beyond the Standard Model.The latest ingredient included in the Standard Model prediction of N eff originates from QED corrections to the rate at which neutrinos interact with the electron-positron plasma at the time of their decoupling, T ∼ (1...3)MeV [4].This result was based on previous work in a different (astrophysical) environment [5].Concretely, the observable evaluated is Q ≡ δρ δt e + e − →νν(γ) ≡ (ǫ e + + ǫ e − ) |M| 2 e + e − →νν(γ) phase space , ( where ρ 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.sec.2.1).
Our computation entails some practical differences compared with refs.[4,5].First of all, e + e − → νν(γ) represents only a subclass of the processes that we consider (cf.fig. 2 on page 14 for the full set); we add reactions like eν → eν including its 1-loop virtual corrections, or eγ → eνν, or a logarithm ∼ ln(T /m e ) that can be associated with a loop-induced ννγ 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 ∼ ln(1/α 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. sec.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 sec.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 sec.3, where our NLO results are given and elaborated upon.Conclusions are offered in sec.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.

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).

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 ≡ |k|, where w † a and w b are creation and annihilation operators, respectively, and a, b are flavour indices. 1 The phase factor φ 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 Here the time derivative reads (...) ≡ (∂ t − Hk∂ k )(...), where H is the Hubble rate; ∆E k is a matrix involving energy differences; and ρ 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 Γ k . 2 We denote its diagonal components (in the interaction basis) by Γ k;ν 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 disentagle the full interaction rate into partial ones, 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].
2 Another important effect are medium modifications to ∆E k [7].They are of O(G F ), where where i labels the different ensembles.We employ this representation in sec.3, however assume the presence of a universal temperature in all the ensembles.

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 Minkowskian (t) and Euclidean (τ ) time coordinates are formally related via the Wick rotation τ ↔ it, whereas for the corresponding momenta we have k n ↔ −ik 0 , where the Matsubara frequencies read k n = 2nπT for bosons, and k n = (2n + 1)πT for fermions, with n ∈ .Whenever in Minkowskian spacetime, the signature (+−−−) is employed.Minkowskian fourmomenta are denoted by K = (k 0 , k), Euclidean ones by K = (k n , k).We use an implicit notation where the symbol used indicates whether analytic continuation has been carried out, for example Euclidean Dirac matrices, identified by the fact that all indices are down, are defined as It follows that The Fermi coupling reads at tree level, or leading order (LO), and the Weinberg angle is parametrized by 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 2 W need to be revisited later on.The numerical values used for their MS renormalized versions at the scale μ = m e are specified in the next section.(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 fig. 2 (including thermal corrections).

Interactions according to a Fermi effective theory
At temperatures T ∼ MeV, all hadrons, as well as the µ and τ 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 ) 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 νµ is the QED field strength and e is the electromagnetic coupling, is generated at 1-loop level, through the diagrams shown in fig. 1 (cf., e.g., refs.[8,9]).The coefficient C a includes a computable divergence and anomalous dimension, but its absolute value cannot be determined perturbatively, as it gets a contribution from hadronic effects through the Zγ 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, This structure is antisymmetric in j ↔ l.This is compensated for by another minus sign from the anticommutation of fermionic fields, whereby eq.(2.10) can finally be turned into This is the form that we employ for our practical computations.
Let us now return to the operator δ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 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γ 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 4fermion 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], where a R ≡ (1 + γ 5 )/2 and a L ≡ (1 − γ 5 )/2 are the right and left-handed projectors.The coefficients g eR , ..., g (µ, τ )L correspond to a notation employed in the literature (cf., e.g., ref. [11]), where often no distinction is made between a = µ and a = τ .After renormalization, the finite parts of the coefficients g aR and g aL become running parameters.We fix them through their values at the MS scheme renormalization scale μ = 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 √ 2G F g aR,L from eq. (2.16) agree with the coefficients c listed in table 4 of ref. [9].Relevant for us is the case that only the charged flavour ℓ ′ = ℓ e is light.In fact, the comparison is overconstrained, as we have parametrized six independent coefficients c νaℓ e R,L , a ∈ {e, µ, τ }, 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 2386, and (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 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 = e, A virtual electron-positron loop couples to a photon only through g aV .Therefore the flavours a = 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 ν e and ν a =e .

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 ∆ (full) If we furthermore define the electron currents and their fermion-line connected ("c") correlation functions then the O(G 2 F ) contribution to the self-energy becomes is the tree-level neutrino propagator, ∆ −1 P ;ρσ 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 fig.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 ∼ πT .If we consider temperatures T > ∼ 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 γ 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 − γ 5 )'s of the second term of eq.(2.26) can be transported into the projectors in eq. ( 2.25), we thereby obtain lim (2.28)

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 Π P .We write it in a spectral representation, where p n is bosonic.Then we are faced with the structure where k n and r n are fermionic, and ǫ 2 kp ≡ (k + p) 2 .We may now write δ 0,k n +p n −r n = T 1/T 0 dτ e i(k n +p n −r n )τ , and carry out the sums over p n and r n .Furthermore the integral over τ is doable, yielding where n B and n F denote the Bose and Fermi distributions, respectively.As a next step, we continue to Minkowskian frequency, ik n → k 0 + i0 + , and take the imaginary part.We also substitute p 0 → −p 0 , p → −p, and note that bosonic spectral functions are odd under this substitution.Finally, ǫ kp is redefined as ǫ kp ≡ |k − p|, and we recall that n (2.32) Subsequently, we may integrate over the angles of p, removing the Dirac-δ's.At this point we also go to the light cone, k 0 = k, relevant for on-shell neutrinos.The first Dirac-δ can be seen to be realized for |p 0 | < p, i.e. it corresponds to a t-channel process.The second Dirac-δ is realized for p 0 > p, and corresponds to an s-channel process.For both channels, the Dirac-δ's imply (2.33) eq. ( 2.32) can be turned into (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 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).

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 sumintegrals have been carried out, going over to Minkowskian signature (whereby i / K → / K ), and redundantly displaying chiral projectors in accordance with eq.(2.20), the result has the form [12] Σ where u ≡ (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, (2.37) The real part b r , which gives a medium modification to the effective Hamiltonian ∆E k appearing in eq.(2.2), is of O(G F ) [7].The imaginary part, which is of O(G 2 F ) and interests us here, can be projected out as 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 where we noted that the part involving γ 5 drops out, if Ë a P ;µν is symmetric in µ ↔ ν, and that the parts proportional to P µ and P ν drop out, if Ë a P ;µν is transverse.
Subsequently, we make use of eq. ( 2.35), with Ë a P ;µν playing the role of Π P .Inserting eq.(2.33), this leads to To complete the transition to Minkowskian spacetime, we recall the rules from eq. (2.7), resulting in

.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 sec.3, let us focus on Im Ë a;µν (2.42) The vector channel spectral function can be decomposed as where the projectors have been defined as (2.44) Inserting this representation, and making use of eq. ( 2.33) in order to resolve the scalar products k • p originating from K µ K ν È µν T , the contraction in eq. ( 2.41) takes the form (2.45)This representation will be employed in sec.3, and can also be crosschecked numerically, as demonstrated in sec.A.2.

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, ∆ −1 P ;ρσ in eq.(2.28), is "soft", |p 0 | ≤ p ≪ πT .The technical reason is that the vector correlators turn into "Hard Thermal Loops" (HTLs) in this domain, V LO P ∼ T 2 .Then they do not regulate that IR regime any more, as would be the case in vacuum, where they behave as ∼ P 2 .In addition, soft photons are Bose-enhanced, with n B (p 0 ) ≈ T /p 0 .
HTLs appear not only in the vertices V LO P , but they also influence the self-energies of the photons.When all effects scaling as ∼ 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 LO P .The imaginary parts are obtained from eqs. (A.15) and (A.17),

.47)
The real parts can be derived from eqs. (A.22) and (A.24), With the help of eqs.(2.46)-(2.49),we can construct the resummed photon propagator, denoted by ∆ −1 * P ;µν .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 R * T,L and I * T,L , (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 R * T,L , but not from I * T,L , since the latter leads to a contribution of O(e 4 ). 3 In the s-channel, we set R * T,L → R T,L ≡ 1/P 2 , since the integral is IR-finite.Furthermore, we omit the effect from I * T,L in the s-channel as well.It would lead to the processes illustrated in set (h) of fig. 2. These correspond to the famous plasmon contributions (cf., e.g., ref. [5]).However, they are formally of O(e 3 ) rather than O(e 2 ): apart from the vertices, the result of a 1 ↔ 2 reaction is phase-space suppressed by the plasmon mass ∼ eT .
After the inclusion of HTL resummation in the t-channel, the contribution from soft momentum exchange is finite.However, it is logaritmically 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).

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, (a) the same neutrino flavour, different neutrino flavours, and electrons and positrons, yield the contributions Γ respectively, where p ± are defined according to eq. (2.34); P 2 ≡ (p 0 ) 2 − p 2 = 4p + p − ; and n F and n B denote the Fermi and Bose distribution functions.The functions ρ LO,NLO T,L originate from the imaginary part of the vector current correlator in eq.(2.22), and are discussed below; the functions χ LO T,L 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 ρ NLO T,L .These represent the diagrams in fig.2(c)-(f): all of these processes, and the interference terms between them, are contained in the evaluations of ρ NLO T,L 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 2 C F → e 2 ≡ 4πα 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 fig.2(g), once these interfere with LO contributions.In the notation of the last term of fig.2(b), the open box is replaced by a closed loop here; it is the closed loop which produces the functions χ LO T,L .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 R * T,L , as explained in sec.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 fig.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 χ LO T,L , 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γ bubble, cf.sec.2.3.
Apart from the cancellation of the divergence, the concrete role of C a is to replace the renormalization scale μ 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 The couplings are evaluated at the scale μ = 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) = − πT 9k .
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 fig. 3.

Conclusions
The general context of our investigation is the observable known as N 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 eff , ref. [4] recently carried out an approximate NLO computation of the energy transfer rate defined in eq.(1.1).They found a ∼ −5% QED correction to the energy transfer rate in the temperature range T ∈ (1...3)MeV.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 eff [1][2][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 momentumdependent 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.fig.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 fig.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 fig. 3. We find a ∼ −(0...2)% contribution, dominated by the t-channel diagrams in fig.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 fig.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 where "V" indicates that the QED vertices are purely vectorlike.As discussed around eq. (2.19), the flavours a = 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 fig.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 ∼ 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 nonequilibrium 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γ mixing.Our study made use of the range in eq.(2.17), extracted from refs.[9,10], however the error, visible as bands in fig.3, could conceivably be reduced through further theoretical and phenomenological work.the full integral representation in eq.(2.45) (cf.sec.A.2). Finally, we specify the real part of the correlator from eq. (2.22), as this is needed in eq.(3.3) (cf.sec.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, 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 )) The first simplification originates from the substitution Q → −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 γ 5 -matrices drop out.Subsequently, taking the Dirac traces, we are left with various scalar products of fourmomenta.Some of these can be immediately completed into squares, 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, .4)Both leave the denominator invariant, whereas the "non-canonical" scalar products in the numerator get transferred as Here, for n ∈ {T,L}, eq. ( 2.51) can be expressed as where we have defined We note that after the substitution p 0 = x p, the self-energies χ IR n and ρ IR n from eqs. (2.46)-(2.49),and consequently Π n from eq. (B.3), are functions only of x.Therefore, we can carry out the integral over p,

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 fig. 4.  values at a number of selected momenta can be found in table 1, and the coefficient 1 is given analytically in eq.(3.5).The coefficients A and C have been multiplied by 1/(6π 2 ) ≈ 0.0169, so that their magnitudes can be meaningfully compared with those of B and D, cf.eq.(3.4).

Figure 2 :
Figure2: (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).

a
≠ e, k = T a ≠ e, k = 3T a ≠ e, k = 10T percentual effect, a ≠ e

Figure 3 :
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 α 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 = e (right).Bottom row: dependence on T /MeV for fixed k/T , for a = e (left) and a = e (right).

Figure 4 :
Figure 4: The coefficients from eqs. (3.4) and (3.5), as a function of the momentum k/T .The precise