Analytic Techniques for Solving the Transport Equations in Electroweak Baryogenesis

We develop an efficient method for solving transport equations, particularly in the context of electroweak baryogenesis. It provides fully-analytical results under mild approximations and can also test semi-analytical results, which are applicable in more general cases. Key elements of our method include the reduction of the second-order differential equations to first order, representing the set of coupled equations as a block matrix of the particle densities and their derivatives, identification of zero modes, and block decomposition of the matrix. We apply our method to calculate the baryon asymmetry of the Universe (BAU) in a Standard Model effective field theory framework of complex Yukawa couplings to determine the sensitivity of the resulting BAU to modifications of various model parameters and rates, and to estimate the effect of the commonly-used thin-wall approximation.


I. INTRODUCTION
A long-standing challenge of particle cosmology is to understand the mechanism by which the baryon asymmetry of the Universe (BAU) is generated. The Standard Model (SM) prediction [1,2] is many orders of magnitude smaller than the observed value of Y obs B ≈ 8.6 × 10 −11 measured by PLANCK [3]. The requirement for the dynamical process that generates the asymmetry to occur out of thermal equilibrium implies a particular structure for the particle dynamics. In electroweak baryogenesis (EWBG, for reviews, see e.g. Refs. [4][5][6]), one calculates the asymmetry that is produced during the electroweak phase transition, as bubbles of non-vanishing vacuum expectation value (VEV) of the Higgs field form and expand to fill the Universe [4,5,[7][8][9][10]. The important dynamics in such a scenario arise from the CP-violating interactions, which occur across the bubble walls and lead to a chiral asymmetry. Weak sphalerons then convert this chiral asymmetry into a baryon asymmetry by acting only on left-handed fermions and changing the baryon number. The importance of diffusion and the role of leptons was identified in Refs. [8,[11][12][13][14][15]. Since the strong sphalerons only wash out the quark asymmetries, and the diffusion into the symmetric phase is larger for leptons, the τ as a lepton with a sizable Yukawa coupling becomes an efficient source for CP violation [12,14,16].
Typically, the calculation is performed semi-classically, such that the particle dynamics is encoded in transport equations -a set of coupled, linear, non-homogeneous differential equations. The solution to these equations determines the eventual densities of each particle species, yielding a prediction for the baryon asymmetry. The current state-of-the-art approaches for solving these transport equations are the following: • Making a set of approximations that simplify the transport equations into a single equation that is analytically solvable and qualitatively understandable [12,14]; • Solving the full set numerically, which is more accurate but makes it difficult to gain physical insight into the solution [14]; • Solving the full set semi-analytically through a process of redefinitions that allow singling out equations to be solved individually as a recursive process [17].
We propose a new, semi-analytic method, which is similar to the latter approach, but simpler in several respects. Its implementation and usage are clear, and the understanding of algebraic features provide an intuitive picture of the physical process. Moreover, under mild approximations, this method allows for a fully-analytic solution, which is useful for estimating the accuracy of the corresponding semi-analytic calculation. Because the approximations are mild, a good agreement between the semi-analytic and the exact solution suggests that the semi-analytic results are reliable also in the original form of the equations and can be extended to more general scenarios.
The paper is organized as follows. In Sections II and III we solve a general set of transport equations, and impose the suitable boundary conditions. In Section IV we discuss the importance of zero modes and illustrate a way to treat them in a numerically stable way.
Section V describes techniques that can be applied to produce a fully-analytic solution in applicable cases. In Section VI we apply our method to calculate the baryon asymmetry in several scenarios within the SM effective field theory (SMEFT) framework of complex dimension-six Yukawa terms, testing the sensitivity of the produced asymmetry to modifications by factors of O (10) to model parameters such as the bubble wall parameters and the rates that are an input to the transport equations and have sizable uncertainties. We summarize and discuss our results in Section VII. The Appendices A -C provide details of derivations, definitions and benchmark parameters, as well as several consistency checks.

II. CONSTRUCTION AND GENERAL SOLUTION
In the following two sections, we will work in what is known as a two-step approach [11,18], where the particle dynamics are approximated as a two-step process: In the first step, CP-violating interactions generate a chiral asymmetry, and in the second step, the weak sphaleron process acts on the chiral density and converts it into a baryon density 1 . This decoupling is possible because the weak sphaleron rate is typically slow compared to other processes (see App. A 3). In App. B 4 we show a comparison between the two-step approach and the one-step approach, where the weak sphaleron is incorporated to the transport equations directly.
The second step consists of solving a single differential equation for the baryon density, and is described in detail in Appendix C. Solving the transport equations of the first step generalizes the solution of a single equation to a set of equations, one for each particle, and is the focus of this paper.
Taking the diffusion approximation [7,8] for the particle density f with the notation where v w is the wall velocity and D f the diffusion coefficient, a typical set in the two-step approach is the following (2.1) The CP-violating sources S i , the k i -functions and the rates Γ i are calculated by standard methods [10,[19][20][21] and their values in our framework appear in App. A. The chemical potentials are related to number densities via If we absorb the factor T 2 /6 in the definition of the effective chemical potentials for each process, their values are given by [14,22] The up quark is a representative of the other light quarks (d, s and c): since they interact only via the strong sphaleron to a good approximation, they are linearly dependent and hence redundant [14].
The sources peak in the broken phase, and for simplicity we approximate the bubble wall as a step function at z = 0, the center of the bubble wall (see Sec. VI C for further discussion on this choice). We consider the rates to be constant at each phase (possibly with different values), while for the sources we maintain their z-dependence in the broken phase, and eliminate them in the symmetric phase. We thus obtain a set of linear equations with constant coefficients for each phase. With N denoting the number of species appearing in the transport equations (for the set in Eq. (2.1), N = 7), we replace these N equations of second order with 2N equations of first order by defining Here Γ is a matrix of couplings between different particles, where each entry is of the form The general solution to the homogeneous part for each species is a linear combination The weights A f i j are determined, up to an overall normalization factor, by the eigenvectors of K. We can thus writeχ in vector form as follows:χ where (f i ,ḡ i ) T are the eigenvectors of K, and C i are integration constants. We organize the eigenfunctions in a z-dependent matrixΦ(z). Using variation of parameters, the full solution in the broken phase is We provide the numerical agreement between g i and f i of the solution in App. B 2. The impact of including more particles species in the set of transport equations is investigated in App. B 3. Furthermore, in App. B 1 we show the conservation of B − L numerically.

III. BOUNDARY CONDITIONS
In each phase, half of the modes decay and the others diverge or are constant. We choose boundary conditions as follows: • In the symmetric phase (z < 0), the integration constants of both the divergent and zero modes are set to 0, complying with the assumption that no baryon asymmetry is present before the electroweak phase transition.
• In the broken phase (z > 0), the integration constants of divergent modes are used to cancel the divergent integrals coming from the non-homogeneous terms in pairs.
• The remaining modes are determined by imposing continuity ofχ at z = 0. Sinceχ contains the vector of derivativesḡ, this is equivalent to requiring continuity of each particle density and its derivative at z = 0.
An important observation is that all modes either decay or are chosen to vanish at infinity, except for the zero modes. These are the only ones to survive deep in the broken phase z → ∞. Therefore, their existence is crucial for the success of EWBG (and is indeed guaranteed by the linear dependencies in Eq. (2.1)).
The solution of Eq. (2.5) in the broken phase for the i'th component ofχ is whereλ is a diagonal matrix constructed from the eigenvalues λ j andφ ij is a matrix of the corresponding eigenvectors. We denote integration constants of positive (negative) eigenvalues by +(−), and a B(S) superscript indicates the broken (symmetric) phase. For positive eigenvalues in the broken phase we choose This choice guarantees convergence at infinity. The continuity conditions are treated as follows. In the broken phase at the phase boundary, Eq. (3.1) reads whereb is a constant vector with entries b i obtained from Eq. (3.2).
In the symmetric phase, we set the integration constants associated with negative eigenvalues to zero, such thatχ Continuity at z = 0 is then To reach the final expressions, we need to solve a linear set of equations for the remaining integration constants. We can organize these constants in a vectorc ≡ (C S+ , C B− ) T and the corresponding modes as columns of a matrixÂ ≡ φ S | − φ B , such that finding the remaining integration constantsc amounts to solving the equationÂc =b. We can then collect the relevant densities, which in the two-step approach involves summing over the densities of the left-handed multiplets in the symmetric phase. In the case of Eq. (2.1), we recall that u acts as a representative of the light quarks. To obtain the densities of the left-handed multiplets of the first two generations, we relate them to u via q 1 = q 2 = −2u [14]. The chiral density is n L = q + l − 4u, which we plug into Eq. (C10) to solve for the baryon asymmetry. Note that since we only need the zero modes for our final result of Y B , Eqs. (3.2) of λ j = 0, (3.5) and (C10) imply that Y B is exactly linear in the integrated CP-violating sources S f .

A. One step and two step approaches
To obtain the baryon asymmetry in the one-step approach, we need to add to Eq. (2.1) the following terms: In this case, the degeneracy among light quarks in Eq. (2.1) is explicitly broken. Therefore we must reintroduce at least one left-handed quark multiplet. We may keep one quark generation implicit as long as we add its contribution to Y B in the end. The baryon density is obtained by summing over the zero modes of each species in the broken phase, and multiplying the quark densities by 1/3. The convergence of the two-step approach towards the one-step solution for small Γ ws is shown in App. B 4.

IV. ZERO MODES AND NUMERICAL REGULARIZATION
We have seen that zero modes are crucial for the generation of a baryon asymmetry, The block Γ corresponds to couplings in the transport equations, which we know are not all linearly independent: In the two-step approach, the couplings of left-handed multiplets are the negatives of the corresponding right-handed ones (e.g. ∂l = −∂τ ). Thus each generation produces a zero mode. In the one-step scenario, the degeneracy is broken between left and right, but reintroduced across species. For example, When incorporating many particle species in the transport equations, finding the eigenvalues of K is an intrinsically numerical task, equivalent to finding roots of high-order polynomials. The zero modes, which necessarily exist, may cause numerical instabilities if not treated carefully. A way to circumvent the problem is to first perform a partial diagonalization of K to extract the zero eigenvalues, and then solve for the rest of the system independently. Let us outline the procedure. Suppose we have a matrix M for which we know only a subset j of its m eigenvalues. We would like to find a matrix U such that whereD is diagonal and consists of the j known eigenvalues of M , andM is arbitrary. If M is diagonalizable, then in particular it is partially-diagonalizable. In the case of Eq. (2.1), of the right-eigenvectors that were already found, and V r remains to be determined. We can where U l are the left-eigenvectors and V l the remainder. We have From this, we see that V r and U l span orthogonal spaces, such that The upshot in our case is that we found a way to reduce the original problem of finding the eigenvalues of the singular matrix K to finding the eigenvalues of a regular matrixK, which should be numerically stable. Going back to the general case, we now need to match the eigensystem of the transformed matrix M to that of the original matrix M . The eigenvalues are the same, as can be seen from For the eigenvectors, suppose y is an eigenvector of the transformed matrix, and denote x = U y. Then, We find that if y is an eigenvector of M corresponding to an eigenvalue λ, then x is an eigenvector of M , corresponding to the same eigenvalue λ.
To summarize the procedure, we start by finding the eigenvectors of the zeros of K to obtain U , which we use to partially diagonalize K. We then find the eigensystem of K , and transform the eigenvectors to obtain the eigensystem of the original matrix.

V. BLOCK DECOMPOSITION FOR ANALYTICAL SOLUTION
In this section we show that, under certain approximations, we can obtain a fully analytic solution. This is useful for checking the semi-analytic method, where the eigenvalue problem is solved numerically, and consequently all downstream calculations are numeric as well.
Since the approximations we are going to use are mild, finding that the results are in good agreement means we should expect the semi-analytic method to be reliable also for the exact equations.
Using the general structure of K, we get We obtained equations that are, first, independent of ϕ i , and second, close to representing an eigenvalue problem for an N × N matrix instead of 2N × 2N . If we assume the diffusion coefficients are all the same, then V becomes a scalar matrix, and we obtain an actual eigenvalue problem for the matrix Γ, with eigenvalues given by This can also be seen from determinant properties of block matrices. For a general matrix, In our case, such that if all the diffusion coefficients are the same (i.e. ∀i, j : Neglecting the Higgs density h ≈ 0 and decoupling the weak sphaleron (two-step approach) allows us to solve the eigenvalue problem (5.4) fully analytically. Doing so allows us to obtain the eigenvectors of K by solving the quadratic equations (5.3) for λ ± i and using the relation (5.2) to constructΦ i .
We can obtain an analytic solution also without assuming the diffusion coefficients are all equal, and instead approximate them as equal only among fields from the same family: This approximation allows us to arrange Γ in blocks of equal D's for quarks and leptons separately. Then, each block is a subproblem of the original eigenvalue problem, which will be solved separately. If we look at Eq. (2.1) under the above approximations, we have q, t, b, u, l, τ , which naively form a 12 × 12 matrix, but reduces to separate 4 × 4 and 2 × 2 blocks, which are easily solvable. Of course, this decomposition works also under the more aggressive approximation of equal diffusion coefficients. If we do not make the approximation D l ≈ D τ , then the quark block still forms an eigenvalue problem with an effective eigenvaluẽ λ q ≡ λ 2 − λV q , but the lepton block does not. Instead, it is just a set of two equations in 3 variables: λ, φ l , φ τ , where we denote for simplicity the latter two by l, τ , respectively. The equations are thus Setting τ = 0 immediately implies l = 0, trivializing the solution. We can therefore choose τ = 1, which gives Plugging in the values of Γ ij eliminates the constant term, reproducing the expected zero eigenvalue, and leaving us with Note that the zero eigenvalue determines the eigenvector to be The other eigenvalues are given numerically by where λ = 0.044 leads again to the eigenvector 13) and the other two eigenvalues both produce the eigenvector (l, τ ) = (−1, 1) . (5.14) The apparent degeneracy in the eigenvectors is resolved when we construct the full eigenvectors usingφ i = λ iφ i , where all the quark entries are 0.
The next step in the process is to invert the matrixΦ corresponding to the eigenvectors of K. Instead of directly invertingΦ, which is computationally taxing, we will follow a similar path to the regularization procedure, and find the eigenvectors of K T as we did for K. These will be the left eigenvectors of K, and when properly normalized construct the inverse ofΦ. Denote the left eigenvectors by (ū,v), such that The effective eigenvalues are again the same, and the eigenvectorsv i are solved for and used to obtainū i . We organize the left eigenvectors as rows in a matrix A, and choose their normalization such that A =Φ −1 . From here on we simply follow with the semi-analytic procedure, and eventually plug in the numbers for the baryon asymmetry with arbitrary precision. The agreement between the fully and the semi-analytical solution is numerically investigated for two different sets of assumptions in App. B 5.

VI. PARAMETER DEPENDENCE
In this section we discuss how various model parameters and rates affect the baryon asymmetry. Our detailed calculations are performed in the framework of a Standard Model Effective Field Theory (SMEFT) with dimension-six complex Yukawa terms. This framework thus introduces new sources of CP violation, but does not enhance the electroweak phase transition which is assumed to be addressed separately. The Lagrangian for dimension 4 and dimension 6 Yukawa-type terms is given by: where v 0 is the Higgs VEV at zero temperature. The definitions of relevant quantities and the benchmark values for the numerical calculations are given in Appendix A. In particular, the benchmark values for T f I appear in Table II. In the examples shown here, we set T f R = 0 for all species. The phenomenology of the muon and third-generation fermions, including the interplay of T f R and T f I , is analyzed in detail in Refs. [16,23]. For T f R = 0, see also Refs. [14,21].

A. Relaxation and Yukawa rates
Consider the relaxation rates Γ M , Γ Y that appear in Eq.(2.1) and explicitly defined in Eq. (A13). These are CP-conserving terms that for large values tend to produce chemical equilibrium and dampen the asymmetry. They are calculated to leading order in perturbation theory. Higher-order corrections and terms beyond the underlying approximations are expected to modify these rates, see e.g. Ref. [20,24], and consequently have an impact on the calculated baryon asymmetry. Here we do not include these higher-order terms. Instead, we study the sensitivity of the baryon asymmetry to modifications of Γ M and Γ Y . In Figure 1 we replace and plot Y B as a function of the modifiers κ f M/Y , allowing for large deviations from the leading-order value.  Figure 1: Y B as a function of the modifier κ f M/Y of the relaxation and Yukawa rates shown for the τ and t sources.
For the tau, changes of O(10) to the rates translate to only O(1) changes in Y B . The top is much more sensitive to changes in the relaxation rate due to its large mass: an On the other hand, the larger Γ Y , the larger Y B . This may be an effect of avoiding the washout due to Γ M by transferring some density to other species with slower rates. To illustrate this point further, we integrate the number densities of each particle species in the symmetric phase, prior to the weak sphaleron action. We denote the integrated density of particle f in the symmetric phase by N f = 0 −∞ dz n f (z). In Fig. 2, we show for each source how the integrated densities are affected by modification to the Yukawa rate. For a τ source, we see that the densities for τ, l (right-handed tau and left-handed third generation lepton doublet) are mostly dominant, but decrease as the the Yukawa rate for the tau is increased, while other particle species increase in density. For a b source, it is b, q (right-handed bottom and left-handed third generation quark doublet) which are dominant, again showing a mild increase in other particle densities at their own expense as Γ b Y increases. We also have a slight decrease in the density of u, the representative of the light quarks, as these get sourced predominantly by the strong sphaleron, considering the smallness of their Yukawa couplings.
Thus a decrease in the bottom density results in less chemical potential for strong sphaleron interactions and less accumulation of light quarks. Finally, for a t source, we see an increase in the density of every particle species. Interestingly, it is not the left-handed quark doublet that contributes most to the baryon asymmetry via the weak sphaleron, because the strong sphaleron quickly spreads the quark density among the quarks, and q is almost canceled against q 2 + q 1 = −4u. Rather, it is the left-handed leptons, enhanced by large Yukawa interactions of the top, that drive the weak sphaleron into increasing the baryon asymmetry.
The reason all densities increase in the top case is that the top relaxation rate is the strongest source of washout, and we see here that by increasing the Yukawa rate, all other species, which experience much less washout, increase in density. To show that the relaxation rate of the top is responsible for this behavior, we show in Fig. 3  We also note that turning off the Yukawa rate in the symmetric phase and neglecting the Higgs density reverses this behavior, as well as flips the overall chiral excess and hence the produced baryon asymmetry. In this case, we would require a CPV operator with a coefficient of opposite sign. This emphasizes the impact of the kinetic redistribution of densities that occurs in the transport equations. In Table I we provide a summary of the effects seen in Fig. 1 for such typical modifications that may occur given more precise calculation of the relaxation rates.  with the corresponding active source term. Effects of modifications to relaxation rates of species with no active source is smaller than for the particle with the active source.

B. Sphaleron rates
The sphaleron rates are similarly subject to uncertainties [25][26][27]. It is interesting to compare the sensitivities to these parameters between the case of a t source and a τ source.
Introducing similar modifiers, κ ss and κ ws , Fig. 4 shows that the top-sourced BAU is suppressed when the strong sphaleron rate is decreased. The tau, in comparison, is virtually unaffected by modifications to the strong sphaleron rate: an O(10) modification to Γ ss with a tau source changes Y B by about about 0.1% (not shown in the figure). This is because the strong sphaleron acts solely on quarks, which are only weakly coupled to the lepton sector via the Higgs, and therefore have little impact in the case of a lepton source. Changes in the weak sphaleron rate impact the baryon asymmetry similarly for both τ and t, as seen in the right plot of Fig. 4.

C. Ultra-thin wall approximation
Approximating the relaxation rate Γ M as a step function requires choosing the point where it is turned on/off, which is essentially choosing the position of the bubble-wall. This is the ultra-thin wall approximation, and is a necessary step in the matrix formalism (see [17,21], and also [14] for a direct comparison between the characteristic bubble wall width L w and other typical length scales). This choice is somewhat arbitrary, since the actual bubble-wall has a smooth profile characterized by φ b (z) (see Eq. (A11)). Two sensible choices would be placing the wall at z = 0, the center of the bubble profile, and shifting it by its characteristic width to z = −L w . In this section, we estimate the impact of this choice. In Figures 5 and 6, we plot the baryon asymmetry obtained by shifting the point chosen for the step function. We overlay the plot of Y B as a function of the wall shift with the shape of the source, which is maintained in this approximation, and with the shape of Γ M .  We can see that shifting the wall to the right quickly eliminates the generated baryon asymmetry. This is because the source is truncated: at a shift of +0.5 GeV −1 , there is virtually no source left in the broken phase (recall that the source is taken with its z dependence, but taken as active only in the broken phase), and hence no baryon asymmetry. For negative shifts, the source is fully present in the broken phase, but we also overestimate the relaxation rates by taking the approximating step functions to be active in regions where the corresponding Γ M 's are in fact already highly suppressed. This explains the decrease in Y B for negative shifts. The exact position of the peak is set by the competition between the inhibitory effect of overestimating Γ M and the enhancement by including more of the source. We find that the variation in the predicted Y B between placing the wall at z = 0 and z = −L w is ∼ 5% for τ , ∼ 20% for t and ∼ 50% for b.

D. Bubble wall thickness and velocity
Successful EWBG requires a strong first order phase transition. The details of the phase transition and the subsequent bubble nucleation and growth are important features that for each specific model will determine important parameters such as the wall thickness and wall velocity. Such a study is beyond the scope of the present paper, we refer the reader to recent analyses [28,29]. In our approach, we estimate the impact of modifying the bubble wall parameters: its wall velocity v w and thickness L w . The wall velocity can directly impact the diffusion time scale for successful baryogenesis and the validity of the two-step approach. In Fig. 7, we plot the baryon asymmetry as a function of the bubble wall velocity for each source, while in Fig. 8,  at v w = 0.05, while T b I is arbitrarily normalized since it cannot produce the observed asymmetry. is the same as in Fig. 7.
changes of the predicted asymmetry in response to changes in the wall velocity and width.
The asymmetry from a tau source is less affected by L w , varying only mildly from Y τ whereas the topand bottom-sourced asymmetry depend more strongly on L w , with a similar slope for t and b. The change of sign in Y t B (v w ) is yet another aspect of the sensitivity of the top source to model parameters. While the benchmark value of v w = 0.05 is near-optimal for the t-and τ -sources (cf. also Ref. [14]), the formalism of Ref. [29] beyond the small-v w approximation shows that high yields of Y B are also possible for larger v w . For large L w , the ultra-thin wall approximation (taking Γ f M as step functions) might also become less accurate, although important length scales as migration, diffusion and interaction lengths, as defined and discussed in Ref. [14], are still larger than L max w = 1 above.

VII. CONCLUSIONS AND DISCUSSION
For the calculation of the baryon asymmetry of the Universe in electroweak baryogenesis, we developed a simple and useful method for solving the transport equations semianalytically as well as fully analytically by exploiting various aspects of the structure of the set of these differential equations which couple the participating particle species. We obtained a physical picture of diverging and converging modes and identified the zero-modes as crucial components for the possibility of generating a non-zero baryon asymmetry. Maintain-ing the analytical form allowed us to identify important features and analytic dependence of the baryon asymmetry on model parameters.
While the derivation of our method is general, for the numerical evaluation we calculated the baryon asymmetry within the SMEFT framework with complex Yukawa couplings of the third-generation fermions and the muon. We analyzed how modifications of model parameters and rates affect the resulting baryon asymmetry. This allowed us to estimate the sensitivity of the baryon asymmetry to changes by a factor of O(10) that may result from more precise calculations of these parameters and rates. This large factor is chosen as a conservative example of modifications.
An important feature of our method is that it is straightforward to implement and avoids possible instabilities by the analytical reduction of the system before numerical evaluations are performed. We confirmed the robustness of our method by the following consistency checks: • Robustness to small changes in model parameters, such as the velocity and thickness of the bubble wall, as well as variations of the relaxation, Yukawa and sphaleron rates, with sensible dependence on the parameters. For reasonable values of the model parameters, we find no pathological behaviors. Furthermore, we investigated the impact of the ultra-thin wall approximation by varying the assumed position of the bubble wall.
• Convergence of the one-and two-step approaches (that differ by the inclusion of the weak sphaleron rate in the transport equations) in the limit of a small weak sphaleron rate; with a relative difference of ∼ 4, 15% for the τ , t, respectively, at the nominal weak sphaleron rate.
• Good agreement between the semi-analytic and fully-analytic results in all the scenarios that can be tested with the fully-analytic method. The relative deviation remains below O (10 −11 ) for approximating all diffusion coefficients equal, and below O (10 −4 ) for distinguishing between a quark and a lepton diffusion constant.
• Derivatives of particle densities receive the correct coefficients in the eigenvectors: precisely an extra factor of the eigenvalue, as expected by exponential solutions, up to relative differences of O (10 −5 ).
• Summing over particle densities confirms B − L conservation up to relative deviations of O (10 −5 ) or better.
• Our method produces consistent results (within less than 1%) whether we incorporate or neglect light particles, as physically expected. This implies that it does not suffer from the increase in computational complexity when enlarging the K matrix. We checked this consistency by explicitly solving the transport equations for various setups of the full SM fermionic sector, which we used to produce the muon results in Ref. [23].
We conclude that the main conclusions presented in our previous works [16,23]  In this Appendix we present the expressions and values for all parameters required to fully reproduce the final results.

Benchmark parameters
We take the nucleation temperature to be T N = 88 GeV. At this temperature, the gauge couplings and Higgs VEV are given by [21] g = 0.36 , g = 0.65 , g s = 1.23 , v N = 152 GeV .
The entropy density, written in terms of the temperature and the entropy degrees of freedom g * , is given by [21] s = 2π 2 45 g * T 3 N , g * = 106.75. (A2) The bubble wall velocity and width are taken from [14], with values v w = 0.05 , L w = 0.11 GeV −1 .
The diffusion coefficients are approximately given by [11,15] D l R = 380/T ,

Thermal properties
The real part of the thermal mass of a particle f is of the form where g i are the gauge couplings and c i are combinatorial coefficients. We denote a left (right) handed lepton by l L(R) , a left-handed quark doublet by q and a right-handed up (down) type quark by u (d). The thermal masses are given by [30] Re[δm 2 l R (T )] = The k-functions related to the chemical potentials in Eq. (2.2) are calculated as [21] k wherek f counts the physical degrees of freedom in the multiplet (e.g.k q = 6,k H = 4), c F (B) = 6(3), and +(−) is chosen for fermions (bosons).
Next, we define and (A10) These are used to calculate the CPV source and the CP-conserving rates. We use the kink solution as a typical ansatz for the space-dependent Higgs VEV:

Source and CP-conserving rates
The CP-violating source is proportional to the relative phase between the mass and its spatial derivative. Explicitly, the source is given by the expression [19,20] For the relaxation and Yukawa rates of the CP-conserving processes, we neglect hole modes to get where , In Table II Table II, get corrected according to Here The sphaleron rates are estimated via lattice calculations, and are given by [32,33] Appendix B: Consistency checks

B-L conservation
A simple and important check using the one-step approach is to verify that B − L is conserved. We define the relative difference between the baryon and lepton numbers, as In Table III we show R B−L for each of the four fermions of interest, setting T f R = 0, T f I = ±0.05 (+ for t, -for b, τ, µ) in each case, and the rest of the dim-6 operators to zero. We find that across the parameter space the relative difference does not exceed ∼ 10 −5 .

Derivative test
We construct our solution as a set of 1st order differential equations. Thus half of the entries are the first derivatives of the various particle densities. Recalling that the solutions are exponents, the entries of the derivative terms in each eigenvector should be the same as those of the corresponding particles, multiplied by the appropriate eigenvalue. In the fully analytic case, the equality is exact. In the semi-analytic case, we define the relative difference between a derivative entry and the particle entry times the appropriate eigenvalue as Here f i denotes the entry of the i'th eigenvector corresponding to particle f . In Table III we show R f ,f for each dim-6 operator, in the broken phase. In each case, R f,f denotes the largest value of all particles and all eigenvectors.

Number of particles
Increasing the number of participating particles in the transport equations may result in numerical instabilities. Physically, very light particles should not affect the resulting prediction for the baryon asymmetry, and are typically neglected. Light quarks participate in strong sphaleron interactions, which are efficient at high temperatures. However, in the approximation that first and second generation quarks are massless and weakly interacting, they are degenerate in the transport equations, and we may choose a single representative to capture their contribution (see Sect. II for the explicit treatment; see also Ref. [14]). We verified that our method is robust to changing the number of participating fermions. We quote in Tab. III the resulting Y B in two scenarios: 1. A set containing t, b, τ, u, h, where u is a representative of the light quarks, which we used to produce the results in [16]. Note that the muon does not appear here. We denote this scenario as Y t,b,τ B .
2. The full SM set, as used to produce the results for the muon [23]. We denote this as We set T f R = 0, T f I = ±0.05 (+ for t, − for b, τ, µ), and the rest of the dim-6 operators to zero.  Table III: For each fermion: Relative difference between baryon and lepton number (first row) and between an eigenvector of a particle density and its derivative (second row); prediction of Y B with t, b, τ, u, h (third row) participating in the transport equations, and with all SM particles including µ (fourth row), for T f R = 0, T f I = ±0.05.

Comparing the one step and two step approaches
Varying the weak sphaleron rate can be used to compare the one-and two-step approaches. For the full SM set of transport equations, we parameterize the difference in their predicted baryon asymmetry as a function of the weak sphaleron modifier κ ws , by where Y 1 B , Y 2 B are the predicted baryon asymmetries as obtained in the one-and two-step approaches, respectively. Since the one-and two step solutions differ by the inclusion of the weak sphaleron rate in the transport equations, we expect the two approaches to converge as we decrease the rate of the weak sphaleron, i.e. A 21 Γws→0 − −−− → 0, and indeed we see this behavior in Fig. 9.  For the benchmark value Γ ws (κ ws = 1) we obtain which corresponds to a deviation of ∼ 4, 15%, respectively. At κ ws = 0.01, we find The relative difference is small at the literature value corresponding to κ ws = 1. It grows for large values of Γ ws , but remains below 30% for the t and 60% for the τ at κ ws = 10. Hence the two-step approach still reproduces the order of magnitude of Y 1 B even in the extreme case of such a large correction factor of the weak sphaleron rate.

Comparing the semi-analytic and fully analytic methods
Under certain approximations, we can solve the transport equations analytically. Below is a comparison of the semi analytic method to the fully analytic method in scenarios where it is applicable. We consider the two scenarios described in Section V: • Case 1 : We neglect the Higgs density, and approximate the diffusion coefficients as equal among all fermions, left and right, setting D f = 100/T .
• Case 2 : We neglect the Higgs density, but the diffusion coefficients are taken to be equal separately among quarks and leptons, i.e.
In Table IV we show the largest deviations in the eigenvalues λ i between the semi-analytic and the fully analytic solution, eigenvector entries φ ij and baryon asymmetry Y B . All these quantities are calculated for a tau source 2 . We define R λ , R φ , R Y τ B , respectively, similarly to the definition of the B − L conservation and derivative test.  We see that the error remains small throughout the calculation, and the resulting baryon asymmetry is in good agreement between the semi-and fully-analytic procedures. We also checked the case D l = 380/T and D τ = 100/T , following the direct calculation outlined in Section V. The results are very similar to Case 2, thus we do not show them here explicitly.

Appendix C: The baryon asymmetry
For completeness, we present the derivation of the expression for the baryon asymmetry, assuming the chiral density n L has already been solved for. The solution for the transport equations is then a straightforward generalization to this procedure.
We approximate the dynamics of the baryon density, n b , by a one-dimensional differential equation in the bubble wall frame, placing a planar wall at z = 0, with the broken phase chosen to be z > 0. Using the diffusion approximation, similarly to the transport equations, the equation for the baryon density n b is where v w is the bubble wall velocity, D q is the quark diffusion coefficient, Γ ws is the weak sphaleron rate and R = 15/4 is the so-called SM relaxation term. The sphaleron process is efficient only in the symmetric phase [27,32,33](assuming a strongly first order phase transition), and we therefore take the sphaleron rate to be a step function Γ ws → Γ ws θ(−z), where Γ ws is constant. All other coefficients are constant numbers as well, and the chiral density acts as an external source for the baryon number density, which, due to the sphaleron rate, is active only in the symmetric phase.
In the broken phase, the solution to Eq. (C1) is of the simple form The particular solution is obtained by variation of parameters n p b (z) = K 1 (z)u 1 (z) + K 2 (z)u 2 (z) .
Using the Wronskian W = e α − z e α + z α − e α − z α + e α + z = (α + − α − )e (α + +α − )z we solve The particular solution is thus given by Let us impose boundary conditions. In the broken phase, the second term in Eq. (C2) diverges as z → ∞, and we set A 2 = 0. The baryon number density is thus completely determined by A 1 and therefore by the continuity condition at z = 0. In the symmetric phase, the second term in Eq. (C3) vanishes as z → −∞, but the first one diverges. In Eq.
(C7), the first term vanishes while the second terms diverges. The divergence of the second term is manifest, and in App. C 1 we show that the first one does indeed vanish. For the divergent term, we set The second line is identically zero, and the first vanishes at z = 0. From this we obtain the continuity condition for n b , n b at z = 0: and finally we obtain