Bubble wall dynamics at the electroweak phase transition

First order phase transitions could play a major role in the early universe, providing important phenomenological consequences, such as the production of gravitational waves and the generation of baryon asymmetry. An important aspect that determines the properties of the phase transition is the dynamics of the true-vacuum bubbles, which is controlled by the density perturbations in the hot plasma. We study this aspect presenting, for the first time, the full solution of the linearized Boltzmann equation for the top quark species coupled to the Higgs field during a first-order electroweak phase transition. Our approach, differently from the traditional one based on the fluid approximation, does not rely on any ansatz and can fully capture the density perturbations in the plasma. We find that our results significantly differ from the ones obtained in the fluid approximation (including its extensions and modifications), both at the qualitative and quantitative level. In particular sizable differences are found for the friction acting on the bubble wall.


Introduction
First order phase transitions (PhTs) in the early Universe proceed through the nucleation of bubbles of a stable phase within a metastable background. Afterwards bubbles expand in the hot plasma and coalesce, filling the whole space. This sequence of processes is characterised by a huge amount of energy stored in the gradients of the scalar field controlling the transition, in sound waves and turbulence in the plasma, all of them sourcing a stochastic background of gravitational waves.
The recent observation of gravitational waves has renewed a vivid interest in the study of the dynamics of such transitions. Indeed, the sensitivity regions of future experiments, such as the European interferometer LISA [1,2], the Japanese project DECIGO [3,4] and the Chinese Taiji [5,6] and TianQin [7] proposals, will probe a range of the expected peak frequencies of PhTs at the electroweak (EW) scale. These interferometers will provide us with a new tool that can support collider experiments in the quest for the physics beyond the Standard Model (BSM), in particular for theories potentially affecting the dynamics of the EW symmetry breaking.
Furthermore, a stochastic gravitational wave background is not the only cosmological relic left after the completion of a PhT. A matter-antimatter asymmetry, dark matter remnants, primordial black holes, magnetic fields and other topological defects can also be produced. A quantitative determination of these quantities obviously requires an accurate modelling of the PhT dynamics. This is controlled, among few other parameters, by the propagation velocity of the bubble wall.
In the steady state regime, the speed of the wall is a result of the balance of the internal pressure, due to the potential difference between the two phases, and the external friction exerted by the plasma particles impinging on the wall. In fact, the motion of a bubble drives the plasma out of equilibrium inducing a backreaction that slows down its propagation. Despite its relevance and the huge amount of literature on the topic, see for instance  for an incomplete list, this is one of the parameters on which we have less theoretical control.
The first computation of the bubble speed in the SM can be found in the seminal work of Moore and Prokopec [8,9] where the authors explicitly determined the friction induced by the plasma on the wall from a microphysics calculation, namely by evaluating all the relevant interactions between the plasma particles and the bubble. This requires the determination of the deviations from the equilibrium distributions of the different species in the plasma through the solution of the corresponding Boltzmann equations.
The formalism introduced in refs. [8,9], that we will denote as the "old formalism", necessarily requires the use of an ansatz for the distribution functions. This is needed to parametrize their momentum dependence and, then, to compute the local collision integral of the Boltzmann equation. In ref. [9] the fluid approximation was employed assuming that the deviation from the equilibrium distributions is entirely described by only three perturbations for each species in the plasma: the chemical potential, the temperature and the velocity fluctuations. The perturbations are then extracted by taking moments of the Boltzmann equation with suitable weights. In practice, the integro-differential Boltzmann equation is converted into a much simpler system of ordinary differential equations. By construction, the fluid approximation is equivalent to a first order expansion in the momenta of the deviation from the equilibrium distribution functions.
A peculiar feature of the fluid approximation is that the Liouville operator of the Boltzmann equation develops a zero eigenvalue at the speed of sound c s and, for larger velocities, all the perturbations trail the source term. This implies that any non-equilibrium dynamics is suppressed in front of the bubble wall [12] with significant consequences especially for non-local EW baryogenesis which would result to be extremely inefficient for bubble walls faster than c s . This has been the common lore for many years. But recently in refs. [15,16] it has been argued that the singularity is only an artifact of the first order truncation in momenta and of the particular set of weights chosen to extract the perturbations. Indeed, different choices of weights can shift the position of the singularity, suggesting that the speed of sound should not be a critical value for the particle diffusion as described by the fluid equations. In ref. [16] this problem was overcome by introducing a "new formalism", as dubbed by the authors, which relies on a different parameterization of the non-equilibrium distributions (specifically for the velocity perturbation), different weights and a factorization ansatz [41]. As a result, the new formalism wipes off the discontinuity at the speed of sound while still providing, for small velocities, quantitatively similar results to the fluid approximation.
The same issue has also been recently revisited in ref. [30] for the computation of the baryon asymmetry and in ref. [31] for the computation of the friction on the bubble wall, two problems that share many similarities. In these works the fluid approximation has been generalized by including higher orders in the small momenta expansion and the absence of the singularity for the perturbations of the heavy species has been corroborated. Besides the issue of the singularity, large differences in the friction arise, with respect to the old formalism, when higher orders are included. This confirms that the use of the fluid approximation, other than being not fully justified, is not particularly reliable, neither qualitatively nor quantitatively. A major consequence is that EW-baryogenesis is indeed achievable for supersonic bubbles opening up the parameter space of many BSM models, in which the observed baryon asymmetry can be reproduced while enhancing, at the same time, the strength of the stochastic gravitational wave background.
Even though the absence of the singularity for the heavy massive species could already be inferred in ref. [9], the speed of sound in the plasma still played a peculiar role in the old formalism as it provides a peak in the integrated friction for v c s . 1 The same peak (possibly accompanied by others, one for each vanishing eigenvalue of the Liouville operator) remains even if higher orders in the momenta expansion are included. As we will clarify with our analysis, such behavior is absent from the actual solution, confirming that the speed of sound is not a critical threshold of the friction for massive species. 2 The approaches discussed above are clearly affected by ambiguities. First of all, they all rely on an ansatz for the shape of non-equilibrium distribution functions which, both in the new and old formalisms (extended or not), is unavoidable in order to compute the collision integrals. Moreover, the choice of the basis and of the weights is not unique and different ansatzes have important qualitative and quantitative impacts on the resulting distribution functions. As such, a full solution of the Boltzmann equation that does not impose any specific momentum dependence is necessary to provide reliable quantitative predictions for both the non-equilibrium distribution functions and the friction exerted on the bubble wall, and to clarify the issue of the presence of a singularity. In fact, by feeding these new results into the equation of motion of the Higgs field, one will be able to carry out a precise computation of the wall speed and of the actual profile of the domain wall (DW). This is a necessary step towards a quantitative and reliable method to asses the potential of a given BSM extension to yield interesting predictions for the relics mentioned above. This is the goal of the present work.
In this paper we will present, for the first time, a fully quantitative solution of the Boltzmann equation. Since the absence of an ansatz prevents the direct computation of the collision integral, in order to extract the solution we will adopt an iterative method. In particular, as it will be detailed below, the collision integral can be split into two parts, one proportional to the solution itself and another effectively treated as a source term. At each iteration, the latter is evaluated using the solution obtained at the previous step. With a clever choice of the starting solution, convergence can be reached within a very small number of steps.
For the purpose of presenting the methodology and to quantitatively asses the differences among the aforementioned formalisms, we will consider the EWPhT and we will focus on the study of top quark species, the one with the strongest coupling, among the SM particles, to the Higgs profile and, as such, the one that provides the largest contribution to the friction. We leave for a future work the inclusion of the electroweak gauge bosons and of the background species.
The paper is organized as follows: in section 2 we describe our method while in section 3 we present and discuss the numerical results. In section 4 we give our conclusions and discuss future directions. All the technicalities related to the computation of the collision integrals are discussed in Appendix A.

The Boltzmann equation
As discussed in the Introduction, our goal is to determine the solution of the Botzmann equation for the distribution function of the plasma in the presence of an expanding bubble of true vacuum. Once the bubble reaches a radius much larger than the thickness of its wall, to a good approximation we can adopt the planar limit, considering a flat DW with a velocity parallel to its normal vector.
Assuming that (for long enough time) a steady state is reached, it is convenient to write the Boltzmann equation in the wall frame (i.e. the frame in which the DW is at rest) in which the solutions are stationary. Orienting the z axis along the velocity of the DW, the equation for the distribution function f of a particle species in the plasma is where m(z) is the mass of the particle, which in general depends on the position z, and (m 2 ) ≡ dm 2 /dz. The term C appearing in the right hand side of the equation is the collision integral, describing local microscopic interactions among the plasma particles, while L is the Liouville operator.
The collision term ensures that far from the DW, where the forces acting on the system basically vanish, each particle species approaches local thermal equilibrium. In the presence of a background fluid with a large number of degrees of freedom (in our case given by gluons and light quarks, which are not much affected by the Higgs phase transition), we can assume that the local thermal equilibrium is described by the standard Fermi or Bose-Einstein distributions for a fluid moving with velocity v along the z axis, 3 namely with β = 1/T and γ = 1/ √ 1 − v 2 . Deviations with respect to the local equilibrium distribution are present mostly close to the DW and are expected to vanish for z → ±∞. For small perturbations, the distribution function can be written as f = f v + δf and the Boltzmann equation can be linearized in δf : where we defined and C[δf ] denotes the collision integral linearized in δf . Notice that the only source term in the linearized Boltzmann equation comes from the Liouville operator L applied to the local equilibrium distribution. The collision integral, on the contrary, vanishes when computed on Sizable values for the source term are therefore present only close to the DW, where the non-trivial Higgs profile generates a non-negligible z dependence in m(z). Away from the DW, the Higgs profile is instead almost constant, thus giving (m 2 ) 0 and suppressing the source term. This behavior is in agreement with the naive expectation that deviations from local thermal equilibrium are only present close to the DW and should decrease to zero away from it.

Flow paths and the Liouville operator
As a first step towards finding a solution of the Boltzmann equation, we need to rewrite the Liouville differential operator in a simpler form. 4 It is straightforward to check that, along the paths on which both the transverse momentum 5 p ⊥ and the quantity p 2 z + m 2 (z) are constant, the differential operator simply reduces to a total derivative with respect to z: The physical interpretation of the paths is quite intuitive. They correspond to the trajectories of the particles in the (p ⊥ , p z , z) phase space in the collisionless limit. In this limit the energy of the particles and their momentum parallel to the DW are conserved (due to time invariance and translation invariance along the DW), therefore the trajectories of the particles are given by The condition p 2 z + m 2 (z) = const gives rise to different classes of flow paths. Since the mass of the particle species receives a contribution from the Higgs VEV, we expect it to smoothly increase going from the symmetric phase outside the bubble to the symmetry-broken one inside it. In particular, if we are interested in particles whose mass comes entirely from EW symmetry breaking (as it happens for the top quark and for the W and Z bosons), we can assume that m(z) → 0 for z → −∞, while it approaches a constant value m(z) → m 0 > 0 for z → +∞. In this case three types of flow paths are present: iii) for p z (−∞) ≤ −m 0 the path goes from z = +∞ to z = −∞ and has always p z < 0.
The three classes of curves are shown schematically in fig. 1 for the choice m(z) ∝ 1+tanh(z/L), with L denoting the wall thickness. The paths of type i, ii and iii correspond to the red, green and purple curves respectively. 6 Exploiting the flow paths we can straightforwardly solve any differential equation of the form where Q and S are generic functions of E, p z and z, and the factors 1/E and p z /E have been chosen for convenience. Rewriting the above equation along the flow paths we find whose general solution is where W is given by and all the integrals are evaluated along the flow paths. Notice that the lower integration boundary in the definition of W can be freely chosen (for each flow path) without affecting the result in eq. (9).
The function B(p ⊥ , p 2 z + m(z) 2 ), which is constant along the flow paths, is arbitrary and can be fixed by enforcing the required boundary conditions. Let us focus separately on the three classes of flow paths.
i) The first type of paths describes particles that travel in the positive z direction, and eventually enter into the bubble. It is natural to choose the boundary conditions in such way that δf vanishes at z → −∞, that is well before the particle hits the DW. This can be enforced by choosing ii) The second type of paths describes particles that initially travel in the positive z direction, hit the DW and are reflected. It is natural to choose the boundary conditions similarly to what we did for the previous type of paths. Therefore we have where the up arrow in the lower integration boundary indicates that the integration is performed starting from z → −∞ in the half path with p z > 0.
iii) The third type of paths describes particles that travel in the negative z direction, and eventually exit from the bubble. We can choose the boundary conditions in such way that δf vanishes at z → +∞, that is well before the particles exit from the bubble. This can be obtained by choosing The consistency of all these solutions requires Q < 0. We verified numerically that this condition is satisfied for the equations we are considering. The form of the solution clearly shows the role of the term (Q/E)δf in driving the system towards the local thermal equilibrium, i.e. in decreasing the value of δf . In fact, due to the exponential factors, the impact of the source term S is exponentially suppressed with the distance. The decay length is of order ∼ p z /Q and, as expected, decreases for larger values of the collision term.

Finding a solution for the Boltzmann equation
Although the full Boltzmann equation is not of the form of eq. (7), we can use the latter to implement an approximation by steps. The basic idea is to split the collision integral C[δf ] in two pieces: a term analogous to (Q/E)δf in eq. (7), and a second term that is included in the source term S and is used to correct the solution through iterations.
Let us now analyze in details the collision integral. For simplicity we consider the collision term for the 2 → 2 processes of a single particle species, but the general case can be treated in an analogous way. The collision integral is given by with where the sum is performed over all the relevant scattering processes, whose squared scattering amplitude is |M i | 2 . In the above formula N p is the number of degrees of freedom of the incoming particle with momentum p, k is the momentum of the second incoming particle, while p and k are the momenta of the outgoing particles. The ± signs are + for bosons and − for fermions. From the above expression we can easily derive the collision integral for the linearized Boltzmann equation. As a consequence of the conservation of the total 4-momentum in the collision processes we have that for the local equilibrium distribution C[f v ] = 0. Moreover the linear terms in δf can be expressed as where the −(+) sign in the sum applies to incoming (outgoing) particles.
The C[δf ] collision integral can therefore be split in two parts. One of them depends only on δf (p) and is given by (17) This expression is clearly analogous to the term (Q/E)δf in eq. (7). The second part of the collision integral includes the terms in which δf depends on k, p or k and thus appears under the integral sign. We collectively denote these terms by δf . The numerical determination of the various contributions to the collision integral can be drastically simplified through a clever choice of integration variables. The explicit procedure is explained in Appendix A.
In order to numerically solve the Boltzmann equation, a possible strategy is to formally rewrite it as eq. (7) by including δf in the source term S. The solution can then be found by iteration, inserting into the equation the value of δf obtained by using the solution at the previous step.

Numerical analysis
In this section we apply the iterative method explained above to numerically solve the Boltzmann equation. For simplicity we focus on a single species in the plasma, the top quark, which is the state with largest coupling to the Higgs and is thus expected to provide one of the most relevant effects controlling the DW dynamics. The analysis of the top quark distribution should be sufficient to provide a robust assessment of the plasma dynamics and to obtain an indication of how much the weighted method used in the literature to solve the Boltzmann equation is qualitatively and quantitatively accurate. We leave for future work the inclusion of the contributions from the W and Z bosons, which are expected to be roughly of the same size as the top quark ones.
The iterative approach explained in the previous section could be straightforwardly applied to determine the solution of the Boltzmann equation. However, in order to improve the convergence of the iterative steps, a slight modified procedure proves more convenient. Since the separation of the collision integral into a contribution to (Q/E)δf and a contribution to S is to a large extent arbitrary, we can devise a splitting that helps in reducing as much as possible the source term.

Annihilation only
Focusing on the top quark case, it can be shown that the main contribution to the collision integral comes from the annihilation process tt → gg, whereas the scattering of tops on gluons and light quarks gives smaller contributions. In our numerical analysis we will therefore consider at first only the contribution from annihilation, including scattering effects afterwards.
In the annihilation case, the linear terms in the perturbation δf appear in the following combination (see eq. (16)) where f g v denotes the equilibrium distribution for the background gluons (which is approximately unperturbed since the number of degrees of freedom in the background species is large). In the above formula the δf (p) and δf (k) terms play an analogous role, but their effects become distinct in the collision integral since an integration over k is performed. It is nevertheless evident that the impact of the δf (k) term in the Boltzmann equation is not particularly suppressed, as can be understood averaging the equation by integration over p, in which case the δf (p) and δf (k) become exactly equal. This line of reasoning suggests that treating the δf (k) contribution as source term, while including the δf (p) term in (Q/p z )δf might lead to slow convergence. To overcome this difficulty we will use a slightly modified procedure. We rewrite P as and we interpret the first contribution as (Q/p z )δf , while the second one (the one in round parentheses) is treated as a source. In this way the contribution to the source is partially canceled and faster convergence is achieved. 7 To determine the numerical solution of the Boltzmann equation we used a dedicated C++ code, validating the results with Mathematica [45]. The solution was computed on a threedimensional grid in the variables z, p ⊥ and p z restricted to the intervals z/L ∈ [−7, 7], 8 p ⊥ /T ∈ [0, 15], and p z /T ∈ [−15, +15]. The solution was computed on a grid with 50 × 300 × 100 points, which was further refined in the region p ⊥ /T < 1 and |p z | ∼ m 0 , where the solution showed a fast-varying behavior. Convergence of the solution (at the ∼ 0.1% level) was achieved within three iterative steps for all values of the wall velocity. We modeled the bubble wall assuming that the Higgs profile has the following functional dependence on z [46]: where L = 5/T is the thickness of the bubble wall and φ 0 = 150 GeV is the Higgs VEV in the broken phase. We fixed the phase transition temperature to T = 100 GeV. This choice of parameters, as we will see, determines the presence of friction peaks in the old formalism solution. It is thus well suited for differentiating the various formalisms and highlights the differences among them. An important quantity that can be derived from the numerical solution is the friction acting on the domain wall, which corresponds to the expression [9]  where N denotes the number of degrees of freedom (N = 12 for the top/antitop quark system).
In the left panel of fig. 2 we show the friction integrated over z as a function of the wall velocity (solid black line). The total friction shows a smooth behavior with a (nearly) linear growth as a function of the wall velocity.
In the same plot we compare our result with the ones obtained with the weighted methods. In particular the green lines correspond to the total friction computed in the old formalism (OF) of ref. [9], taking also into account higher-order terms in the fluid approximation [30]. The old formalism results at order 1, 2 and 3 are given by the dotted, dashed and solid lines respectively. The solid red line, instead, is obtained using the new formalism (NF) of ref. [16].
Our result for small and intermediate velocities, v 0.5 is in fair numerical agreement with the old formalism ones, which show a minor dependence on the order used for the computation. At higher velocities, instead, the old formalism develops some peaks related to the speed of sound in the plasma and to any other zero eigenvalue of the Liouville operator. The number of peaks and their shape crucially depend on the approximation order, denoting an intrinsic instability of the old formalism method. 9 Our results for the full solution of the Boltzmann equation show that the peaks are an artifact of the old formalism approach and that no strong effect is present in the top contributions for velocities close to the sound speed one.
On the other hand, the new formalism correctly predicts a smooth behavior for the total friction for all domain wall velocities. A roughly linear dependence on v is obtained up to v 0.8, while for larger values a faster growth is found, in contrast with the behavior of the full solution (FS) result. The quantitative agreement with the full solution is good only for very low velocities, v 0.1, while order 50% differences can be seen for higher velocities.
For a more refined comparison of the results we show in fig. 3 the behavior of the friction F (z) as a function of the position. The plots clearly show that the overall shape of the friction is very similar in all approaches, the main difference being the height of the peak. This property is not unexpected, since the size of the perturbation δf is controlled by the source term in the Boltzmann equation, whose z dependence is given by dm 2 /dz. One can easily check that the shape of all the curves in the plots roughly agree with the function dφ 2 (z)/dz.
A more detailed comparison of δf as a function of z, p ⊥ and p z shows drastic differences among all the approaches. Although the overall size of δf is comparable in all formalisms (being controlled by the source term), the various solutions significantly differ even at the qualitative level in most of the kinematic regions. We conclude from this comparison that the fluid approximation is not reliable if we include in the Boltzmann equation only the top annihilation channel. We will see in the following that, introducing the top scattering processes, a better agreement is found.

Full solution
We now consider the Boltzmann equation for the top quark distribution, including in the collision term also the main top scattering processes, namely the ones onto gluons tg → tg and onto light quarks tq → tq. We found convenient to include these additional contributions treating them as source terms in the iterative steps.
To determine the numerical solution we used a grid analogous to the one described in the annihilation-only case. The convergence of the iterative procedure is somewhat slower when top scattering processes are taken into account. For v ≤ 0.6 we used the solution of the annihilationonly case as starting ansatz and we performed six iterative steps to reach a good convergence. For higher velocities the annihilation-only solution is not a convenient choice for the first iterative step, thus we started from the full solution determined for a lower value of v. Also in this case six iterations were sufficient to achieve convergence.
We found that the scattering processes significantly modify the solution of the Boltzmann equation, especially for large values of the domain wall velocity (v 0.5). The impact on the total friction acting on the domain wall is shown in the right panel of fig. 2. Analogously to the annihilation-only case, an almost linear dependence on the wall velocity is present for small and intermediate v values, but a flattening is present at higher velocities. Quantitatively, the scattering processes induce only minor corrections to the total friction for v 0.6, while a decrease of order 25% is found for v ∼ 0.8.
The impact of the scattering processes on the solution obtained through the weighted methods is, on the contrary, much more pronounced. The old formalism approach (green lines in fig. 2) including only the lowest-order perturbations predicts a strong peak for v 0.55. The peak however gets substantially reduced once higher-order perturbations are included in the expansion, with a milder additional peak forming for v 0.75. We expect that including additional higher-order perturbations could smoothen the curve, giving a qualitative behavior similar to the one we get with the full solution, with a linear behavior for v 0.6. At the quantitative level, however, the old formalism solution differs from the one we found by order 10 − 25%.
The result obtained through the new formalism (red line in fig. 2) is also substantially modified by the scattering contributions. In particular the increase in the friction for v 0.8 is removed and a maximum followed by a mild decrease is now found for v 0.6. The new formalism prediction is now in good quantitative agreement with our result for v 0.2, while differences up to order 50% are found for larger velocities.
The friction as a function of the position z for some benchmark wall velocities is shown in fig. 4. Analogously to what we found for the total friction, the results we obtain with our method are only mildly modified by the scattering contributions. In particular the shape of the friction remains almost unchanged with only minor modifications in the overall normalization. Similar considerations apply for the shape of the z dependence of the friction in the old and new formalism. In this case, however, significant changes in the overall normalization are found, as expected from the above discussion on the total friction. Finally we show in fig. 5 the perturbation δf for the benchmark velocity v = 0.2. The results for different velocities are qualitatively analogous, the main difference being an overall rescaling with a limited change in shape. The plots show the full solution of the Boltzmann equation we got in our analysis, along with the results obtained applying the old and new weighted approaches. Notice that the new formalism does not fully determine the velocity perturbation, whose impact can only be computed averaging over the momentum through a factorization ansatz [41]. To plot the solution in the new formalism we chose to identify the distribution perturbation with following eq. (B5) of ref. [16]. This identification tends to produce a divergent behavior for small p z , which however has no impact on the determination of the friction since it is odd in p z . In the first (second) row of the figure we show the dependence of the even (odd) part of δf on the momentum along the z axis, p z . The plots are obtained fixing p ⊥ = 0, but similar results are found for p ⊥ /T 2 (for larger p ⊥ the solution is significantly suppressed and its impact on the domain wall dynamics is subleading). The plots in the first row show that, at the qualitative level, the old and new formalisms fairly reproduce the overall shape of the even part of the solution, although a somewhat different behavior is found for small p z (|p z /T | 0.5). This difference is most probably due to the fact that in the weighted approach the term is neglected. This term, although subleading in most of the kinematic space, dominates close to p z = 0, where the (p z /E)∂ z term vanishes. In spite of the fair qualitative agreement, large quantitative differences are present between the full solution and the ones obtained with the weighted approaches.
The agreement of the various formalisms in the determination of the odd part of δf proves quite poor. In particular, marked differences are found for z 0, in which case the highp z behavior of the solution is not captured by the old formalism, even including higher-order corrections. A mildly better agreement is found for z/L 1. The new formalism tends to reproduce the correct shape for p z /T 1, but presents large differences for small p z . It must be noticed that the odd part of δf does not contribute to the friction, thus the large differences found among the various solutions do not show up in the determination of F (z).
On the third row of fig. 5 we show δf as a function of z for p ⊥ = 0 and for the benchmark values p z /T = −1, −0.1, 0.1, 1. The old formalism reproduces the overall qualitative behavior of the solution. Higher-order terms tend to improve the agreement, although failing to fully reproduce the full solution, especially at the quantitative level. The new formalism is in qualitative agreement with the full solution for |p z /T | 1, while completely fails to reproduce the correct shape for small p z .

Conclusions and Outlook
In this paper we presented for the first time the fully quantitative solution of the Boltzmann equation that describes particle diffusion in the presence of a moving domain wall. Contrary to the existing approaches, we did not rely on any ansatz nor we imposed any momentum dependence on the non-equlibrium distribution functions. This clearly represents a necessary step towards a reliable understanding of the bubble wall dynamics. Using the friction obtained with the numerical method developed in this work, one can solve the equation of motion of the Higgs profile (or of any other scalar field driving a first order PhT) and extract the velocity of the domain wall as well as the features of its shape, such as the wall thickness. These parameters crucially impact on the prospects of any BSM theory to predict interesting cosmological signals, such as a gravitational wave background and the amount of matter-antimatter asymmetry.
We critically compared our results with the ones obtained using the formalisms developed so far in the literature, namely the fluid approximation originally developed in ref. [9], its extended version [30,31] (we dubbed both approaches 'old formalism', following ref. [16]), and the "new formalism" [16].
To establish our approach, we focused on a slightly simplified set-up in which only the top quark contribution to the DW dynamics is taken into account. Other species can however be included in a straightforward way. We computed numerically the distribution function for the top species and we obtained the friction F that the plasma exerts on the DW. The latter quantity is shown in fig. 4 as a function of the position z for the three different setups, namely: the old formalism (OF), the new one (NF) and our full solution (FS). The spatial dependence of the friction is quite similar in the three cases because the overall shape is mainly determined by the source term in the Boltzmann equation, namely dm 2 (z)/dz. There is however a significant disagreement at the quantitative level, as can be seen in the right panel of fig. 2, where the integrated friction is plotted as a function of the DW speed v. For small velocities v 0.2 a good agreement between the new formalism prediction and the full solution is found, while the old formalism (both original and extended) shows minor differences, of order 10%. The agreement significantly worsens at larger velocities. In this case the new formalism predicts a significantly smaller total friction, reaching a maximum much earlier than the full solution. The old formalism, on the other hand shows a different qualitative behavior, with a series of peaks related to the zero eigenvalues of the Liouville operator. 10 These features, which strongly depend on the order at which one fluid approximation is truncated, do not seem physical and are not present in the full solution, which shows a completely smooth shape (linear behavior for v 0.7 and a flattening for larger velocities).
In fig. 5 we also compare the distribution perturbation δf for the various approaches. Although in some kinematic regions a qualitative agreement can be seen, the differences among all approaches are quite strong. In particular the new formalism shows large differences in the odd part (with respect to p z ) of the perturbation. This difference does not show up in the friction result, since only the even part contributes to F (z).
We add that, as an intermediate step in our procedure, we considered the set-up in which only the annihilation channel for the top quarks is included in the collision integral, excluding the scattering processes. The friction for the full solution proves remarkably similar to the one in the complete set-up (see the left panel of fig. 2), apart from the fact that the maximum is reached for larger DW velocities. The old and new formalisms, on the other hand, show drastically different behavior. In particular the old formalism predicts sharp peaks connected to the sound speed. The inclusion of higher orders in the fluid approximation does not seem to achieve convergence in a reliable way.
As we mentioned, for the purpose of presenting the methodology and setting up the stage for a determination of the velocity of the bubble wall, in the present paper we only considered the top quark contribution to the DW dynamics. The inclusion of the electroweak gauge bosons and of the background species is clearly important to obtain quantitatively reliable predictions. We leave the investigation of this aspect for future work.
We also exploited another minor simplification in the computation of the collision integrals, ignoring the space dependence of the collisional kernels (see Appendix A), which appears through the top mass in the integrated equilibrium distribution functions. This approximation is also used in the old and new formalisms, and is expected to induce only minor corrections to the results. Within our approach the full space dependence could be taken into account, at the cost of increasing the computation time.
process Table 1: Amplitudes for the scattering processes relevant for the top quark in the leading log approximation. In the tq → tq process we summed over all massless quarks and antiquarks.

A Evaluation of the collision integrals
A.1 The term proportional to δf (p) We focus, at first, on the term of the collisional integral proportional to δf (p), which, for a single matrix element, reads (24) To evaluate the integral it is convenient to change variables through a boost, going to the plasma frame, in which the Boltzmann distribution is the standard equilibrium one f v . We denote by a bar the momenta in the plasma frame, namelȳ and analogously for k, p and k . We thus get (notice that the integration measure d 3 p/E p is invariant under boost) In order to evaluate the integrals, we follow the approach of ref. [9], including only leading log contributions. In this approximation we can also neglect the masses of the particles involved in the scattering. This approximation significantly simplifies the numerical evaluation, since it removes any explicit dependence on the z coordinate in the integrals. Closer inspection of the integral appearing in eq. (26) shows that it is invariant under rotation of the three-momentum components ofp, thus it is just a function of Ep. 11 The evaluation of the integral can be simplified by exploiting the delta function and the symmetries of the integrand. In this way one can perform analytically five of the nine integrals. An efficient parametrization for performing the integration is presented in ref. [47].
In the leading log approximation, only t-channel and u-channel scattering amplitudes are relevant (see table 1). So we can focus on these two types of contributions and neglect s-channel processes (and interference terms).

t-channel parametrization
We start by considering amplitudes coming from t-channel diagrams. The integration over d 3k can be easily performed exploiting the δ-function. The remaining integrals can be handled through a change of variables. As in ref. [47], we introduce the three-momentum q ≡p −p = k −k . Rotational invariance allows us to trivially integrate on the orientation of q. Fixing q to be along a z axis, we can express the orientation of thep andk momenta in terms of the polar angles θp q and θk q and the azimuthal angle φ between thep-q and thek-q plane.
The remaining delta function can be handled by introducing an additional variable ω linked to the t Mandelstam variable as t ≡ ω 2 − q 2 , where q ≡ |q|. In this way the integrations on the angles θk q and θp q can be performed analytically and one is left with the final expression for the integral on the second line of eq. (26): As alternative parametrization, which can help in the numerical integration and in studying the behavior of the integral, one can define in terms of which The χ ± parametrization can be also useful to leave as last integration the one on Ek: This choice of integration order clearly shows the symmetric role of Ep and Ek in the collisional integral.
The expressions for the s and u Mandelstam variables as a function of ω, q and Ek are given by while the relative angles between the three-momenta are given in ref. [47] (see Appendix A.2, eqs. (A21a)-(A21e)), among which , .

u-channel parametrization
Analogous formulae can be found for the u-channel parametrization, by exchangingp andk in the t-channel parametrization. In this way the integral becomes with q ≡k −p =k −p and The expressions for the s and u Mandelstam variables are given by

Structure of the contribution
From the above formulae we can easily infer the global structure of the collisional term proportional to δf (p). The quantity (we consider the t-channel parametrization for definiteness) only depends on Ep, as we already anticipated. Therefore we get We can now go back to the wall frame, obtaining Notice that the massless-limit approximation introduced a small 'mismatch' in this expression, since we chose f v (p) in the prefactors to have the full mass dependence (from the definition of E p ). In the approach to the solution via the use of weights, instead, f v is treated in the massless limit for all the factors in the collisional integrals. This problem could be solved by also considering the massive form for all the f v factors inside the collisional integral. This however is computationally more demanding, since it introduces an explicit z dependence in the integrand, so that the kernel should be evaluated also as a function of z. 12 The numerical analysis shows a behavior K(Ep) ∼ log Ep + const (44) which, as expected, has a logarithmic divergence for mass and thermal mass going to zero. We can thus infer the rough behavior (at least for small v) plasma reference frame and to the angle θpk between the momentap andk. Putting everything together we find which can be rewritten as 2π 0 dφk f 0 (k) δf (k ⊥ , γ(k z + vEk), z) (−f 0 (k)) .
(52) The collision integral for the scattering processes includes an additional set contributions in which δf (p ) or δf (k ) appears. In analogy to the previous case, for the δf (p ) terms, we can first perform the integrals over k and k , obtaining the following expression which can also be rewritten as The contributions from δf (k ) can be treated in an analogous way.

A.2.1 Evaluation of the K 1 kernel
The evaluation of the kernel K 1 can be performed as in ref. [48]. As a first step we perform the integration over k exploiting the Dirac delta: Notice that in the above expression we expressed the energies E p and E k in the Lorentzinvariant form u µ p µ and u µ k µ . As we will see, this is useful to keep track of the changes of reference frame. As a second step, we rewrite the Dirac delta (in the center-of-mass (COM) frame) as with s = (p + k) 2 the usual Mandelstam variable. The integration over p can be performed by rewriting d 3 p = E 2 p dE p d cos θ dφ, where θ is the angle between p and p in the COM frame of the scattering process: As a last step we need to compute u µ p µ and u µ k µ in the COM frame. We conveniently choose the orientation of the COM frame axes such that u y = 0 leading to We then introduce the four-vectors P µ and Q µ defined as In the COM frame we find that We can get a further simplification by choosing the frame such that Q lies along the z axis. Since, in the massless case, |P| = |Q| = √ s, the vectors P µ / √ s and Q µ / √ s coincide with the versors along the first and fourth Minkowski directions.
The u 0 and u z components can be easily computed in terms of the momenta of the particles in the plasma frame : The u x component can be determined from the condition u µ u µ = 1: Putting everything together we find that u µ p µ is given by Similarly we find u µ k µ = 1 2 Ep(1 − cos θ) + Ek(1 + cos θ) + 2EpEk(1 + cos θpk) sin θ cos φ .

A.2.2 Evaluation of the K 2 kernel
We now discuss the evaluation of the K 2 kernel: Also in this case we follow ref. [48]. We introduce the four-vectors Finally, from u µ u µ = 1, one gets In order to make the numerical evaluation of the kernel more stable, we used the following coordinate change |K| = √ −t tan θ, and then we defined 1/ cos θ = x. The expression for K 2 becomes with u µ k µ = 1 2 (Ep + Ep )x − x 2 − 1 2EpEp (1 + cos θpp ) cos φ + (Ep − Ep ) ,