{'fraction_non_alphanumeric': 0.02542534160545046, 'fraction_numerical': 0.00477676702317969, 'mean_word_length': 4.618798118049615, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the wall velocity dependence of electroweak baryogenesis

We re-evaluate the status of supersonic electroweak baryogenesis using a generalized fluid Ansatz for the non-equilibrium distribution functions. Instead of truncating the expansion to first order in momentum, we allow for higher order terms as well, including up to 21 fluctuations. The collision terms are computed analytically at leading-log accuracy. We also point out inconsistencies in the standard treatments of transport in electroweak baryogenesis, arguing that one cannot do without specifying an Ansatz for the distribution function. We present the first analysis of baryogenesis using the fluid approximation to higher orders. Our results support the recent findings that baryogenesis may indeed be possible even in the presence of supersonic wall velocities.


Introduction
Despite the immense empirical successes of the Standard Model, the origin of the matter-antimatter asymmetry in the Universe remains one of the outstanding open problems in particle physics. The issue is particularly tantalizing because our current theories do seem to contain all the necessary ingredients for solving the problem, namely baryon number non-conservation, charge and charge-parity violation, and non-equilibrium dynamics. Indeed, in the Standard Model baryon number is violated by processes which become very efficient at high temperatures, such as those found in the early Universe, where one also has non-equilibrium dynamics such as the Hubble expansion and possibly phase transitions, while C and CP violation are present in the weak and fermionic sectors. However, the requirement that these ingredients appear simultaneously and in a sufficient amount turns out to be a severe constraint, making it difficult for a concrete mechanism to be devised that correctly predicts the observed asymmetry.
A very appealing approach is non-local electroweak baryogenesis [1][2][3], which has attracted much attention from the community in the past decades for its phenomenological and cosmological implication at current and near-future experiments. The mechanism relies on charge transport along bubbles of the Higgs field in a first-order electroweak phase transition in the early Universe. The basic idea is that particle collisions in front of the expanding bubble may lead to a net CP asymmetry, which can be converted to a baryonic asymmetry by baryon-number-violating sphaleron processes in front of the wall. When the asymmetry is diffused into the bubble, where the sphalerons are inefficient, the reverse baryon-wash-out reactions are suppressed, ensuring that the Universe remains asymmetric to the present day.
Because this mechanism relies on diffusion of particles in the primordial plasma, it has been common lore that it could not work for supersonic wall speeds, since in this case the time for non-equilibrium processes to take place in the plasma before the bubble sweeps past would supposedly be insufficient. This intuition was confirmed in the first works that studied the velocity dependence of electroweak baryogenesis in detail [4].
However, this claim has been disputed in a recent work by Cline and Kainulainen [5]. The main argument is that diffusion, being an essentially microscopic non-equilibrium process, heavily dependent on particle interactions, is not necessarily related to a collective phenomenon such as the propagation of sound in the plasma. Indeed, the Bose-Einstein and Fermi-Dirac equilibrium distribution functions show no regard for the sound speed c s 1/ √ 3, and some of the particles do have larger velocities and could diffuse into a supersonic bubble. This is confirmed by their analysis of the relevant non-equilibrium transport equations, which shows that the baryon asymmetry behaves completely regularly at the speed of sound.
However, these results crucially depend on the Ansatz for the shape of the non-equilibrium particle distribution function. In this recent study [5] the working hypothesis was (based on previous baryogenesis analyses [4]) that these functions are fully characterized by two perturbations only: a chemical potential and a velocity parameter encapsulating how the particle's velocity deviate locally from the plasma velocity u µ . All other possible perturbations are assumed to be linearly related to these velocity perturbations. Apart from being a restrictive and largely unjustified assumption, this factorization hypothesis does not sufficiently specify the shape of the distribution functions to allow for a computation of the relevant collision terms. For these, one has to assume a specific Ansatz, e.g. the fluid Ansatz of the form and expand the perturbations δ in powers of momenta, truncating the expansion at some appropriate order [6,7]. Thus, the working hypothesis of reference [5] is actually inappropriate for practical purposes, since it forces one to use different Ansätze for different terms of the Boltzmann equation.
When the transport analysis is performed thoroughly and consistently with this Ansatz, keeping only terms up to first order in momenta (a so-called fluid Ansatz which amounts to including a chemical potential, a velocity perturbation and a temperature fluctuation as well), the speed of sound emerges as a singularity of the Liouville operator, corresponding to the kinetic term in the Boltzmann equations. This implies that, for supersonic walls, all fluctuations trail the source [8]: there are no non-equilibrium dynamics in front of the wall, and the baryon asymmetry should therefore vanish. This is in sharp contrast to recent findings in [5], but still does not constitute a definitive position on the issue, because one must still ask whether this first order truncation is justified at all, and how the system would change if higher order terms were included.
The goal of the present paper is to extend this analysis of the transport equations using an extended fluid Ansatz, including more perturbations appearing at higher orders in momenta. To keep consistency, collision terms from annihilation and scattering processes associated to these new fluctuations are computed numerically and analytically. Since the issue under investigation is encapsulated in the behavior of the kinetic term, we choose to simplify the analysis by neglecting complications arising from CP-violation as well as neglecting collision terms that would be relevant to a full baryogenesis study, such as chirality flips, W -boson interactions and strong sphalerons.
Our treatment is rigorously adequate for the study of friction in front of the wall, a problem which has been tackled in a fully relativistic approach both with the firstorder fluid approximation [8] and in the formalism adopted in [5]. But baryogenesis and wall friction are very similar phenomena, since both rely on out-of-equilibrium distributions of the particles in the plasma close to the propagating wall, and in both cases the source that drives the distribution functions away from equilibrium is the interaction between the particles and the Higgs. The main difference is that in baryogenesis one is mostly interested in the CP-violating component of the deviation from equilibrium, but for the purposes of discussing the behavior of the solutions around the speed of sound this is not relevant. Interestingly, our findings support the results in [5] that the fluctuations behave continuously across the sound speed, albeit for different reasons.
The paper is organized as follows. In section 2 we review the standard fluid approximation for the non-equilibrium distribution function [6,8], and how the Boltzmann equation can be used to determine the three fluctuations in chemical potential, local particle velocity and local temperature. We will see how the speed of sound emerges naturally in this setup as the singularity of the kinetic term. In section 3 we present some arguments leaning towards a criticism of this simple approximation and motivate a generalization with more fluctuations. We also present here some specific criticisms on the standard approach for transport baryogenesis found in the literature [4,5,9], arguing that one cannot escape from making a specific Ansatz for the shape of the distribution functions. In section 4 we present our generalized fluid Ansatz, showing explicit expressions for the kinetic matrix, the collision terms and the source up to six fluctuations. Our results are presented in section 5, and we conclude in section 6. An argument from hydrodynamics for the relevance of the speed of sound is presented in Appendix A. General expressions for the collision terms to leading-log accuracy up to arbitrary orders in the generalized Ansatz can be found in Appendix B.

The fluid approximation
Ultimately, we want to study the Boltzmann equation, that in its relativistic form reads where p µ is the four-momentum of the particles (evaluated on-shell p 0 = E = p 2 + m 2 ), f is the particle distribution function, F denotes the forces in the system and C is the collision term.
The Boltzmann equation is a non-linear partial differential equation (due to the collision term), and in order to make progress different approximations are typically used. The first approximation is to assume that the system is close to equilibrium. The equilibrium distribution is given in terms of the four-velocity of the fluid u µ (relative to the global restframe of the plasma) and the temperature T . In the current setup, the forces will drive the system out-of-equilibrium while the collision terms will relax the system to some local equilibrium that in principle can be different on both sides of the wall (in fact it has to be so, as demanded by energy momentum conservation, see Appendix A). In the analysis of baryogenesis, this problem can be eliminated by introducing a background fluid (made up mostly from the gluons and quarks that are not much affected by the Higgs) that will reflect this change as demonstrated in [6].
Since we assume the system close to equilibrium, one can split the distribution function into a deviation and an equilibrium value After this one arrives at a system of equations that reads where the source term S contains all the forces and depends at leading order only on some forces acting on the equilibrium distributions.
Here we want to note that any peculiarities at the speed of sound are due to a vanishing (or very small) eigenvalue in the first term of this equation (the so-called Liouville term). In particular, it has nothing to do with the specific form of the sources or collision terms. This is also why the relativistic analysis for the wall friction can be carried over to baryogenesis without much effort.
In the following we summarize the steps in ref. [8] to present the analysis for the fluid approximation in Lorentz covariant form. The fluid approximation assumes the different species have individual temperatures and fluid velocities and also chemical potentials. Assuming that the fluctuations are small relative to the background, one obtains δf − µ/T + p µ δu µ /T − p µ u µ δT /T 2 f eq (p µ u µ /T ) .
The fluctuations µ, δT and δu µ describe (for every species) the deviation from equilibrium. Note that δu µ obeys u µ δu µ = 0 so that, in the planar wall case we will consider, it constitutes only one degree of freedom. Therefore in this approximation the system is described by three fluctuations in total. Before we come to the usual approach of solving these equations, let us discuss whether this Ansatz makes sense. The main argument here is that (depending on the system under consideration), there often is a hierarchy between different classes of interactions. The strongest interactions are often gauge scatterings that do not change the particle content of the plasma. In the Standard Model, these are e.g. 2-by-2 scatterings mediated by the gluons. Notice that the t-channel interactions even feature a logarithmic enhancement on top of the rather large gauge coupling in the strong sector [7,10]. When inspecting the collision terms, these interactions will equilibrate the quarks among themselves. Very quickly, all quarks should have a distribution function close to the fluid approximation and share the same temperature and velocity but not chemical potential, since this requires interactions that change particle number.
This leads to the conventional wisdom that kinetic equilibrium is attained faster than chemical equilibrium. (Notice that there is a subtle difference in the definition of the chemical potentials of particles and anti-particles when it comes to the study of the friction compared to baryogenesis, see section 4). Still, ultimately, there is no proof that the fluid approximation is holding in the baryogenesis and/or friction setup and that differences between the fluctuations in the temperatures and flow velocities should equilibrate slower than any other fluctuations. This leaves the question of how important these fluctuations are for a correct determination of the friction or baryon asymmetry. This is an open question that was already raised in the seminal work by Moore and Prokopec [6], and is discussed to some extent in their Appendix B.
Even though one cannot make a strong argument in favor of the fluid approximation with only three perturbations, let us briefly discuss the standard procedure in order to understand the generalization which will follow, as well as to understand where the strong dependence on the wall velocity in the setup comes from. Since the fluid approximation has three degrees of freedom per species, one can take three moments of the Boltzmann equations to obtain a (non-partial) system of differential equations. The most natural choice is to simply multiply by four-momenta and integrate, so that the Boltzmann equation (4) reads Note that the denominator includes only one power of E from the Lorentz invariant integration measure. At this point our treatment differs from [5] (and also from [4,9]) who advocate different moments in their analyses, obtained by multiplying the Boltzmann equations with extra factors of 1/E as well. Our choice of moments seems more natural, since the resulting equations can be easily interpreted as the divergence of the charge current and the energy-momentum tensor. As we will see, this is one important reason why the speed of sound is of relevance in the analysis of [6], while it is not in [5]. Obviously, the analysis presented here does not rely on small wall velocities. We can then linearize the system using the fluctuations in (5) and focus on the case of an approximately planar wall, which is justified after the bubble reaches a steady-state, and also because we are interested in the behaviour of fluctuations close to the wall where its curvature becomes negligible. This means that there is only one directionū µ perpendicular to the fluid velocity, i.e. withū µ u µ = 0 andū µū µ = −1. In particular the velocity fluctuations can be written as δu µ =ū µ δv. Moreover in the steady-state situation the fluctuations depend only on ξ = x µ v µ , where v µ = γ(v w u µ − u µ ) characterizes the four-velocity of the wall, which then leads to ∂ µ = v µ ∂ ξ . Finally, we note that equation (7) contains two components, which we can project onto u ν andū ν , leading to a total of three equations for three fluctuations. Putting everything together, we see that the Liouville part of the Boltzmann equations, corresponding to the kinetic component of the flow, involves derivatives ∂ ξ (...) of the fluctuations multiplied by coefficients of the form where, in the massless case and for n ≥ 2, with f BE and f FD the Bose-Einstein and Fermi-Dirac distribution functions, respectively. For the collision terms associated to annihilations, it will also be convenient to define which we quote here for reference. In Appendix B, the collision terms are derived to linear order, so this momentum decomposition of the Boltzmann equation leads to a system of the form with the vector q = (δµ/T, −δT /T, δv) describing the fluctuations of the system. The matrix Γ stems from the linearized collision term while S contains the sources. On the other hand, the matrix A comes from the Liouville term and has the form (for a relativistic plasma in equilibrium) Now, the authors of [6] made the interesting observation that for v w = 1/ √ 3 the determinant of this matrix vanishes and one eigenvalue changes sign. The corresponding eigenvector corresponds to v w δv = −δT /T . For supersonic wall velocities all eigenvalues have the same sign such that all fluctuations trail the wall. Hence the only contribution to baryogenesis comes from a local source in the wall and the resulting asymmetry is strongly suppressed. Notice that this discussion does not rely on the special form of the collision term or the source. Only the analysis of the Liouville term in the Boltzmann equation is relevant to make this argument. Moreover, the effect should occur in every species individually and does not arise from the interplay of different species.
This fits rather nicely with the observations of how the system behaves on the largest scales, as inferred from hydrodynamics (see Appendix A). In particular, the behavior of the system changes qualitatively once the wall velocity surpasses the speed of sound.

Why higher moments are essential
Since the validity of the fluid Ansatz with only three momentum-independent fluctuations is somewhat controversial, let us now venture towards a generalisation of the previous argument. As a consequence of the discussion in this section, we will propose a generalised fluid Ansatz with higher order perturbations, analogous to including a momentum dependence of the fluctuations in the previous section.
Let us start by noticing that sound waves are a collective phenomenon that obviously needs at least some notion of equilibrium. Assuming that the plasma has a rest frame where the distribution functions only depend on the energy, in a general frame the distribution function is of the form f (p µ u µ /T ) where u µ is the fluid four-velocity and the temperature T is introduced for dimensional reasons. Otherwise, the function f is arbitrary. Now consider two fluctuations that mimic the temperature and velocity fluctuations from the last section, and let us again focus on the steady-state situation, with ∂ µ = γ(v w u µ −ū µ )∂ ξ . The Liouville term then contains three different Lorentz structures, As in the fluid approximation, for v w δv = −δT /T the first term vanishes but the others do apparently not. The remaining terms combine into However, for a relativistic plasma, in the plasma frame p µ u µ = E and p µū µ = p z . Moreover, for massless particles E n p 2 z f = E n+2 f c 2 s and odd moments of p z vanish (this follows e.g. from the fact that the energy momentum should be traceless). This means that the lowest moment that can potentially spoil the vanishing eigenvalue in the kinetic term contains four factors of p µū µ . This is one of the motivations to generalize the fluid approximation, to introduce a larger class of fluctuations and to use higher moments of the Boltzmann hierarchy to study the system more rigorously.
Before jumping into this analysis, let us briefly comment on the alternative approach found in the literature, which claims to perform the analysis without an explicit Ansatz for the perturbations. One can quite generally write the non-equilibrium distribution function as letting δf encode all fluctuations away from equilibrium except for the momentumindependent part of the chemical potential, which is encapsulated by µ. In standard transport baryogenesis treatments [4,5,9] one then take moments by multiplying the Boltzmann equation by powers of p z /E in the plasma frame, resulting in fluctuations defined as 1 Of course, this description is so general that it includes the fluid approximation as a particular case. But, without imposing any further restriction on the form of δf , one cannot make any progress towards determining the shape of δf . So some Ansatz is unavoidable. This difficulty is then typically sidestepped by imposing a factorization condition which essentially states that all fluctuations u for ≥ 2 are proportional to u 1 . This entire approach is problematic for a few reasons. First, the choice of momenta is very peculiar. While u 1 is the z-component of the four-current in the wall frame, u 2 is not simply related to the relevant conserved quantities like the energy momentum tensor, and consequently the physical meaning of the resulting equations in this formalism is less transparent. Moreover, using the energy-momentum tensor conservation as a moment, as in equation (7), has the advantage that the corresponding collision terms vanish (not for individual species but collectively).
Next, this approach is supposedly the most general since it does not rely on an Ansatz. But in fact the specific choice of factorisation, relating all higher-order fluctuations to one velocity perturbation u 1 , is equivalent to an Ansatz, with the disadvantage that its physical significance is not very transparent. Besides, this Ansatz is not specific enough to allow for the computation of the collision terms. For that, one needs to know the explicit form of the non-equilibrium distribution, such as we do have in the fluid approximation. Even though this seems unrelated to our issue at hand, this makes the method impractical in actual calculation. In previous baryogenesis analyses [4,5,9] the collision terms have been ultimately derived assuming the fluid approximation.
Finally, as discussed in the last section, for the peculiar fluctuations (14) the kinetic term vanishes up to fourth order in p z at the speed of sound, so by construction u 1 = u 2 = 0 and the system is oblivious to the eigenvector with zero eigenvalue. In other words, the functions u 1 and u 2 do not encode this eigenvector and only higher moments would be sensitive to it.
In summary, we cannot escape from using a specific Ansatz for the non-equilibrium distributions, yet there is little justification in keeping only three momentum independent fluctuations in the fluid approximation. In the following, we will generalize this fluid Ansatz to include terms at higher order in momenta. We will see that the results for baryogenesis across the speed of sound are modified even at a qualitative level.

Generalized fluid Ansatz
A general non-equilibrium distribution can always be written in the form shown in equation (1). But, in much the same way as argued above, this is too generic to allow for the computation of the resulting perturbations from the Boltzmann equation. Some progress can be made by expanding the perturbations δ in powers of momenta. The standard fluid approximation is then obtained by truncating this expansion at first order, as in (5). Here, we would like to extend the set of fluctuations and study the behaviour of the solutions. This is not the first time that higher moments have been considered in treatments of the Boltzmann equation, see e.g. [6,11,12] for an interpretation of the higher-order equations in terms of entropy production due to dissipative work and heat flow. However, the present work is the first that uses this extended approach in the context of baryogenesis.
To be specific, we parameterise the fluctuations as This generalizes the fluid approximation, and due to the symmetries only n+1 degrees are relevant in w (n) . For instance, expanding up to n = 2 will lead to 6 fluctuations in total. In order to truncate the Boltzmann hierarchy we then multiply the Boltzmann equations by the factors (p µ u µ ) a (p νū ν ) b with a, b ≥ 0 and (a + b) ≤ n and take the corresponding moments to arrive at a linear equation system for the fluctuations.
The result is again a system of the form shown in equation (12), but now with enlarged (n + 1)(n + 2)/2 square matrices, whose entries will involve the coefficients defined in (10) and (11). Indeed, for n = 2 the kinetic matrix has the form The source term is defined in equation (4) as and the force can be divided into a CP-even and a CP-odd term as mF µ = ∂ µ m 2 /2 + mFC P µ . As mentioned before, baryogenesis relies on C and CP violation, so the CPodd component of the source is essential in this calculation. However, in this work we are interested in studying the behaviour of the solutions across the speed of sound, which depends exclusively on the kinetic term. In this sense the calculation of the baryon asymmetry and of friction are similar, and we will therefore consider here a CP-even source for simplicity. Then, after multiplying the Boltzmann equation by factors of p µ u µ and p µū µ and integrating, we again get integrals as in equation (8), resulting in for six perturbations. Here c f 1 = log 2 and c b 1 = log(2T /m), whereas c 2 and c 3 are given in equations (10).
Focusing on the CP-even terms has an extra advantage in simplification. For baryogenesis one would have to compute the CP-odd components of the chemical potential and of the other fluctuations. In particular, the chemical potential in baryogenesis calculations is a proper chemical potential, while the fluctuation µ in friction calculations parametrizes a tilt in the distribution function that is equal for particles and anti-particles. Hence, in the collision terms the fluctuation µ in friction calculation will be damped by annihilation processes. In baryogenesis, on the other hand, damping of the chemical potential obviously requires particle number changing interactions that are far weaker and less abundant. This is because, due to the CP violating source of baryogenesis, the chemical potentials have opposite signs for different chiralities and are reduced by the mass terms of the fermions in the broken phase. Moreover, the strong sphaleron can flip the chirality of the quarks and erase their chemical potentials. Therefore, in focusing on the CP-even part of the fluctuations only, we neglect these complications and simplify the analysis of the collision terms as well.
In summary, we consider collision terms from top-quark annihilations and scatterings off other quarks and off gluons. The collision matrix for annihilations is where A = 1 for scatterings by gluons and A = 5/3 for scatterings by quarks. Analytic expressions for the coefficients γ ann and γ scatt at leading-log can be found in equations (58) and (67) in Appendix B. Now, a further complication arises when we attempt to include the sphalerons in this generalized framework. In standard baryogenesis calculations, based on two fluctuations only, the sphalerons couple only to the zeroth-order fluctuation, i.e. that which is not multiplied by any power of momenta, and which is interpreted as the chemical potential of the corresponding field. In principle, in our framework one would have to recalculate how the effective interaction of the strong sphaleron damps the fluctuations in our generalized Ansatz (18). This task is beyond the scope of what we want to achieve here and we will mimic the true collision terms of the sphalerons in different ways.
The strong sphaleron will be mimicked using the collision term from the friction calculation. This will lead to a damping of the chemical potential similar to what the strong sphaleron achieves in baryogenesis. At the same time it is a set of physical collision terms that just occur in a different setting. So, strictly speaking, the convergence properties that we will observe are the ones of the friction network and not the ones of baryogenesis. In order to facilitate transport, we choose the corresponding gauge couplings rather small, such that the charges can travel into the symmetric phase as in the case of non-local baryogenesis.
To mimic the weak sphaleron, we use the standard result for the sphaleron rate, where a ≈ 37 and φ(z) = φ 0 2 (1 − tanh z Lw ) is the bubble profile with wall thickness L w and vev φ 0 at the critical (or nucleation) temperature. We couple the sphaleron to J µ u µ , where J µ = p µ /E is the current of the plasma, since this is after all the moment associated to particle number density and should therefore be interpreted as the full chemical potential [12] 2 . Coupling the weak sphaleron to the zeroth-order fluctuation would be rather arbitrary and would make the result highly dependent on the precise basis of fluctuations. In any case, the coupling to the weak sphaleron is of course not relevant to study convergence of the out-of-equilibrium distribution functions but rather to reproduce the qualitative behavior of baryogenesis for a low number of fluctuations and supersonic wall velocities where the resulting BAU is suppressed.

Results
With this setup, we have all we need to solve the resulting Boltzmann system and find the fluctuations and the toy baryon asymmetry. The latter is obtained by solving with A = 15/2, and κ = 3/(4πc 2 ) is a factor chosen so that the result will agree with the traditional coupling to the chemical potential at zeroth-order. In order to find the current J µ we solve the system If λ i , χ i are respectively the eigenvalues and eigenvectors of A −1 · Γ, the Green's function is 2 The interpretation of the zeroth-order perturbation w (0) as the chemical potential is only valid under the condition that [12] d 3 p E p µ u µ δf = 0, which is not automatically satisfied in our Ansatz. The coincidence is only granted in a two-fluid approximation, which has often been employed in previous baryogenesis studies. and the fluctuations are We see that for positive eigenvalues the fluctuations at z get contributions from the source at z > z, and we say that the fluctuations trail the source. The opposite occurs for negative λ i , in which case the fluctuations are ahead of the source. Because the baryon number is obtained from a convolution of the fluctuations with the sphaleron rate Γ ws , and because the latter is only active in front of the wall, as can be seen from equation (24), the trailing fluctuations produce a negligible baryon asymmetry. Thus, when all eigenvalues are positive, the resulting asymmetry should vanish. This is exactly what happens for supersonic wall velocities in the fluid approximation with three perturbations, which leads to the interpretation of the speed of sound as an upper limit for viable transport baryogenesis. But when we add additional fluctuations, not all eigenvalues become positive beyond this threshold. Some eigenvalues may remain negative and yield a contribution to the BAU, albeit smaller than for lower velocities because some eigenvalue sign flipping does invariably take place. This is illustrated in figure 1, with the resulting toy baryon asymmetry for two values of the relevant coupling entering the collision terms (in this case, the strong coupling α s ). The red line represents the case of three fluctuations, where it is clear that there is no asymmetry beyond the speed of sound. However, once we add more perturbations, this picture changes and a resulting asymmetry does become possible. There also appear other thresholds, related to new singularities of the larger kinetic matrix. Curiously, as one adds more and more perturbations, the sharp drops in these threshold values seem to smoothen out and one approaches a continuous curve, similar to the result obtained in [5].
We highlight, however, that the reasoning for this similar behaviour is fundamentally different. Here the smoothness is an emergent asymptotic behaviour obtained from a well-defined expansion in momenta, rather than from an ad hoc factorization assumption. Furthermore, the speed of sound does not constitute a sharp discontinuity, but it remains a feature of the system, even if hidden in the first few momenta only. Our approach is also thoroughly consistent, in the sense that we use the same Ansatz to compute all the terms in the Boltzmann equation, including the collision terms.
It is also interesting to point out that the result clearly converges as we add more fluctuations. As can be seen by comparing the two plots in figure 1, the parameter determining this convergence is the inverse coupling appearing in the collision terms. This is not unexpected: the stronger the interactions are, the quicker and  more effective the thermalisation processes will be, so we would effectively need less perturbations to describe the distribution functions well.
In figure 2 we extend the velocity dependence of the BAU down to very small values, v w 0.001. As expected, the result becomes highly suppressed and also converges rather quickly in this region. This is also in agreement with the findings in [5]. We note also that the first-order approximation (with three perturbations, corresponding to the red curves) may be an over-or an underestimate of the full result, depending on the wall velocity and also on the value of the coupling in the collision terms.
Finally, the dependence with the wall thickness is displayed in figure 3. As one might expect, the shape of the curve does not change drastically by the addition of new fluctuations, but it is displaced as already seen in figures 1 and 2 above. The convergence of the series is also clearly highlighted in this plot.

Conclusions
In order to analyze electroweak baryogenesis or the bubble wall friction during a cosmological phase transition, a Boltzmann equation has to be solved for the particles in the plasma. Since a full solution of this system of partial differential equations is often (even numerically) not attainable, the main way to progress is to take moments which leads to the Boltzmann hierarchy.
Some assumptions have to be made to decouple the Boltzmann hierarchy and the evaluation of the collision terms basically forces one to choose an Ansatz for the dis-  tribution functions. It turns out that if only a few moments and a few basis elements are chosen, the outcome critically depends on these choices. For some choices, the speed of sound is an important quantity and baryogenesis for wall velocities beyond the speed of sound is insufficient, for others not.
In this work, we use extensive basis sets (that generalize the fluid approximation) and many moments of the Boltzmann equation to study a toy model that resembles the most important features of electroweak baryogenesis calculations. The main stumbling block for this kind of analysis is the evaluation of the collision term that we detailed in Appendix B. We find that for the fluid Ansatz and a low number of moments, supersonic baryogenesis is indeed suppressed. Then again, for a large number of moments, we find that the outcome behaves smoothly in wall velocity and supersonic baryogenesis becomes possible, in support of recent findings by Cline and Kainulainen [5].
The reason for the suppression in case of a few moments beyond the speed of sound is that the eigenvalues in the Liouville operator change sign. For large wall velocities all eigenvalues share the same sign, fluctuations only exist behind the bubble wall and baryogenesis becomes local. Qualitatively, it is the same for higher moments, but the velocity at which all eigenvalues share the same sign progressively moves to v w → 1. This is not too surprising since for a time-like wall velocity vector v µ all fluctuations have to be damped in the positive time direction. So it is not too surprising that all eigenvalues have the same sign when v µ approaches the light cone.
In order to leverage our results in a realistic baryogenesis calculation some ingredients are still missing. The main improvement would be to determine how the strong and weak sphalerons couple to the fluctuations in this generalized Ansatz. These are essential in the baryogenesis calculations since they break chirality and B +L number. While we model these effect in our toy setup, representing these processes in a full analysis would require the evaluation of the corresponding collision/interaction rates for the extended fluid Ansatz.
Another important issue concerns the systematics governing the momentum expansion in the fluid Ansatz, and the criteria for deciding the appropriate order for truncation. On general grounds one can expect the fluid to behave as nearly perfect if the mean free path is much smaller than the relevant macroscopic parameter, namely the wall width L w . This indicates an expansion parameter of the order (ΓL w ) −1 , with Γ denoting some combination of the relevant interaction ratios. This is in line with our findings that the convergence of the BAU gets better as the coupling constant (hence the interaction rates) become larger. As the mean free path increases, non-perfect behaviour such as viscosity and conductivity starts to become relevant, which are modelled by higher moments in the distribution function. All that said, it should be emphasized that a rigorous and a prioristic proof of the convergence of this moment expansion, or of the accuracy of any specific truncation, is still lacking [12]. However, the convergent behaviour in numerous applications of this method, including the results of this work, corroborate the validity of this approach as long as non-linear phenomena (such as shock waves) are not present. A related issue is the accuracy of a fluid-like Ansatz in describing intermediate and low energy collisions. An expansion in powers of p µ will naturally be dominated by the high energy regime, but may be less sensitive to complex angular behaviour of collisions at lower energies. This is not a problem in the leading-log approximation considered here, but may become an issue when attempting a leading order estimate. In this case a different Ansatz may become necessary, perhaps expanding in the energy E and momentum direction p µ /E. Be it as it may, our numerical results show that, to leading-log accuracy, the conventional fluid approximation with three fluctuations seems to work fine at subsonic wall speeds when aiming at an O(10%) accuracy.
Finally, we comment how the present calculation relates to the evaluation of the bubble wall friction. Unlike the baryogenesis calculation, the friction calculation leads to a change in the collective temperature and velocity of the fluid through energy injection. This is captured by fluctuations in a background field that also become singular when the Liouville operator develops zero eigenvalues, leading to a blow-up in friction as explicitly shown in reference [8]. In the baryogenesis calculation, this singularity in the source is counteracted by a divergent damping and no singularities occur in the final baryon asymmetry. In the friction calculation, however , such divergences can be expected on the ground that there is a build up of particles in front of the wall precisely at the speed of sound, similar to a sonic boom effect, which drastically increases the friction in this case. In [17] the behaviour of friction was analysed in the framework laid out in reference [5], and because the speed of sound plays no role in such formalism the findings point to a continuous non-singular behaviour for the friction at all values for the wall velocity. However, in our generalized fluid Ansatz with higher moments it is well possible that the singularities remain and yet others appear, indicating that, in the limit of very high moments, the problem essentially becomes non-linear for wall velocities beyond the speed of sound.

A An argument from hydrodynamics
Even though hydrodynamics does not describe the phenomenon of diffusion, the hydrodynamic analysis of bubble nucleation and expansion can shed some light on the relevance of the speed of sound. This is an old topic already presented in textbooks [13] and we will only summarize the main relevant points to the present discussion.
While the wall is expanding, the interactions of the particles in the plasma with the Higgs field will drive the system out-of-equilibrium. This will lead to the fact that the equilibrium attained before and behind the wall are not the same. On length scales much larger than the bubble wall or diffusion length but smaller than the bubble size, the system is in local equilibrium. The energy-momentum tensor of the combined system (plasma+Higgs) is conserved and the system can be described by hydrodynamics. Moreover, since there is no inherent length scale, the system behaves self-similarly and the hydrodynamic equations read [14] ( where v is the fluid velocity in the plasma frame, γ is the corresponding Lorentz factor, e and p are the energy and pressure densities and ξ is the self-similar coordinate ξ = r/t. At the same time, in the vicinity of the wall, the Higgs field injects energy (and pressure) into the plasma which leads to the matching equations where u ± denote the plasma velocities (in the wall frame) and T ± the temperature in front/behind the wall.
The interesting point about these equations is that only certain classes are valid globally. For example, if the wall velocity is supersonic (and the phase transition is not too strong), the fluid velocity in front of the wall has to vanish, since otherwise (29) would imply a singularity in the fluid velocities somewhere in front of the wall. So the two matching equations abide to v + = v w and T + is the phase transition temperature. This leads to detonations (v − < v + ) with a rarefaction wave developing behind the wall.
On the other hand, if the wall velocity is subsonic, the fluid velocity behind the wall has to vanish since otherwise the solution will encounter a singularity according to (29). So the only valid solution in this regime is that of a shock before of the wall which are called deflagrations (v − > v + ). The fluid will then drop to the equilibrium configuration in the shock front (see ref. [14] for details).
So, these solutions fulfill all the naive expectations. For supersonic wall velocities, the plasma changes only behind the wall while it is still in equilibrium in front of the wall. For subsonic wall velocities, particles have to be reflected which leads to a snowplow effect that ultimately will build the shock in front of the wall. Qualitatively, the solutions for subsonic and supersonic walls behave quite differently and one would expect that this also has to be reflected by the behavior close to the wall once one looks into the details of diffusion.
Actually, there is also an intermediate regime where so-called hybrids develop with rarefaction waves and a shock. However, for very weak phase transition, this case becomes less and less relevant. Moreover, the strongest detonations (with the smallest wall velocity) are so-called Jouguet detonations. In this case, the fluid profile becomes very steep behind the wall (dv/dξ → ∞) due to the fact that the local fluid velocity (in the wall frame) is the speed of sound, so the Mach number is 1.
The main takeaway from this section is that the macroscopic behavior of the fluid indeed changes from subsonic to supersonic wall velocities, which suggests that also microscopically the picture has to change qualitatively in this transition. Still, it is only very indirect evidence that baryogenesis is not possible for supersonic walls.

B Collision terms
As explained in the last section, we will use the collision terms of the friction network to mimic the damping from the strong sphaleron in baryogenesis. We follow closely the analysis in [6] and also use for most parts their notation and conventions. The collision integrals are of the form 3 with and the shorthand k = d 3 k/2E k . The incoming momenta are denoted k and p and the outgoing k and p . M is the matrix element of the process and f i are the particle distribution functions of the particles (that are fermionic/bosonic depending on process).
We will evaluate the collision terms to linear order in the fluctuations (see [6]), i.e. where and 1 ± f eq p = exp(p µ u µ /T )f eq p . (37) Notice that this relation together with energy-momentum conservation in the process implies that the last factor in P[f i ] is actually symmetric under exchange of k with p, exchange of k with p and exchange of k,p with k ,p as long as the statistics of these particles are the same. In effect the symmetries of P[f i ] depend then on the first factor δ i . Following [6] we are only interested in contributions that are logarithmically enhanced due to IR sensitivity, which only can arise from the t and u channels. The mass dependence of the particles in the matrix element regulate the IR sensitivity of these integrals but we assume the particles to be massless otherwise. There are two types of diagrams we need to evaluate: annihilation diagrams and scattering diagrams. For the incoming particle with momentum p we consider only fermions (quarks) and the scattering can happen off gauge bosons or other fermions.
Our Ansatz for the fluctuations we call collectively χ and expand and so on. In order to obtain the various moments of the Boltzmann equation, we multiply the Boltzmann equation with some factors χ a (p) and integrate over p. The outcome is Annihilations: The fluctuations are in the incoming particles while the produced particles are assumed to be in equilibrium. In case of the friction analysis the dominant contributions are CP and C conserving such that the incoming quarks The matrix element for annihilations is (we use the QCD rates of the Standard Model) We can symmetrize the expression by exchanging simultaneously p with k and p with k what essentially gives Remember that the factors like χ(p) will contain some products of the momentum with an tensor structure for the momenta. This will make the evaluation of these integrals somewhat cumbersome. Below we discuss three methods to deal with this Lorentz structure.
Scattering processes: We can follow the same steps and obtain expression like where the matrix element is for scattering off a quark and for scattering off a gluon. Symmetrisation leads to the replacement In the following we discuss three different ways to solve these integrals. The first one uses a symbolic solver to resolve and invert the Lorentz structure of the integrals. The second method is to evaluate the integrals numerically. The last one is to evaluate all integrals explicitly.

B.1 Symbolic solver
In this subsection we explain how to resolve the Lorentz structure arising from the basis functions χ a (p) using a symbolic solver. We wrote a python code to automatically follow these steps. At some point the procedure becomes to expensive but we obtained results up to fourth order that we compared with ref. [6] as well as with the explicit solution found below.
As noted in [6], one can further simplify the integrals in (44) and (41) by only capturing the leading log. In particular, in the limit t → 0 the difference (p µ −p µ ) and (k µ − k µ ) are of order √ t. This means that since the Matrix elements in combination with the basis functions χ behave as 1/t, one can send k µ → k and p µ → p in the remainder of the integrals. Notice that this is for the scattering integrals only possible after the symmetrization (47).
Annihilations: At this stage, the only dependence on the momenta p and k appears in the matrix element and the Dirac delta function. Without loss of generality one can boost into the frame where p + k = 0, and integrate over the spatial part of the delta function. This enforces p + k = 0. Since we assume the particles to be massless (which implies E k = E p = | p | = | k | = p and E k = E p = | p| = | k| = p), the remaining integral is of the form (see [6]) Here cos θ denotes the angle between p and p and the final result was written in a Lorentz invariant way using p µ k µ = 2pk, which holds in the center-of-mass frame. The next step is to evaluate the Lorentz structure. Ultimately, the integrals only depend parametrically on the vector u µ , so the final result can only involve the vector u µ and the metric η µν . We construct the most general basis that has the correct symmetries and is build from these two ingredients. We then contract the Ansatz as well as the integral with the different basis elements and invert the system.
The resulting integrals can be evaluated in the plasma frame. The expression then involves the energies E p = p and E k = k as well as the Mandelstam variable s = 2kp(1 − cos θ), where cos θ parametrizes the angel between k and p. The integral over the angle reads in leading log approximation. The remaining integrals factorize and using the approximation [6] dp p n log p T f p (1 ± f p ) log(n + 1/2) dp p n f p (1 ± f p ) , and the final integrals can be evaluated depending on the statistics of the particles, yielding the coefficients in equations (10) and (11).
Scattering: The scattering diagrams are a bit harder to determine. The basis functions χ depend also on p such that the integrals over p and k are non-trivial. Still, also this problem can be solved by choosing an appropriate basis and inverting the problem by taking contractions of the integrals with this basis.
The most general outcome after integrating p and k can contain a tensor structure involving k µ , p µ and the metric η µν . Consider the contractions of this basis with terms of the form (51). The contraction is at least of order t such that the overall integrand scales as t −1 just as the annihilation contributions. The Mandelstam variables are given by 2k µ p µ = s , 2k µ p µ = −u = s + t , 2p µ p µ = −t .
First consider basis elements that contain factors p and k but not η. The corresponding contractions with (51) are at least of order t 2 and can be discarded. Basis elements with more than one η are also higher order in t. The only contributions linear in t arise from basis elements with one factor η and the remaining factors are all k. We call the set of these basis elements C.
To invert the problem, one has to consider again the most general basis elements involving k, p and η, and considers the class of basis elements that can give a nonvanishing contraction with elements of C (meaning not suppressed by any factors t). One obvious class are elements with one factor η and the remaining factors are p. But there is a second class: elements with one factor k and the remaining factors are involve only p.
The fundamental integrals over p and k that have to be solved are then of the same form as for annihilation processes and the resulting expressions have a Lorentz structure in k and p that can be resolved in the same way as for annihilations. One novelty is that the scattering on gauge bosons involves integrals with two bosonic distribution functions as given in (10).

B.2 Numerical integration
Another method to evaluate the collision integrals is to do it numerically, preferably with a Monte-Carlo integrator like CUBA [15]. One might think that this even allows to evaluate the integrals beyond the leading log approximation, but this is a mirage since the true result going beyond leading log requires to incorporate hard thermal loop correction in the scattering processes [7,10,16].
The only nontrivial problem in the numerical evaluation is how to represent the four-dimensional Dirac delta function. We choose to align the wall along the zaxis and then sample the six momenta k and p. The vectors k and p are then parameterized as k = P + Q , p = P − Q , where P = ( k + p)/2 and we choose Q = qq = (0, q sin α, q cos α). The remaining constraint on the energy then gives with 2Ē = E k + E p . Due to the symmetries of the integral and since we sample all signs of the six momenta k and p, one can restrict oneself to the positive branch of q. Notice that this construction leads to an additional factor from the Jacobian determinant when the integral of the delta function is performed, We checked that the leading log result corresponds to the one obtained with the other two methods. The full numerical result can differ from the leading log up to a factor 2 in extreme cases (depending on the involved masses that regulate the IR behavior). We checked that this has no effect on our main results.

B.3 Explicit solutions
It turns out that, to leading log order, the collision integrals for annihilation and scattering can be solved analytically to a simple closed form. odd, while for r + s even one has γ mnrs scatt = log 5(2m + 2n + 2r + 2s + 1)T 2 m 2 g × × ms + mn + nr r + s + 1 The term c f 2 enters in scatterings by fermions and c b 2 appears in scatterings by gluons.