Complete leading-order standard model corrections to quantum leptogenesis

Thermal leptogenesis, in the framework of the standard model with three additional heavy Majorana neutrinos, provides an attractive scenario to explain the observed baryon asymmetry in the universe. It is based on the out-of-equilibrium decay of Majorana neutrinos in a thermal bath of standard model particles, which in a fully quantum field theoretical formalism is obtained by solving Kadanoff-Baym equations. So far, the leading two-loop contributions from leptons and Higgs particles are included, but not yet gauge corrections. These enter at three-loop level but, in certain kinematical regimes, require a resummation to infinite loop order for a result to leading order in the gauge coupling. In this work, we apply such a resummation to the calculation of the lepton number density. The full result for the simplest"vanilla leptogenesis"scenario is by $\mathcal{O}(1)$ increased compared to that of quantum Boltzmann equations, and for the first time permits an estimate of all theoretical uncertainties. This step completes the quantum theory of leptogenesis and forms the basis for quantitative evaluations, as well as extensions to other scenarios.


Introduction
The origin of the observed baryon asymmetry in the universe is an as yet unsolved problem of physics. Although the standard model of particle physics (SM) features baryon plus lepton number (B + L) and CP -violating processes before the electroweak phase transition [1], the amount of CP violation is too small to arrive at the observed asymmetry [2]. Moreover, for the experimentally observed Higgs mass of ∼ 125 GeV [3], the electroweak phase transition is merely an analytic crossover [4][5][6], whereas a sufficiently strong first-order transition is required to provide a departure from equilibrium for baryogenesis. Viable models for baryogenesis are thus based on extensions of the SM, with, e.g., additional particles on a GUT scale, which then generate a finite asymmetry via CP -violating out-of-equilibrium decays. A particularly attractive model in this context is thermal leptogenesis [7]. In this case, the SM is extended by three additional right-handed Majorana neutrinos, with Yukawa couplings to the SM Higgs field and left-handed leptons. In the standard scenario [8], these Majorana neutrinos are produced in the early universe at temperatures beyond their mass scale, T > M . When the temperature drops below their mass, they decay out-ofequilibrium violating CP , such that a lepton asymmetry is generated. Additional washout processes diminish this asymmetry, until they fall out of equilibrium and a finite lepton asymmetry is frozen in. Any lepton asymmetry is partially converted into a baryon asymmetry via SM sphaleron transitions, which violate B + L while keeping B − L constant [9,10]. Besides giving robust predictions for baryogenesis, this simple extension of the SM would also explain the smallness of the light neutrino masses via the see-saw mechanism [11]. For recent reviews, also including other scenarios for leptogenesis, see [12][13][14][15][16][17].
Since leptogenesis is a non-equilibrium process, it is frequently studied using Boltzmann equations, which are inherently classical. For the collision terms, zero temperature matrix elements are used, supplied with thermal momentum distribution functions which introduce quantum effects. However, also off-shell and memory effects become important for non-equilibrium processes and have to be treated appropriately. Further conceptual difficulties arise in gauge theories, where gauge boson numbers are neither gauge invariant nor conserved. The problem of including quantum effects by hand into a classical framework is avoided by a formal, quantum field theoretical treatment based on Kadanoff-Baym equations, which contain everything once computed to sufficient depth in perturbation theory. In [18], quantum corrections to the Majorana neutrino reaction rates due to Higgs and lepton loops have been included, and a detailed comparison with the solutions of Boltzmann equations was given. Qualitative agreement was achieved by additionally introducing thermal damping widths for Higgs and lepton propagators, which are expected to be caused by interactions with the weak gauge bosons and may significantly modify the result at a quantitative level [19].
In this work we complete the quantum mechanical treatment of leptogenesis by calculating gauge corrections systematically. First entering at three-loop level, their evaluation is complicated by the dynamical generation of scales, which necessitate infinite resummations in order to complete even a leading-order result in the gauge coupling. Such resummations have been performed for the Majorana neutrino self-energy and the associated production rate [20]. Here we extend these calculations to the lepton number density in a first complete quantum mechanical calculation of leptogenesis. Our formulation using Kadanoff-Baym equations [18] differs from another field-theoretical approach, in which linearised rate equations are obtained by exploiting hierarchies of scales in the non-relativistic regime, T < M [21][22][23].
In order to render this paper self-contained, sections 2 and 3.1 recall the formalism and previous results, which form the basis for our calculation. In section 3.2, we devise a resummation technique to systematically include gauge corrections to the lepton number density, which then is complete to leading order in all SM corrections. In section 4, some approximations are motivated which allow for a simple evaluation. Finally, we compare to previous results in the literature and discuss the significance of gauge corrections.

The leptogenesis scenario
Starting from the SM, a standard approach to leptogenesis is to add three additional righthanded Majorana neutrinos (which are electroweak singlets) to the Lagrangian L SM , Here the flavour indices i, j = 1, 2, 3 count the families in the gauge eigenbasis, and we denote ν c R = Cν T R , with C = iγ 2 γ 0 the charge conjugation matrix, andφ = iσ 2 φ * . The additional neutrinos are weakly coupled to SM Higgs fields φ and massless left-handed leptons l Li via Yukawa couplings λ ij . Without loss of generality, we consider here the simplest scenario, also termed "vanilla leptogenesis", in which the final asymmetry is assumed to be independent of its flavour composition and, in a diagonal mass basis, the Majorana masses M i are hierarchically ordered, M i>1 M 1 := M . For a review including details, viable parameter ranges and other scenarios, see [8]. In this case, the two heavier neutrinos can be integrated out, leading to the effective Lagrangian for the lightest neutrino with an effective vertex for lepton and Higgs fields and a combination of Majorana neutrino couplings, respectively, The heavy neutrino mass is in the range 10 9 GeV M 10 15 GeV [8], and its coupling is assumed very weak compared to SM couplings g SM , λ g SM . Explicitly, the relevant SM couplings for our work are the SU (2) and U (1) gauge group couplings g and g , the Yukawa coupling of the top quark h t , and the Higgs self coupling λ φ , i.e. g SM ∈ {g, g , h t , λ φ }.
During the generation of the lepton asymmetry, the heat bath is kept in thermal equilibrium via SM interactions. These act on a time scale τ SM ∼ 1/(g 2 SM T ), which is much shorter than the equilibration time of the heavy neutrino, τ N ∼ 1/(λ 2 M ) for T M , which thus is out of equilibrium. During this first stage, heavy Majorana neutrinos are produced from zero initial abundance by scattering of leptons and Higgs bosons. This out-of-equilibrium process generates an initial lepton asymmetry, which is later diminished by washout processes. At temperatures T ∼ M , part of the Majorana neutrinos decay out of equilibrium to produce a final lepton asymmetry, or B − L > 0, which is of the same size as the initial one [24]. 1 Since one can neglect the Hubble expansion, i.e. assume a Figure 2: Keldysh contour C in the complex time plane for the calculation of nonequilibrium correlation functions. It starts at some initial time x 0 = t i + i✏ and runs parallel to the real axis up to some final time t f + i✏ and returns to t i i✏. Physical correlation functions are then obtained in the limit t f ! 1 and ✏ ! 0.
contour C in the complex time plane shown in Fig. 1, leading to the following results for massless left-handed leptons and Higgs fields All other propagators can be obtained from linear combinations of the propagators given above.
The nonequilibrium case is much more involved. The calculation of nonequilibrium correlation functions requires physical time making the use of the real-time formalism necessary. The calculation of nonequilibrium propagators are initial value problems. The S-Matrix approach familiar from zero temperature quantum field theory fails as it is not possible to send the initial time t i ! 1 and the final time to t f ! +1 to obtain asymptotic 'in' and 'out' states. For a nonequilibrium system prepared in an initial state at initial time t i , no information on the system at later times is known. As a result, a time integration contour C for nonequilibrium correlation functions can only start and end at the same initial time t i . The time integration contour satisfying this condition is given by the Keldysh contour Fig. 2.
The nonequilibrium correlation functions can be calculated along this time contour C by making use of the so-called Schwinger-Keldysh formalism. When defining the Green's function of the Majorana neutrino with the time argument living on C, we have to impose the correct path ordering with the ✓-function enforcing time ordering along the path. Because the Green's function is calculated along the Keldysh contour it depends on the location of the time arguments on the integration contour C. If both time arguments lie 7 Re t number of each state. For many calculations in statistical quantum field theory, it is useful to define another set of propagators known as Wightman functions Both types of propagators are connected by Analogously, it is possible to define these kinds of propagators for the massless lefthanded leptons S ± Lij and S >,< Lij appearing in the Lagrangian (2.2) of our e↵ective model. The index L denotes the projection to left-handed fields S ± L = P L S ± with the projector P L = (1 5 )/2. Keep in mind that for fermions, all commutators have to replaced by anticommutators and vice versa. The equilibrium SM propagators are well known, see e.g. [18]. In most cases, it is convenient to consider the correlation functions in momentum space by performing a Fourier transformation For later purpose, remember that propagators in momentum space evaluated in statistical quantum field theory fullfill the Kubo-Martin-Schwinger (KMS) relations [18] < q (!) = e ! > q (!) , (2.13) with the inverse temperature . The equilibrium Green's functions can be evaluated in the real-time formalism using the 6 Figure 1. Non-equilibrium (left) and equilibrium (right) contours C in the complex time plane.
constant temperature, during the first stage, we will for technical convenience follow [18] in calculating this initial lepton asymmetry at fixed T starting from zero initial Majorana neutrino abundance. In reality, any lepton asymmetry gets continually converted to a baryon asymmetry by sphaleron processes, which are in equilibrium for temperatures above the electroweak phase transition. Our first calculation of gauge corrections focuses on the lepton aymmetry only and neglects these and other spectator processes, which may affect the final baryon asymmetry by ∼ 50% [25][26][27].

Correlation functions and Kadanoff-Baym equations
In a quantum field theoretical description of leptogenesis, the desired information is encoded in the non-equilibrium correlation function of the Majorana neutrino, which is weakly coupled to a thermal bath of SM particles. This requires the real time formalism of thermal quantum field theory [28], in which two-point Green functions live on a contour C in the complex time plane, whose form depends on whether the field is in equilibrium or not, cf. figure 1. Our formulation and notation follows [18], where details of the calculation can be found. Green functions can be expressed in terms of Wightman functions, which for a real scalar field read with the Heaviside θ C -function enforcing path ordering along the time contour C. These in turn are related to the spectral function and statistical propagator, respectively, Analogous quantities are defined for the Green functions of leptons and the heavy Majorana neutrino, with the role of (anti-)commutators exchanged for fermions, In the following, Π C denotes the self-energy of the leptons and Σ C that of the Majorana neutrino. These are equally defined on the time contour and can be similarly decomposed into Wightman functions, According to the physics scenario given in section 2.1, the non-equilibrium Green function of the Majorana neutrino has to be evaluated on the Schwinger-Keldysh contour, figure 1 (left), while the lepton and Higgs fields are in the thermal bath of SM particles described by the equilibrium contour, figure 1 (right), with inverse temperature β = 1/T . The full time evolution of the Majorana neutrino correlation function is given by the coupled set of Kadanoff-Baym equations, with the spatial Fourier transforms assuming spatial homogeneity ( These equations are exact and in particular contain all quantum effects. Interactions with the plasma are automatically included via the self-energies Σ ± (x 1 , x 2 ), as can be seen via generalised cutting rules. The leading-order contribution to the self-energy is given by the diagram in figure 2. An analytic solution can be obtained if the Yukawa coupling of the Majorana neutrino to the thermal bath is weak, so that back-reactions, being of higher order and additionally suppressed by the large number of degrees in freedom of the bath, can be neglected, and a narrow-width approximation Γ ∼ λ 2 T T is justified [18]. 2 Since the particles in the loop are thermal, the heavy neutrino spectral function is time-translation invariant, Using a zero abundance of the Majorana neutrinos at initial time t i = 0 and denoting with ω p = p 2 + M 2 and the Fermi-Dirac distribution function f F (ω p ).
In general, Γ p is the Majorana neutrino decay width resulting from its self-energy via Due to its appearance in the Majorana neutrino propagators, we will mostly need the on-shell value Γ p (ω p ), which satisfies (2.20)

The lepton number matrix
A net lepton number density is obtained from the flavour-diagonal lepton number current The lepton number matrix is defined as the (spatial) Fourier transform of the current's zero component, It has to be evaluated calculating the leading order correction to the lepton statistical propagator, contributing the necessary CP -violation. The corresponding Feynman diagrams are shown in figure 3, giving [18] with the lepton self-energy from the first graph 3 in figure 3, Here we introduced the notation p ≡ d 3 p/(2π) 3 andG p is the scalar non-equilibrium part of the Majorana neutrino propagator, connected to the full propagator via projections with the solutionG Note that the equilibrium part of the neutrino propagator does not contribute. We do not need and hence do not give the explicit form of the Higgs and lepton propagators ∆ ≶ q and S ≶ k , which can however be found in [18]. Using the relations between the different kinds of propagators [28] and identifying it is possible to express the lepton number matrix in terms of the retarded Majorana neutrino self-energy Σ ret . Using Im(λ * i1 (ηλ * ) i1 ) = 16π ii /(3M ), with ii parametrising the strength of CP violation [29], we have . This expression contains quantum mechanical off-shell and memory effects. However, it does not yet contain corrections from interactions with SM gauge bosons, which introduce thermal damping and the associated widths for the Higgs and lepton fields. As explained in a detailed discussion in [18], this fact is responsible for the apparent difference between quantum Boltzmann solutions and the present solution of Kadanoff-Baym equations. Moreover, for T M these damping widths are of order γ l,φ ∼ g 2 SM T λ 2 T ∼ Γ, i.e. larger than the heavy neutrino rate, and hence can modify the result significantly [19]. Any quantitative theory of quantum leptogenesis therefore requires a systematic inclusion of these effects. Figure 3. Two-loop corrections to the lepton self-energy Π ± leading to non-zero lepton number densities.

Gauge corrections to leptogenesis
The inclusion of gauge corrections to thermal leptogenesis is complicated by the fact that it happens before the electroweak transition, i.e. gauge bosons are massless and the nonabelian weak dynamics are QCD-like [30][31][32]. This in particular implies the dynamical generation of mass scales ∼ T, gT, g 2 T , of which the ultra-soft ∼ g 2 T is entirely nonperturbative. For processes with gauge particles on the soft scale ∼ gT , hard thermal loops (HTL) contribute the same power in the coupling constant to the self-energy to any loop order, and hence have to be resummed when counting orders in the gauge couplings [28,33]. Here we consider fermion processes, whose external momenta in a plasma are typically hard, k ∼ T . Nevertheless, for hard momenta near the light cone, k 2 ∼ g 2 SM T 2 , loop momenta collinear with the external ones require a similar resummation of collinear thermal loops (CTL) [34]. For power counting purposes, all SM couplings are considered as parametrically the same.

CTL resummation of the Majorana neutrino self-energy
As a first step, we need the gauge corrections to the Majorana neutrino self-energy and its proper resummation shown in figure 4 to determine the Majorana neutrino production rate [20]. One effect of interactions with the plasma is to modify the dispersion relations of the thermalised leptons and Higgs particles by the asymptotic thermal masses 4 Next, we consider nearly collinear momenta k close to the light cone, pointing approximately in the three-direction of the light-like four-vector v = (1, v), v 2 = 1. Three-momenta are then specified by components parallel and perpendicular to v, The leptonic and Higgs loop momenta thus satisfy and thermal masses are important. Further, the gauge loop momenta are considered soft, In this kinematical region the diagrams without and with gauge corrections, figures 2 and 4, respectively, are parametrically the same, where we used the power counting rules from [20,34]. Thus, adding soft internal gauge boson lines does not increase the order in the gauge coupling and necessitates a resummation of the entire ladder to infinite loop order, figure 4, while contributions from crossed ladder rungs and vertex corrections to loop particles are further suppressed [20,34]. The problem corresponds to the non-abelian version of the Landau-Pomeranchuk-Migdal effect and has been solved for gluon radiation in a QCD plasma [35]. Its application to the heavy Majorana neutrino was treated in [20], where the ladder-resummed self-energy, including thermals widths of the Higgs and leptons, is written as with d(r) = 2 the dimension of the SU (2) fundamental representation and the combination of Fermi-Dirac and Bose-Einstein distributions The functions ψ and f = (f 1 , f 2 ) are defined as solutions of the integral equations with the kernel The kernel is obtained from the gauge field propagator and contains the ladder resummation, with the SU (2) Casimir operator C 2 (r) = 3/4, the hypercharge y l = −1/2, and the HTL-resummed Debye masses [36] The quantity (p, k) is given in a frame, where v p and p ⊥ = 0, as (3.14) The integral eqs. (3.7), (3.8) are best solved in Fourier space, , and the equations are converted into ordinary differential equations,  where the c 2 are asymptotic solutions The off-diagonal elements of eq. (3.19) receive contributions vanishing in the limit b → 0, and O(1/b) divergent contributions. These are removed by renormalization, since they appear in the temperature independent part of the self-energy, while the temperature dependent part is given by the function F(p , k ). For later use in the lepton number density we define the functions Because of the symmetry between right-and left-handed Dirac components, we have In order to quantify the effect of the ladder resummation, it is useful to also analyse eqs. (3.7), (3.8) with C(|u ⊥ |) = 0, giving a Majorana neutrino self-energy which is only partly corrected by asymptotic masses. The equations in this limit can be solved trivially and one obtains after k ⊥ -integration The integration boundaries in this case are given by In figure 5 the fully corrected σ h,ψ (ω p , p ) are plotted for M = 10 10 GeV, T = 10 11 GeV, and compared to the σ m∞ h,ψ (ω p , p ). The latter are significantly smaller than their fully resummed counterparts in most of the parameter space considered. Note that above a certain p , the σ m∞ h,ψ (ω p , p ) vanish kinematically. For the SM corrected Majorana decay width we start with eq. (2.19). In Weyl basis it is possible to simplify further,

Gauge corrections to lepton self-energies
We now have everything at hand to also include gauge corrections to the CP -violating diagrams in figure 3, which enter the lepton number matrix. Adding one internal gauge boson line modifies those diagrams by either self-energy or vertex corrections of the internal lines, the resulting three-loop diagrams are collected in appendix A. From the last section it is clear that in the collinear kinematical region this will not suffice for a result to leading order in the gauge coupling. Instead, infinitely many insertions have to be resummed. In order to achieve this, it is useful to re-introduce the heavier Majorana neutrinos N 2,3 explicitly, which had been integrated out for the simplest leptogenesis scenario in eq. (2.2). The following procedure, depicted in figure 6, thus readily applies to more general scenarios. This step (1) brings the two-loop diagram to a symmetric form, figure 6. Next we close the external fermion lines (2), which implies an integration over lepton momenta and transforms the lepton number matrix into the lepton number density, leading to a cylinder diagram. Its upper and lower parts now are manifest analogues of the diagram discussed in the last subsection and thus can be resummed accordingly (3).
The resulting diagram will be referred to as resummed cylinder diagram in the following. It is instructive to check by direct comparison with appendix A that indeed all three-loop diagrams shown there are included and resummed.

Lepton number matrix with complete SM corrections
Collecting all results from the previous sections it is possible to give an expression for the SM corrected lepton number matrix. Starting with eq. (2.29) we first carry out the q and q momentum integrals by using the momentum conserving δ-functions, so that the trace in the integrand takes the form = tr ImΣ ret R,p (ω 21 )ImΣ ret R,p (ω 23 ) , (3.34) where we have used the relation between the left-and right-handed self-energy (3.24).
Finally we use the resummed cylinder diagram figure 6 to insert the resummed Majorana neutrino self-energy in the equation above. Then, the trace can be carried out using eq. (3.19), tr ImΣ ret R,p (ω 21 )ImΣ ret R,p (ω 23 ) = σ ψ (ω 21 , p )σ ψ (ω 23 , p ) + σ h (ω 21 , p )σ h (ω 23 , p ) . (3.35) Altogether, we arrive at the following result of the corrected lepton number density: This expression contains all SM corrections to leading order in the respective couplings and for the first time gives a complete description of quantum leptogenesis.

Evaluation of the lepton number density
It now remains to numerically evaluate our final result eq. (3.36). As we shall see shortly, the time integrations can be performed analytically. However, the remaining three-dimensional integral has an integrand with rapidly oscillating parts, rendering an accurate numerical evaluation difficult. Moreover, σ h and σ ψ are specified by the asymptotic solutions of the differential eqs. (3.16) and (3.17), which have to be solved numerically for each call of those functions. A brute force numerical result would thus require a significant amount of computing power.
Here we adopt a different approach. Since the final result is accurate to leading order in the couplings only, we limit the integration range to those regions, which give the parametrically leading contribution. This permits us to do two more integrations analytically, leaving only the p -integral as a final numerical task. Since we are left with only one momentum variable, we simplify the notation p → p in this section.

The time integrals
We begin by factoring out the time integrals, The resulting function, can be integrated explicitly with the help of Mathematica, and the result is spelled out in appendix B. We now use this expression to estimate in which domain of frequencies it gives parametrically dominant contributions to the integral. The denominator of (B.1) reads At typical particle momentum in our system p ∼ O(T ), the Majorana neutrino decay width is of order Γ p O(λ 2 T ), and since we are working at temperatures T M , this results in ω p ∼ O(T ) as well. Hence in the four factors including Γ p , the latter is subdominant compared to the frequency terms. These factors are thus smallest when the two frequencies approach the mass shell. Alternatively, there is a pole due to the first factor. The function T thus has peaks in the following two limits, for which we estimate their parametric contribution: Sufficiently far away from the peaks we can simplify T by considering Γ p → 0, while keeping Γ p t = const., which is possible because of the smallness of Γ p [18], 23 )ω p sin(tω p ) (cos(tω 21 ) − cos(tω 23 )) + (ω 21 ω 23 − ω 2 p )(cos(tω p )(sin(tω 21 ) + sin(tω 23 )) − e Γpt/2 sin(t(ω 21 + ω 23 ))) . 3. ω p |ω 21 |, |ω 23 |: here we approximate (ω 2 For very small values of |ω 21 | we may expand the trigonometric functions in eq. (4.6), leading to where we have kept Γ p to compare with the peak region. For Γ p t ∼ O(1) one finds Comparing all estimates shows that the region around the peak at ω 21 = −ω 23 = ±ω p gives the largest contribution, T (t; ω 21 , ω 23 , p) ∼ O(λ −6 T −3 ). All other contributions are suppressed by a factor of order O(λ 4 ). Note that we later use λ 2 = 10 −8 , which implies O(λ 4 ) ∼ 10 −16 . This conclusion is corroborated by an independent consideration of the large time limit, which we give in appendix C. An exemplary plot of the function T (t; ω 21 , ω 23 , p) illustrates our findings numerically in figure 7.

Approximate integration of the lepton number
We are now going back to the full expression for the lepton number, eq. (3.36), and factor out the frequency integrals, Following the analysis in the previous section, we restrict the frequency integrations to a diagonal strip including the on-shell peak, ω 23 ∈ [−ω 21 − a, −ω 21 + a]. We choose a = const. × Γ p , such that ω p a Γ p (c.f. figure 7). Using the symmetry properties (3.25), we rewrite the ω 21 integration to obtain cos(ω 21 ∆t 21 + ω 23 ∆t 23 ) (4.13) Next, the σ h,ψ -functions can be checked numerically to vary at most at the percent level in the interval ω 21 , −ω 23 ∈ [ω p − a, ω p + a], as depicted in figure 8, with the conservative choice a = 100 Γ p . The variation is largest for small momenta p, which however do not contribute significantly to the integral. For the relevant momenta p ∼ O(T ), deviations from the on-shell values σ h,ψ (ω 21 = ω p ) are negligible. In a second step we may then only keep the leading term in a Taylor expansion around ω 21 = −ω 23 = ω p (c.f. eq. (3.25)), (4.14) allowing us to pull the σ h,ψ out of the integration. Note that it is sufficient to consider +ω p in this step, because we made use of the symmetry in ω 21 before. We can now do the integration in eq. (4.13) analytically leading to Using this on-shell approximated intermediate result, we are able to carry out the time integration of the lepton number density × cos(ω p ∆t 13 )e −Γp t 13 2 sin(a∆t 23 ) where Ei(x) denotes the exponential integral function and Si(x) the sine trigonometric integral function. Our choice a = 100Γ p satisfies a Γ p and allows us to approximately consider at → ∞, while keeping Γ p t fixed, leading to an a-independent result, p) . (4.17) Note that this expression only holds for t 1/Γ p , since for earlier times the approximation leading to eq. (4.6) is not valid, i.e. the region around the on-shell peaks is less suppressed so that the peaks are less pronounced. These memory and off-shell effects could of course be kept by a complete numerical integration without restriction to the regions around ω 21 = −ω 23 = ±ω p . Nevertheless, one makes the remarkable observation that at sufficiently late times our result has the same time dependence as the solution of the corresponding Boltzmann equations [18], but this time based on a full calculation without any assumptions.
Ultimately, we are interested in the value of the generated lepton number once it is thermalised. This is easily obtained by taking the limit t → ∞, resulting in As a valuable crosscheck for our approximate calculation, we change the order of integration in appendix C, i.e. we first do the time integrals, followed by the t → ∞ limit and only then integrate over the frequencies, arriving at the same result. Finally, it is also interesting to quantify the timescale for thermalisation, which is expected to be of order t T ∼ O(1/Γ p ) according to the discussion in section 2.1. We define t T as the time where which gives t T = 1/Γ for a constant Γ p ≡ Γ, cf. eq. (4.17).

Time evolution of the lepton number
The final integration over the Majorana neutrino momentum is performed numerically using integration and ODE algorithms implemented in the GNU scientific library [37] and the BOOST library [38]. Following [18], we consider a Majorana neutrino mass of M = 10 10 GeV and a constant temperature of T = 10 11 GeV for the time-dependent calculation. This choice is consistent with the physical scenario outlined in section 2.1, allowing to neglect the Hubble expansion. For the coupling of the corrected Majorana neutrino decay width we choose λ 2 = 10 −8 , in analogy to [20]. 6 The SM couplings are obtained by solving the renormalization group equations from [39,40] using the SM parameters from [3]. In analogy to the calculation of the asymptotic masses, we have neglected all contributions from the quark sector except for the top quark. All details are given in appendix D.
The resulting time evolution of the lepton number is shown in figure 9 (left) by the red curve. The blue curve corresponds to an alternative evaluation of the integrals, where we start from the time-integrated expression, T (t; ω 21 , ω 23 , p) in appendix B, and numerically integrate over the frequencies with ω 23 ∈ [−ω 21 − a, −ω 21 + a], this time with the larger range a = 1000 Γ p and the full functions σ h,ψ . As the figure shows, this crosscheck fully confirms the validity of our approximations at sufficiently large times. Deviations between In figure 9 (right) we compare our complete result to various other approaches discussed previously [18]. These are the solutions of Boltzmann equations (B), Boltzmann equations using quantum mechanical distribution functions (QB), and Kadanoff-Baym equations without gauge corrections, but with thermal widths for the lepton and Higgs propagators, γ = γ l + γ φ , introduced by hand (KB). 7 In all three cases we work with a constant decay width for the Majorana neutrino propagator, Γ p (ω p ) = Γ ∼ 10 −10 M ∼ O(GeV), corresponding to the dominant contribution in our full calculation. For the SM thermal damping widths we choose γ ∼ 0.1 T as expected from thermal field theory [28]. All results are normalised by the value of the thermalised lepton number density from the full calculation. While all curves have the same qualitative features, the Boltzmann result deviates by almost an order of magnitude, the quantum Boltzmann computation reduces the difference significantly and Kadanoff-Baym with thermal widths is O(1) accurate. We may thus conclude that gauge corrections indeed cause the expected narrow damping widths at late times. However, as expected in [18,19], for quantitative results the full calculation is mandatory.

Temperature dependence of the gauge corrections
In the simplest leptogenesis scenario considered here, the temperature is constant to first approximation. However, gauge corrections apply to other parameter choices and scenarios as well, and it is interesting to investigate their relative importance as a function of tem- perature. In figure 10 (left) we repeat the comparison of the previous section in case of the thermalised lepton number as a function of temperature T . Only including asymptotic masses leads to a kinematical suppression of the Majorana decay and production rate [20], and hence of the lepton number, in some temperature range, which is unphysical. This illustrates the importance of a complete calculation including the resummed gauge corrections. At yet higher temperatures, once T M , the kinematic suppression is switched off and only including asymptotic masses gives a good approximation of the full result. For low temperatures T < M , the results of Boltzmann, quantum Boltzmann and Kadanoff-Baym with thermal damping widths are very close. This agrees with the smallness of corrections observed in the non-relativistic case [23]. Nevertheless, the full computation of gauge corrections in this regime enhances the thermalised lepton number compared to the other solutions, for the chosen parameters Γ and γ. For temperatures T > M , the different solutions split up. Introducing thermal damping widths by hand repairs the qualitative deficit of kinematic suppression (and agrees well with the full result in a part of this region), but overestimates the lepton number for high temperatures. This is because the dominant momenta become more light-like, and hence the resummation more important. Note however that, by altering the scaling of the Majorana neutrino decay width to e.g. Γ ∼ 10 −11 T , better agreement with the full result is achieved.
Finally, figure 10 (right) shows the thermalisation times for the lepton number, extracted from the full result and the one including asymptotic masses only. Apart from the region with kinematic suppression we roughly find a scaling t T ∼ O(1/Γ p ) ∝ 1/T as expected.

Conclusions
In this work the leading-order gauge corrections to quantum leptogenesis were included in a fully quantum field theoretical calculation based on non-equilibrium Kadanoff-Baym equations. In this setting, gauge corrections appear first at the three-loop level of lepton self-energies. In the symmetric electroweak phase, gauge bosons are massless on tree-level and, in the kinematical region of collinear thermal loop momenta, require resummations to infinite loop order to complete a leading-order calculation in the gauge coupling. We have adapted a previous CTL resummation of the Majorana neutrino self-energy in such a way that it can be applied to the lepton self-energy graphs relevant for quantum leptogenesis. Our resulting lepton number density, eq. (3.36), is complete to leading order in all SM couplings and forms the basis for a quantitatively accurate evaluation of leptogenesis.
We have evaluated our expression for the simplest vanilla leptogenesis scenario at late times approaching thermalisation, which allows for an analytic analysis and integration of the lepton number density in Fourier space. One observes how the gauge corrections dynamically cause a dominant peak to grow in the integrand, whose width at late times corresponds to the thermal decay width of the Majorana neutrino. Similarly, thermal widths for leptons and Higgs propagators are automatically and fully included. The resulting late time behaviour qualitatively agrees with that observed from Boltzmann equations. However, for quantitative results, the full quantum calculation is necessary. Finally, we have evaluated the temperature dependence of the gauge corrections and the thermalisation time to estimate effects on other parameter choices or scenarios. Boltzmann equations dressed with quantum distribution functions underestimate the full result by O(1) for T M , introducing an uncertainty of similar size as that of spectator processes. Thus, we have also obtained a first complete estimate of the theoretical uncertainties for leptogenesis.
Our calculation can be easily adopted to other parameter regions or modfied scenarios, such as resonant leptogenesis. For yet better accuracy, the calculation could be further refined by also including washout diagrams as well as the temperature change due to the Hubble expansion.

D SM parameters for the numerical evaluation
In order to evolve the SM parameters to the scale given by the temperature µ = 2πT , we employ the renormalization group equations (RGEs) from [39,40] taking only into account the largest coupling to a quark, i.e. the top quark being the heaviest quark,