Source terms for electroweak baryogenesis in the vev-insertion approximation beyond leading order

In electroweak baryogenesis the baryon asymmetry of the universe is created during the electroweak phase transition. The quantum transport equations governing the dynamics of the plasma particles can be derived in the vev-insertion approximation, which treats the vev-dependent part of the particle masses as a perturbation. We calculate the next-to-leading order (NLO) contribution to the CP-violating source term and CP-conserving relaxation rate, corresponding to Feynman diagrams for the self-energies with four mass insertions. We consider both a pair of Weyl fermions and a pair of complex scalars, that scatter off the bubble wall. We find: (i) The NLO correction becomes large for $\mathcal O(1)$ couplings. If only the Standard Model (SM) Higgs obtains a vev during the phase transition, this implies the vev-insertion approximation breaks down for top quarks. (ii) The resonant enhancement of the source term and relaxation rate, that exists at leading order in the limit of degenerate thermal masses for the fermions/scalars, persists at NLO.


Introduction
In electroweak baryogenesis (EWB) the baryon asymmetry of the universe is created during the electroweak phase transition. The scenario requires new physics at the electroweak scale, in particular an extension of the standard model (SM) scalar sector to obtain a strong firstorder phase transition, and new sources of CP violation beyond those present in the CKM matrix. A major motivation to study EWB in detail is that it can be probed by experiment, for instance by collider searches for new scalars [1][2][3][4], precision Higgs studies [5,6], and CPodd collider observables [7][8][9][10][11]. Particularly tight constraints on new sources of CP violation come from measurements of the electric dipole moment [12][13][14][15][16][17]. In addition, gravitational waves produced during the phase transition may be measured by LISA or other future gravitational wave observatories [18,19].
To draw definite conclusions on the viability of EWB scenarios, accurate theoretical predictions are needed, which are currently somewhat lacking. Although no comprehensive comparison of different theoretical approaches exists, it has been shown for specific models that predictions may vary by more than an order of magnitude [20]. The computation is hampered by the nonequilibrium, non-perturbative, and finite-temperature aspects of the process. The starting point for a quantum treatment of EWB is the closed-time-path formalism at finite temperature. The dynamics of the phase space densities of plasma quanta are described by the Kadanoff-Baym equations [21][22][23]. Unfortunately, these equations are complicated, and a series of approximations -which are not always controlled -are needed to make progress.
We will focus on the vev-insertion approximation (VIA) method [24][25][26]. To derive a set of transport equations for the plasma particles that are simple enough to solve (numerically), the following main approximations are made: (i) The bubble dynamics is treated as slow compared to the typical time-scale of the plasma excitations. (ii) Quantum coherence effects are neglected, which allows to rewrite the Kadanoff-Baym equations in terms of number densities rather than phase-space densities. (iii) It is assumed that Fick's first law can be used to incorporate diffusion. And finally, (iv) the effective mass of the relevant plasma particles, which is spacetime-dependent in the bubble wall background, is treated as a perturbation.
The VIA method derives its name from the last approximation, as it corresponds to an expansion in the number of vev-insertions, that is, insertions of the two-point coupling, in the Feynman diagrams for the particles scattering in the bubble background. In this paper, we will determine the validity of this expansion, by calculating the next-to-leading order (NLO) correction to the CP-odd and -even rates for fermion/scalar interactions with the bubble wall. We will consider both a left-and right-handed fermion pair with a Yukawa-type interaction with the Higgs field, and a pair of complex scalars -also denoted by left and right -with a left-right-mixing coupling to the Higgs.
The transport equation for the right-handed field is of the form Similar equations exist for the left-handed field and for the Higgs boson. Here j µ R is the number current for the right-handed particles (the zeroth component corresponds to the number density of particles minus antiparticles), µ i the chemical potential of particle i, and c depends on the coupling to the Higgs field. The CP-violating term S cp sources the creation of a non-zero number density of right-handed particles, and it originates from interactions with the bubble wall. All other rates on the right-hand side of the equation conserve CP, and they give rise to washout of the number density. The relaxation rates Γ ± encode interactions with the bubble wall, Γ H interactions with the Higgs quanta, and the ellipses stand for all possible other interactions. A net chiral density in front of the bubble wall is transformed into a net baryon asymmetry by weak sphaleron transitions.
We will calculate the CP-even rates Γ ± and the CP-odd source term S cp to NLO in the vevinsertion expansion. They arise from the scatterings of the fermions/scalars with the bubble wall, mediated by the vev-dependent two-point coupling, f f (ϕ b ) and f s (ϕ b ) for fermions and scalars respectively, with ϕ b the bubble wall Higgs background. Let's focus on the fermions. Based on naive dimensional analysis one would expect the NLO correction to be suppressed by a factor |f f | 2 /T 2 N , with T N the temperature at bubble nucleation. However, it is not obvious that this estimate is accurate as there are other scales in the problem, such as the thermal decay widths Γ T of the particles. In fact, at leading-order in VIA the dominant rate and source term are both resonantly enhanced in the limit that the field-independent (thermal) masses for the left-and righthanded particles are degenerate, and then scale as S cp , Γ − ∝ 1/Γ T . Moreover, as we will show, since the correction can be order one for the top quark, it is important to check for possible numerical factors in the expansion parameter. It should be noted that the vev-insertion expansion is not a loop expansion and thus there are no factors of (4π) appearing.
The NLO expressions generically have a complicated form, but simplify enormously in the limit of degenerate masses for the left-and right-handed particles, which is a good approximation for the top quark. We then find that the ratio of the NLO to LO contribution for fermions is (Γ + vanishes in the mass degenerate limit): with m T the thermal mass, and α the QCD (electroweak) coupling for fermions with strong (electroweak) interactions. For scalars the enhancement factor also depends on the thermal width, and is given explicitly in eqs. (94) and (102). For quarks eq. (2) is of the order of the naive estimate |f f | 2 /T 2 N . The couplings f i are model dependent, and the validity of the vevinsertion approximation should be checked per model. In set-ups where only the SM Higgs is non-zero along the bubble wall |f f | ≈ y f ϕ b with y f the SM Yukawa coupling, and thus VIA breaks down for a pair of top quarks, but works well for all other SM fermions.
We further find that in the degenerate mass limit the NLO correction is just a multiplicative factor. The resonant condition of the leading order source term and relaxation rate is not shifted or otherwise affected, and is a robust result. This paper is organized as follows. In the next section, we briefly summarize the relevant formalism. We first introduce the fermionic and scalar set-up and define the two-point coupling f f and f s in section 2.1. Section 2.2 recaps the Feynman rules and propagator relations relevant for the calculation of the source terms in the CTP formalism, and section 2.3 gives the transport equations in the vev-insertion approximation. To calculate the source terms, multiple contour integrals need to be performed. To set the notation, we give the master integrals in section 2.4. In section 3 we review the leading order calculation for both the CP-violating and -conserving source terms. In section 4 the calculation is extended to the next-to-leading order contribution. We will provide all the calculational details. Readers only interested in the final results can go from section 2.3 straight to section 5 where we summarize the outcome of the NLO calculation, and discuss the implications.

Set-up, formalism and notation
In this section we introduce the set-up, and recap the formalism to calculate the source S cp and wash-out rates Γ ± in the transport equation (1) from Feynman diagrams. The leading order and next-to-leading order calculation is then presented in the next sections.

Model
We consider a two-flavor 1 system consisting of either a pair of chiral fermions (ψ L , ψ R ) or a pair of complex scalars (φ L , φ R ), and we are interested in their dynamics during the electroweak phase transition. In the high temperature expansion the thermal masses of the particles are field independent, and they are included in the free Lagrangian, and thus appear in the propagators defined in section 2.2. Scalars can in addition have a flavor diagonal constant mass term; this possibility can be included by substituting everywhere m 2 T → m 2 T + M 2 with M the bare mass. CP-violation resides in left-right mixing couplings to the Higgs field, which in the bubble wall background generates a field-dependent mass for the fermions/scalars. This mass is treated as a perturbation in the vev-insertion approximation (VIA), and included in the interaction Lagrangian. The Higgs vev can be parametrized as the space-time dependent bounce solution (in multi-Higgs models, we are interested in the linear combination of Higgs vevs that enters the bounce solution for tunneling). We can then write the two-point interaction for the fermions and scalars as with f i and i = f, s the CP-violating field-dependent mass term. It describes the scattering of the plasma particles off the bubble wall. The interaction violates CP, and particles and antiparticles scatter differently, provided Im(f * iḟ i ) = 0. In explicit models, this can be achieved if two (or more) different interactions with complex couplings contribute to f i , as then the relative phase cannot be rotated away in the bubble background. For example, CP-violation (CPV) can arise in a two-Higgs doublet model from interference between couplings to the two Higgs fields [28][29][30][31][32][33][34][35]. Alternatively, in an effective field theory approach it can come from interference between SM Yukawa and dimension-six effective interactions [13,16,36,37]. When we discuss the implications of our work in section 5 we will focus on the case where (ψ L , ψ R ) are the chiralities of a single SM fermion, e.g. the left-and right-handed top quark, but the results can be straightforwardly generalized to set-ups where the CP-violating interaction is between particles from different families [38,39]. For the scalar set-up the coupling f s can both be linear in the Higgs field (plus possible higher dimensional terms), as happens in supersymmetric models [26], and quadratic [40].
For our considerations, the origin of the mass term is not important, and we work with the phenomenological parametrization in eq. (3).
The evolution of the plasma quanta during the electroweak phase transition is described using the finite-temperature, non-equilibrium Closed Time Path (CTP) or Schwinger-Keldysh formalism [41][42][43]. For an extensive introduction to the CTP formalism, see e.g. [23].
All integrals and derivatives are performed along a closed path, that can be split into a plus-branch from the initial time (which for initial thermal equilibrium can be taken to minus infinity) to some finite time t, and a minus-branch going backwards. There are then four Green functions, depending on the branch that the time arguments of the fields are located on. The interactions connect fields on the same branch; the Feynman rule is that every mass insertion gives a factor −if j s for the scalar theory, with j = ± denoting the branch, and a factor −if j f / √ 2 for the fermionic theory. Since the minus branch runs backward in time the coupling picks up an additional sign Scalars For scalars the Wightman functions are defined as where we have suppressed the labels L, R on the fields and Wightman functions. Since the thermal masses are flavor diagonal, and the off-diagonal mass term is treated as a perturbation, all Green functions are diagonal in flavor space. The time-and anti-time-ordered propagators are with Θ(x) the Heaviside step function. Under complex conjugation the Green functions transform as where we introduced the notation G xy = G(x, y). Explicitly, the Wightman functions are given by with λ =>, < and k = d 4 k/(2π) 4 . We will also use the notation k = d 3 k/(2π) 3 and with We include the (leading order) thermal corrections in the propagators and take ω 2 k = k 2 +m 2 T (T ) with m T and Γ T the thermal mass and width respectively. We neglect deviations from thermal equilibrium for the thermal distribution functions, and expand in small chemical potential: with the Bose-Einstein distribution and its derivative given by Fermions The fermionic Wightman functions are defined via where we suppressed both spinor and flavor indices. The time-ordered propagators are defined by equivalent relations to the scalar ones eq. (6). The hermiticity properties are The explicit form of the Wightman function is with the spectral density again given by eq. (9), and m T the thermal mass. For future reference we introduce the notation that is, the subscript Tr(m) = 0 indicates that the mass term in the last factor of eq. (15) is set to zero. The thermal distribution functions for fermions are with the Fermi-Dirac distribution and its derivative given by

Transport equations
The transport equations can be derived directly from the Schwinger-Dyson equations [25,26], or equivalently, from the Kadanoff-Baym equations [44]. The latter approach makes transparant that only the leading order terms in the derivative expansion are included, i.e., that it is assumed that the bubble wall dynamics is slow compared to the typical timescale of the plasma excitations which is set by the temperature. The result for scalars is with Π λ xy ≡ Π λ (x, y) the self-energy. In the diffusion approximation j µ = (n, −D ∇n) with D the diffusion constant. Neglecting the bubble wall curvature, the transport equations reduce to ordinary differential equations for the number densities of plasma particles. The transport equation for a Dirac fermion current is with Σ λ xy the self-energy. Apart from the overall sign and the trace over spinor space, this has the same form as the scalar transport eq. (19).
The transport equations can be simplified using the hermiticity properties. The self-energies satisfy the same type of relations as the propagators in eqs. (7) and (14) ( We will verify this explicitly for the LO and NLO contributions in the next two sections.
We will always calculate the transport equation for the right-handed particle, the corresponding one for the left-handed particle then follows from ∂ µ j µ L,i = −∂ µ j µ R,i with i = f, s for fermions respectively scalars. Putting back the flavor indices, and using the hermiticity properties of the self-energy, this becomes with Θ xy = Θ(x 0 − y 0 ). In the expression on the right-hand-side we split the current in a CP-conserving source S cp i ≡ S cp R,i and a CP-violating source S cp i ≡ S cp R,i , and to avoid notational clutter we have dropped the flavor index. The relaxation rate extracted from the CP-conserving source term for the right-handed particles is 2 where we suppressed the scalar/fermion index. Equations (22) and (23) are the master equations to calculate the source terms and rates. We can expand the scalar self-energy in the coupling f s : The LO and NLO Feynman diagrams for the self-energy of the right-handed scalar are shown in fig. 1. Note that it is an expansion in the coupling f s , that is, in the number of vev insertions, not a loop expansion. The LO term is Π λ R,lo = O(f 2 s ), the NLO term Π λ R,nlo = O(f 4 s ), and the ellipses denote O(f 6 s ). An analogous expansion holds for the fermionic self-energies. Consequently, we can also expand the source terms as S I = S I lo + S I nlo + ... with I = CP, cp . 2 The rescaled relaxation rate that usually appears in the transport equations found in the literature is Γ ± M = 6 T 2 Γ ± . If the fermions/scalars are (s)quarks there is an additional Nc color factor.

Contour integral
To calculate the source terms we will have to perform multiple contour integrals. Here we introduce the master integrals, which also serve to set the notation. The integrals are of the form For the scalar model g(k 0 ) = g λ (k 0 ) the thermal distribution function eq. (11), for the fermionic model g(k 0 ) denotes the thermal distribution function eq. (17) times a k 0 -dependent trace factor. The function ξ I in the calculation of the CP-conserving and CP-violating source respectively is For I + the contour is closed in the upper half plane where the spectral function ρ(k 0 ) has two poles U (k). The contour for I − is closed in the lower half plane with poles at D(k), with Then where we introduced the notation Σ U = Σ k 0 =U (k) = Σ k 0 ={U 1 ,U 2 } and similarly for Σ D . Further F = 1 for the k 0 = ±E poles, and F = 0 for the k 0 = ±E * poles, as the former pick up an extra minus sign from the residue. If there are several contour integrals over k 0 i the notation is generalized to Σ DU = Σ k 0 1 =D Σ k 0 2 =U , etc. The sign is then denoted by (−1) F i . The final u 0 -integral converges as Im(k 0 ) cuts off the integral at u 0 → ∞ and The LO and NLO source terms require two and four contour integrals respectively. It will be useful to divide the summation over the poles into "opposite" halves, as these contributions either add or subtract. For example, consider the double sum Σ DU . If one term in the summation is (k 0 1 , k 0 2 ) = (D a , U b ) then we define the opposite term as (Dā, Ub) withā = 1 if a = 2 andā = 2 if a = 1. In this way we split the summation in halves, which we denote by where the 2nd sum on the right-hand-side contains all the opposites of the terms in the first sum (which of the pair of opposites is in DU is arbitrary). We can then use the relations U a = −U * a and similar for D to simplify relations.

Leading order source terms
In this section we review the LO calculation of the CP-conserving and -violating source terms appearing in the transport equation for the right-handed fermion/scalar.

Scalars
At leading order the self-energy of the right-handed scalar, given by the Feynman diagram in fig. 1, is where we used the notation G λ L,xy = G λ L (x, y), f x = f (x) etc, and λ ∈ {>, <}. We suppressed the label for scalars on f s . To get the final expression we used eq. (4). As f x f * y = (f y f * x ) * , and the Green's functions obey eq. (7), it follows that eq. (21) is satisfied as well. From eq. (22), the transport equation at LO is with Fermions The self-energy for the right-handed fermion at leading order is with P L,R the left-and right-handed projection operators. Substituting in the transport eq. (22), and inserting the propagators eq. (15), the trace is of the form Hence, we can neglect the mass-term in the trace and work with the propagators defined in eq. (16). Following the same steps as for the scalar case, the result for the CP and CPV source is

Derivative expansion
Let's start with the CP-conserving source term. At leading order in the derivative expansion, valid when the bubble wall background changes slowly compared to timescales set by the temperature, we can approximate Re[f x f * y ] ≈ |f (x)| 2 and take it out of the integral for the CP-even source in eqs. (33) and (36). Now insert the explicit form of the propagators eqs. (8) and (16). The integration over spatial coordinates can readily be done, and gives a delta function that sets all spatial momenta equal. We further introduce a new time coordinate u = x 0 − y 0 , such that the theta function becomes Θ(u) and we only have to integrate over positive time values. Then with The momentum k 0 1 corresponds to the left-handed particle with chemical potential µ L , the momentum k 0 2 to the right-handed particle with chemical potential µ R . Equation (37) is valid for scalars as well as fermions. For scalars s = 1 and the trace factor is trivial tr lo (k i ) = 1, f = f s , and the distribution functions g λ depend on the Bose-Einstein distributions eq. (11). For fermions s = −1, f = f f , and the distribution functions g λ depend on the Fermi-Dirac distributions eq. (17).
To calculate the CP-violating source term in eqs. (33) and (36) we expand the coupling where in the last step we only included the time-derivative term as the spatial part cancels in the source integral because of spatial isometry, and we defined δ = Im f (x)ḟ (x) * . The source term becomes which again applies to both scalars and fermions.

Contour/k 0 i Integrals
The integrals in the CP-conserving and violating source terms eqs. (37) and (41) and ξ cp = 1 and ξ cp (u) = u as in eq. (26). To do the contour integrals we close the contour for k 0 1 below and for k 0 2 above the real axis. There are two poles per integral, giving rise to four residue terms in total. Using eqs. (28) and (29) the result is with Now expand c +− = c +− 0 + c +− 1 + O(µ 2 i ) in small chemical potentials using eqs. (11) and (17), with the 0th order term c +− 0 independent of µ i and the the 1st order term c +− 1 linear in µ i . Using the Bose-Einstein and Fermi-Dirac distributions for scalars and fermions respectively, we get with s = 1 for scalars and s = −1 for fermions, and µ ± = µ R ± µ L . The summation in eq. (43) can be divided into the poles (a, b) ∈ DU and the opposite pairs (ā,b) ∈ DU , as discussed in section 2.4, where we introduced the notation that (a, b) denotes k 0 1 = D a , k 0 2 = U b . Then where we used g > (−k 0 , 0) = −g < (k 0 , 0) and h(−k 0 ) = h(k 0 ). This can be used to simplify J I lo .

Relaxation rate
For the CP-conserving source it follows that the leading order term in the µ-expansion cancels. The summation in eq. (43) is over At first order in the chemical potential the summation is over (a,b) , and we get instead Putting it all together, we find for the CP conserving source (setting k 0 1 = k 0 L and k 0 2 = k 0 R for the left-and right-handed particle respectively) from which we extract the relaxation rates (see eq. (23)) As the final step, we sum over the poles in DU which we choose (k 0 , where for the first combination (−1) F i = −1 and for the second (−1) F i = 1. The relaxation rate becomes with for scalars tr lo i = 1, and for fermions with tr lo given in eq. (39).

CP-violating source
Because of the non-trivial Ξ cp -factor in eq. (43) for the CP-violating source, the leading order term in the µ-expansion already contributes. Indeed (iA 0 Ξ cp )| (ā,b) = (iA 0 Ξ cp ) * | (a,b) and Putting it all together, then Choosing to sum over (k 0 with tr lo i defined in eq. (52). In the second term between square brackets, the term ∝ s in the numerator diverges, also in the zero-temperature limit; it can be removed by adding a counterterm and is absent in the renormalized theory [45].

Next-to-leading order source terms
In this section we calculate the source terms in the transport equation eq. (22) at NLO. The self-energy at this order is given by the second Feynman diagram in fig. 1. It is still a tree level diagram, but now with four mass insertions. The NLO diagram thus represents the next order contribution in the coupling expansion, and is O(|f | 2 ) suppressed with respect to the LO diagram.
Although the extra mass insertions complicate the calculation of the relevant integrals, since the self-energy diagram remains tree level, all steps in the calculation are essentially the same as for the LO calculation.
Scalars Let's start with the scalar set-up. The self-energy of the right-handed scalar at NLO is where we summed over all possible ±-signs for the internal vertices. Using eq. (7) it follows that (Π −+ R,nlo (x, y)) † = Π −+ R,nlo (y, x). We can then extract the CP-conserving and -violating source terms from eq. (22) Fermions The self-energy for the right-handed fermion at NLO is which has the property that Σ −+ R,nlo (x, y) † = γ 0 Σ −+ R,nlo (y, x)γ 0 . When used in the transport equation eq. (20), the trace of the ΣS and SΣ-terms is respectively The last expression holds since the spatial integration gives a set of delta functions that set k i = k j . We can then use that Tr γ 5 / k 1 / k 2 / k 3 / k 4 = −2 µρκσ k µ 1 k ρ 2 k κ 3 k σ 4 = 0. The source terms for fermions are then analogous to those for the scalars withÎ 1 ,Î 2 as in eq. (58) with the replacement G xy → S xy .

Derivative expansion
Let's start again with the CP-conserving source. To lowest order in the derivative expansion f x f * z 1 f z 2 f * y ≈ |f (x)| 4 , and it can be taken out of the integrals for the CP even source in eqs. (57) and (61). Now we express all time and anti-time ordered propagators in terms of the Wightman functions using eq. (6), and insert the explicit expressions for the latter eqs. (8) and (16). The integration over spatial momenta and coordinates sets all spatial momenta equal. The CP-even source becomes The particles with momenta k 1 , k 3 are left-handed, those with momenta k 2 , k 4 right-handed.
Here c i 1 i 2 i 3 i 4 denotes the combination of distribution functions with i j = ± andī j is + if i j is −, and vice versa. g For scalars s = 1 and tr nlo = 1, while for fermions s = −1 and the trace factor now becomes with k i · k j = k 0 i k 0 j − k 2 . For the CP-violating source we expand lim y,z 1 ,z 2 →x Integrating over spatial momenta, the CPV source terms in eqs. (57) and (61) are and tr nlo , θ i 1 i 2 i 3 i 4 , c i 1 i 2 i 3 i 4 the same as for the CP-even source given in eqs. (63) to (65).

Contour/k 0 i Integrals
The integrals in the CP-conserving and -violating source terms eqs. (62) and (67) can be written as S cp nlo = 2s|f | 4 Re k J cp nlo , and S cp nlo = 2s|f | 2 δIm k J cp nlo , with We will work out two example terms J I For the CP even integral ξ cp = Ξ cp DU U U = 1, for the CP odd integral Ξ cp DU U U is given in eq. (76) below. In the last line of eq. (69) we did the u, v, w integration using eq. (29). For the summation over poles we used the notation defined in section 2.4.
To evaluate all terms in eq. (68) it is sometimes necessary to split the integration interval using theta-functions, and to do further coordinate transformations. Consider J I ++−− , which consists of two terms, see eq. (64). To do the contour integral for the first we use 1 = Θ(u − v + w) + Θ(−u + v − w) to split the k 0 1 integral; the result is J I In the first term we can make the and it follows that the first and third term in eq. (71) cancel. To calculate the remaining term, with Ξ I DDU U given in eq. (76) below. In this way we can evaluate all terms in J I nlo . The result is with F defined in eq. (70) and For the CP conserving source Ξ cp i = 1 are trivial, while for the CP violating source .
(76) The c-factors in eq. (74) can be rewritten in terms of the distribution functions where we used that (g > i − g < i ) = 1. The c +− -terms defined in eq. (38) are expanded in small chemical potential as in eq. (45).

Relaxation rate
For the CP conserving source, the leading order term in the small µ expansion cancels, just as for the LO calculation. To see this, divide the summation in (a 1 , a 2 , a 3 , a 4 ) ∈ DU U U and the opposite half (ā 1 ,ā 2 ,ā 3 ,ā 4 ) ∈ DU U U , as discussed in section 2.4. For zero chemical potential the first term in eq. (74) is proportional to where we used that for zero chemical potential g > (−k 0 ) = −g < (k 0 ). Similarly, the other terms in J cp nlo vanish at leading order. Hence, the first contribution comes from the term linear in the chemical potential. We can again divide the summation into two halves, and using that h(−k 0 ) = h(k 0 ) now they add: Putting it all together we find that with µ 1 = µ 3 = µ L and µ 2 = µ 4 = µ R . We first do the sum over the k 0 i appearing in the h-functions. For the first term, we choose the set DU U U = (1, 1, a 3 , a 4 ), (2, 1, a 3 , a 4 ) = (E * L , E R , a 3 , a 4 ), (−E L , E R , a 3 , a 4 ) with a 3 , a 4 unconstrained. The remaining sum is over the poles in k 0 3 = k 0 L and k 0 4 = k 0 R . Similarly, for the second term we choose DDDU = (a 1 , a 2 , E * L , E R ), (a 1 , a 2 , −E L , E R ), and for the third term DDU U = (a 1 , −E R , −E * L , a 4 ), (a 1 , −E R , E L , a 4 ). We write where we factored out the leading order trace factors tr lo i defined in eq. (52). The coefficients ∆ cp i are with tr nlo i = 1 for scalars, and for fermions with tr nlo given in eq. (65). Finally we do the summation over the remaining poles in k 0 L,R . We find ∆ cp ≡ ∆ cp 1 = ∆ cp 2 for both scalars and fermions with For scalarstr nlo = 1 while for fermions

CP violating source
Since there is an extra factor Ξ cp in J cp nlo in eq. (74), the leading order term in the small chemical potential expansion contributes for the CP-violating source. Just as in the LO calculation we can effectively set µ i = 0. The summation over poles can again be divided into opposite halves, which add. The result is then We sum over k 0 i appearing in the number densities n; we make he same choices as for the CP-conserving source in the previous subsection. We can then write We can do the final summation over poles in ∆ cp i , but unlike the CP conserving case where some terms cancel and the final result is reasonably simple, the result cannot be written in a concise way. What is more, we now find ∆ cp 1 = ∆ cp 2 and also ∆ cp i is generically complex. The NLO contribution can thus no longer be written as a simple additive factor in the integrand.

Results and discussion
In this section we summarize the main results for the NLO calculation, and discuss the implications.

Relaxation rate
The full relaxation rate up to NLO is with the NLO contribution captured by ∆ cp given in eqs. (84) and (85). The rate Γ − is dominated by the first term in the square brackets, which has a resonance in small thermal widths for degenerate masses. In this limit the expression for the rate and the NLO correction simplifies considerably, and we can derive analytical expressions. To see the resonance, we write m 2 L = m 2 R + δm 2 and expand in small mass difference |δm 2 | m 2 L . The first term between square brackets in the expression for Γ − in eq. (89) becomes where we additionally approximated Γ R ≈ Γ L , and in the denominator assumed Γ L T as is appropriate for a perturbatively generated thermal width. If in addition |δm 2 | 4Γ L T , the mass difference in the denominator can be neglected. This is an excellent approximation for the quarks. On the other hand, for leptons the δm 2 -term in the denominator in eq. (90) dominates. Since δm 2 T 2 this term is still enhanced compared to the 2nd term between square brackets in eq. (89), but the parametric dependence of the resonance is different. The second term between brackets in Γ − and both terms in Γ + do not have this resonant structure and are subdominant. Indeed, expanding in small mass difference gives It thus follows that the rate Γ − is resonantly enhanced in the degenerate mass limit. The NLO contribution ∆ cp is an additive factor to the integrand, and thus also to the resonant factor in the integrand; it does not affect or shift the resonance structure itself up to O( Γ L ω L , δm 2 To estimate the NLO contribution to Γ − we will for simplicity consider the exactly degenerate case, a good approximation for quarks, and write with thermal width Γ 2 T T 2 . Then The traces are trivial for scalars, and tr lo 1 = 2tr nlo = (m 2 T + Γ 2 T ) for fermions. Putting it all together we get for scalars and fermions respectively where the ellipses denote O( Γ L ω L , δm 2 ω 2 L ) corrections. When the NLO contribution dominates and |∆ cp i | > 1, the relaxation rate becomes negative, an unphysical result, signalling the breakdown of the vev-insertion approximation. This is expected to be cured when the full tower of the higher order terms in the coupling expansion is added, which amounts to a resummation of the spacetime-dependent mass.
Whether the VIA expansion is valid depends on the details of the model. To get some insight we estimate the thermal mass and width as m 2 T ∼ αT 2 , Γ T ∼ αT , with α = g 2 4π with g the largest relevant (gauge) coupling. The integral in eq. (94) is dominated by momenta ω T ∼ T .
Let's first look at scalars. In supersymmetric set-ups, the interaction with the Higgs is typically trilinear and we parametrize |f s | = 1 √ 2 A s ϕ b . The relaxation rate is maximized in the broken phase with ϕ b = ϕ N , the Higgs vev at nucleation. For a strong first order phase transition we need ϕ N /T N 1, with T N the nucleation temperature. During the phase transition we then find For neutralinos the thermal corrections are set by the weak gauge interactions and α ≈ 0.03. The NLO contribution is important unless the trilinear coupling is significantly smaller than the electroweak scale A s < O(0.1)T N . We conclude that VIA breaks down in these type of SUSY models.
For fermions with a Yukawa type Higgs interaction we parametrize |f f | ≈ỹ f ϕ b . Possible corrections from higher order Higgs interactions, although maybe essential for CP-violation, can be neglected in |f f |. If only the SM Higgs field obtains a vev during the phase transition theñ y f = y f equals the SM Yukawa coupling. We estimate shows R cp for quarks (leptons), and the dashed orange and green lines show R cp for quarks and leptons respectively. The black line shows the estimate of R cp and R cp for quarks in the mass degenerate limit, corresponding to eq. (94) and eq. (102) with ∆ cp , ∆ cp estimated by eq. (96). We take ϕ N = 100 GeV and T = 100 GeV .
In fig. 2 we show the ratio of the NLO and LO contribution to the relaxation rate R cp ≡ |Γ − nlo |/|Γ − lo | deep into the broken phase (where ϕ b → ϕ N ) as a function of the Yukawa couplings for quarks and leptons in red and blue respectively. The estimate of eq. (94), with eq. (96) for the quarks in the mass degenerate limit is shown in black and matches the full result up to a factor O(1.5). The good agreement between the numerical and analytical calculations confirms that the relaxation rate is dominated by the resonance. We indeed find that the (absolute value of the) NLO result for quarks becomes comparable to the LO result forỹ f 1.
For leptons the ratio R cp is approximately a factor 8 larger than for quarks for the same coupling, and as a result the VIA breaks down already at smaller couplingỹ f = O(0.1). Nonperturbativity is never an issue for leptons if the bounce is along the SM direction, but may still be in other set-ups that haveỹ f y f . As discussed above, for leptons the exactly degenerate limit eq. (92) and thus the estimate eq. (94) is not a a good approximation, but we can nevertheless understand why R cp is larger. For quarks the NLO correction eq. (99) is enhanced by α −1 , with α the QCD coupling constant. For leptons we similarly expect that the NLO correction scales with inverse powers of α, now arising from inverse powers of δm 2 rather than the thermal width. For leptons, though, the thermal corrections are set by the weak interactions and α is the weak coupling constant, giving a larger enhancement.

CP violating source
Adding the LO and NLO contribution to the CP-violating source term gives where we removed the divergent term by normal ordering, and with the coefficients ∆ cp i given in eqs. (83) and (88). Just as for the relaxation rate, the first term is resonantly enhanced in the degenerate mass limit Unlike the relaxation rate, ∆ cp 1 = ∆ cp 2 and Im[∆ cp i ] = 0, and the NLO contribution is in general not simply an additive part in the integrand. However, the resonant term only depends on ∆ cp 1 which is a real factor in the mass degenerate limit eq. (92). In particular, we then find where the ellipses denote O( Γ L ω L , δm 2 ω 2 L ) corrections. Apart from the overall sign the corrections are exactly the same as for the relaxation rate eq. (94). For fermions the VIA approximation thus also breaks down for the calculation of the CP-odd source for couplingsỹ = O(1) for quarks. Figures 2 and 3 show the numerical results for the source. In fig. 2 we plot the ratio R cp ≡ |S cp NLO |/|S cp LO | for quarks (dashed orange) and leptons (dashed green) as a function of the coupling y f . Not surprisingly giving our estimate in the degenerate limit, the ratio is nearly identical to R CP for the relaxation rates. Although the source terms is space-time dependent and varies over the bubble wall, this cancels out in the ratio R cp which is evaluated deep in the broken phase. In fig. 3 we plot the LO and NLO contribution across the bubble wall for a quark for two different choices of the Yukawa coupling. For this plot we took the bubble profile ϕ b (z) = ϕ N 2 (1 + tanh z/L w ). The CP-violation stems from a complex dimension-6 Yukawa-like interaction and f f =ỹ f ϕ b (1+iϕ 2 /Λ 2 ) with Λ the cutoff scale. Evaluating the ∆ cp with ϕ b = ϕ N gives a good estimate for the size of the NLO correction in the z-region where the source peaks.

Discussion
To summarize, we have shown that the NLO contribution to the relaxation rates and source terms is small as long as |f s |/T 2 N α for scalars and |f f |/T N √ α for fermions, with α the QCD (electroweak) coupling for scalars/fermions with strong (only electroweak) interactions. For larger effective couplings |f i | the vev-insertion approximation breaks down.
Focussing on specific implementations, in a supersymmetric setup with CP violation in the neutralino sector, this implies that VIA breaks down for trilinear couplings |f s | ≈ A s ϕ b for : z-dependence of the LO (solid lines) and NLO (dashed lines) source terms for a quark withỹ f = 0.1 (red) andỹ f = 1 (blue). We use the bubble wall profile ϕ b (z) = ϕ N /2(1 + tanh z/L w ), with L w = 10/T ϕ N = T N = 100 GeV and |f f | ≈ỹ f ϕ b . CP-violation stems from a dimension-6 operator as in Refs. [13,16,36,37], giving δ = 3ỹ 2 f v w ϕ 3 b ϕ b /Λ 2 , where the prime denotes a derivative with respect to z, we take Λ = 1 TeV and bubble wall velocity v w = 0.05. forỹ f > O(1) for quarks, andỹ f > O(0.1) for leptons. This means in particular that if only the SM Higgs field obtains a vev during the phase transition, the source terms cannot be reliably computed for the top quark, but a baryogenesis scenario with CP violation in the lepton sector [38,39,[46][47][48][49] is under calculational control.
We finally note that the resonant behavior of the relaxation rate and source term is preserved at next-to-leading order.