Friction pressure on relativistic bubble walls

During a cosmological first-order phase transition, particles of the plasma crossing the bubble walls can radiate a gauge boson. The resulting pressure cannot be computed perturbatively for large coupling constant and/or large supercooling. We resum the real and virtual emissions at all leading-log orders, both analytically and numerically using a Monte-Carlo simulation. We find that radiated bosons are dominantly soft and that the resulting retarding pressure on relativistic bubble walls is linear both in the Lorentz boost and in the order parameter, up to a log. We further quantitatively discuss IR cut-offs, wall thickness effects, the impact of various approximations entering the calculation, and comment on the fate of radiated bosons that are reflected.

In this work we make progress on the computation of the friction pressure on relativistic bubble walls. In Sec. 2, as a warm-up we compute the retarding pressure coming from particles acquiring a mass in the broken phase, as initially derived in BM '09 [56], and recover their results. In Sec. 3, we review the possibility for the incoming particle to radiate one gauge boson acquiring a mass m in the broken phase, and again recover the result from BM '17 [57], P NLO ∝ g 2 m γ T 3 nuc , where T nuc is the bubble nucleation temperature, g is the gauge coupling constant and µ is the IR cut-off on the transverse momentum of the emitted boson. In Sec. 3.4, we discuss various possibilities to cure the IR logarithmic divergence. In Tab. 1, we show that for large supercooling m/T nuc 1 and/or large gauge coupling constant α, perturbativity breaks down and we must account for the possibility to radiate multiple vector bosons.
In Sec. 4, we perform an analytical Sudakov resummation of virtual and real emissions at all leading-log orders (LL). The perturbative splitting probability can then be used to express the mean exchanged momentum, and thus the pressure, as a resummed quantity that includes leading-log real and virtual corrections to all orders, see Eq. (74). In the limit where the initial energy is large, E a /m 1, such that the kinematics of the multiple emissions can be considered as independent, we find that the pressure is linear in both the wall Lorentz factor γ and in the order parameter of the transition m, up to a log. In Sec. 5, we simulate a particle shower using a Monte-Carlo algorithm, and we confirm the analytical resummation in the limit E a /m 1, in which energy depletion due to multiple emissions can be neglected. In Sec. 6, we deduce the bubble wall Lorentz factor at collision time and discuss the consequences on the nature of the GW source. We conclude in Sec. 7.
To keep our paper easier to read, we defer a lot of technical details and complementary -yet interesting -calculations to a series of appendices. In App. A we compute the vertex function for boson emission, in App. B we compute the particle mode functions precisely, and discuss the validity of the thin-wall and relativistic-soft-collinear limits, in App. C we compute the 3 → 2 boson scattering rate, in App. D we sketch the effect of particles reflected multiple times, in App. E we treat the case of a massless vector boson, in App F we discuss the discrepancy of our result with the one of [1] that found a scaling P LL ∝ γ 2 T 2 nuc , in App. G we compute the backreaction of successive boson emissons on the kinematics of the parent particle, and in App. H we compute the Lorentz factor of a constantly accelerating bubble wall.

LO friction pressure
In this work, we assume that the wall is moving at relativistic velocities γ 1, such that we can work in the so-called ballistic regime in which we can neglect the interaction between neighboring particles during the time when they cross the wall, see e.g. [54]. The leading-order pressure comes from the particle getting a mass in the broken phase [56].
where g a is the number of internal degrees of freedom. The momentum change is given by where we assumed energy conservation and the relativistic limit E ∆m. We compute ∞ 0 4πp 2 dp (2π) 3 1 e p/Tnuc ± 1 ∆m 2 2p = ∆m 2 T 2 nuc 24 × 1 boson 1 2 fermion (3) where T nuc is the nucleation temperature. Therefore where c a = 1 (1/2) for bosons (fermions). If the negative pressure due to the vacuum energy difference of the phase transition ∆V is larger than the LO retarding pressure in Eq. (2), then the bubble wall is supposed to be accelerated to larger and larger γ factors until either it collides with other walls or until soft particle radiation further contributes to the pressure and eventually stops the wall from accelerating. The rest of the paper is dedicated to the computation of the latter contribution. 3 The splitting probability at first order 3

.1 Transition splitting
In this section we discuss another contribution to the pressure, which arises when a particle entering the wall radiates a gauge boson which gets mass in the broken phase [57]. The resulting NLO pressure reads where dP a→bc is the differential splitting probability, p a and p b are the momenta of the incoming particle before and after the splitting while p c is the momentum of the radiated boson, 1 see Fig. 1. The momentum p a is the momentum in the far past and thus in the symmetric phase, while p b and p c are the momenta in the far future and thus in the broken phase if they are transmitted, or in the symmetric phase if they are reflected by the wall boundary. We summed over all the species a likely to participate in the process, g a being their number of degrees of freedom. The Pauli blocking or Bose enhancing factors 1 ± f b are 1, while 1 ± f c sum to 1 when considering both absorption and emission processes. 2

The momentum exchange
Kinematics. Upon introducing we can write 1 Note that our notation for 'a', 'b' and 'c' is different from [57] where the roles of 'b' and 'c' are interchanged. 2 The interaction Hamiltonian can be written as (see e.g. [87]) where ax are the usual creation operators in Fock space. Then the transition amplitudes for emission and absorption read, respectively where +/− refers to boson/fermion statistic. We deduce the interaction rate accounting for both emission and absorption Hence we see that as long as (fa − f b )fc fa, we have |M a→bc | 2 − |M bc→a | 2 |M0| 2 fa.
As we will see, for non-abelian theories we will cut fc ∼ 1/g 2 , so our calculation in those cases holds good as long as |fa − f b | g 2 . We leave a more quantitative study for future works. where m a (z), m b (z) and m c (z) are the masses of the three particles involved in the vertex as a function of z, and k ⊥ ≡ | k ⊥ |. As the wall breaks z-translation invariance, the momentum alonĝ z is not conserved. We approximate the wall by a Heaviside function at z = 0 where the masses with 's' and 'h' denote the ones at infinity in the symmetric and Higgs phase, respectively. As long as the masses in the symmetric phase are small compared to the ones in the broken phase, we can safely assume that m a,s = m b,s = 0.
We can not make the same simplification for m c because m c,s = 0 implies the existence of a double (soft and collinear) divergence, see Eq. (38). Hence, m c,s plays the role of an IR cut-off which regulates the double singularity, whose possible values are discussed in Sec. 3.4. Note that the parameterization of the kinematics in Eq. (14), (15) and (16), imposes and the corresponding allowed region for x and k ⊥ is shown in Fig. 2. Since all the results derived in this work are UV insensitive, we can simplify the upper boundaries in Eq. (19) and (20) as The associated correction terms are anyway beyond the soft-collinear approximation, that we will assume when deriving the phase of the mode function in Eq. (34) and the vertex function in Eq. (31).

The splitting probability
Relation between splitting probability and matrix element. The differential splitting probability is given by [57] dP a→bc ≡ We have regulated the behavior at infinity by introducing an imaginary momentum ±i , which we can safely set to 0 in the right-hand side of Eq. (35) since A h , A s = 0. In App. B.1.1, we derive the mode function χ c (z) in the presence of a finite wall thickness L w and we show that it can be neglected up to logarithmic corrections in the limit, cf. Eq. (137) Since we expect L w ∼ m −1 c,h , cf. Eq. (138), we conclude that the wall thickness can be neglected as soon as m c,s m c,h . Note that the validity of using the Heaviside function is not given by the comparison between the (inverse) wall thickness and the momentum of individual a, b, or c particles.
In App. B.1.2, we compute the full mode function valid for any p z c . Particularly, χ c (z) contains a second wave propagating in the opposite direction which implies that in the limit p z c m c,h , the particle c is reflected. The full matrix M accounting for reflection and transmission coefficients is given in Eq. (146).
Perturbative splitting probability. Therefore, the matrix element in Eq. (35) becomes The splitting probability in Eq. (29) reduces to where Π(k ⊥ ) contains the IR and UV suppression factors and Here note that we integrated the radial direction of k ⊥ as d 2 k ⊥ = πdk 2 ⊥ . The charge factors C abc in the SM can be found in [1]. In the rest of the paper we assume m c,h m c,s . In Eq. (38), we have replaced p z i with i = a, b, c in Eq. (29) by E i . In App. B.2, we investigate the validity of the relativistic-soft-collinear approximation for the phase space factor 1/p z i , the vertex function V , the phase A of the mode function and the momentum exchange ∆p and we show that this underestimates the final result by a few percents only for T nuc = 10 −2 T start . to usual splitting functions in collider context [88], vanishes in the limit where the symmetry is restored m c,h → m c,s or k 2 ⊥ m 2 c,h . In contrast, as first claimed by [58] the splitting function used by [1] does not go to zero in the limit where the distinction between the two phases disappears. We give more details on the origin of this discrepancy in App. F.

The IR cut-off
The splitting probability Eq. (38) is divergent in the IR, for k 2 ⊥ → 0. That divergence is regulated by values of m c,s > 0, via the factor k 2 ⊥ /(k 2 ⊥ + m 2 c,s ) 2 , or by some other physical process. To encompass both possibilities, we find it convenient to define a general IR cutoff µ and rewrite Eq. (38) as where by definition µ ≥ m c,s .
This definition gives us also the possibility to treat µ as a free parameter, so that one could later account for IR cut-offs which are not relevant for the physical situation of interest for this paper, or which we simply miss. Let us start by discussing the case µ = m c,s , for which Eq. (41) goes back to Eq. (38). Values of m c,s > 0 are guaranteed by screening effects in the plasma [89,90], with the sum running over all species in the plasma i that couple to c, f i being the occupation number of the i particle, g i the number of relativistic degrees of freedom of species i, C i the quadratic Casimir (g 2 C i = charge squared for abelian theories) and d A the dimension of the adjoint representation. We discuss two possible screening effects in the next two paragraph, and later discuss other possible cut-offs for k ⊥ .
Thermal mass. A contribution to f in Eq. (43) is always given by the particles in thermal equilibrium. Taking a U (1) and an SU (N ) gauge theories as examples, we obtain where N f is the number of Dirac fermion flavors in the fundamental representation. The thermal mass in Eq. (44) and (45), which for the sake of simplicity we write as constitutes the minimal IR cut-off for abelian and non-abelian gauge theories.
Phase-space saturation. As the occupation number of emitted vector bosons grows in the IR, it must exist a scale m sat below which perturbation theory breaks down and vector bosons start to act collectively [57]. 4 Due to the soft-collinear divergence of the splitting function, the occupation number f c (p c ) grows like where we have again used dp 0 c dp z c in the first equality and Eq. (38) in the second one, g a is the number of relativistic degrees of freedom of particle a and ν a = 1 (3/4) for bosons (fermions) assuming f a to be a thermal distribution. From Eq. (43) and (47), we can see that in the case of non-abelian gauge theories (without loss of generality, we focus on SU (N ) here and for numerical applications we fix N = 3), in the large γ limit, the emitted vector bosons back-react self-consistently on the screening mass which implies Note that m sat increases as γ 1/3 . Its value at the time of bubble collision depends on whether the wall runs away or reaches a terminal velocity. Upon plugging typical values of the parameters, we obtain 5 where we have introduced parameters that will be discussed later in the paper: T start is the temperature where vacuum domination starts (Eq. (65)), H * is the Hubble parameter at the time of the phase transition, β is the time variation of the nucleation rate (Eq. (246)), γ run is the Lorentz factor of bubble walls that run away evaluated at the time of collision (Eq. (105)), γ LL is that of bubble walls where the external pressure compensates the internal one ∆V (Eq. (104)), c w is the bubble radius at nucleation in unit of T −1 nuc (Eq. (244)), and κ is the exchanged momentum in unit of ζ a m c,h ln(m c,h /µ) (Eq. (80)).
The modification of the dispersion relation of vector bosons, that we just derived, relies on a perturbative particle description. However, as we now sketch, perturbativity breaks down for p c so low to give roughly the same occupation numbers f c that lead to our cut-off m sat above. Eq. (48) can be rewritten as where we define p * as some IR cut-off of the momentum of the c particle, because the integrand in Eq. (48) is peaked in the IR. where p * is the typical momentum of the c particle, which is 5 Note that if the screeening mass due to gluon collective behavior, msat in Eq. (52), becomes larger than the mass of free gluons in the broken phase m c,h , then the logarithmic divergence in the splitting probability in Eq. (38) is replaced by a 1/k 4 ⊥ divergence and the friction pressure becomes [57] msat m c,h =⇒ ∆p 0.3 ζa msat m c,h msat This change of scaling γ → γ 4/7 , which is enough to prevent bubble walls to run-away, should only occurs for large latent heat ∆V φ 4 , small gauge coupling α 0.01, and at the end of the bubble expansion stage in the friction dominated regime γ γLL.
identified as the IR cut-off because the integrand in Eq. (48) is peaked in the IR. Modes with p * m c,s are screened, hence the condition for the vector bosons to back-react on the dispersion relation can be recast as [57] g 2 f c (p * ) > 1.
The occupation number f c (p * ) of vector bosons is related to the vector boson wave function A µ by We deduce that the hierarchy between the 3 terms in the Lagrangian which is essential for perturbation theory to apply, breaks down as soon as we enter the regime of Eq. (54), or equivalently, as soon as p * < m sat in Eq. (52). Therefore, a treatment beyond perturbation theory would be needed to determine the friction pressure for momenta below m sat . We leave such an interesting study for future work, and here for simplicity we conservatively choose to still interpret m sat as an IR cut-off. Let us finally comment on abelian gauge theories. In those cases both previous arguments do not apply, because they rely on the existence of gauge self-interactions, so that one cannot use the cutoff m sat in Eq. (52). This motivates the use of the naive thermal cutoff in Eq. (46), with one potential limitation that we now comment upon. In the abelian case we expect an IR cut-off to arise from the presence of fermions and scalars originating from the splitting of the soft gauge bosons, themselves radiated from the incoming energetic quanta. This splitting of soft gauge bosons into softer fermions or scalar pairs is not enhanced, however one has many soft gauge bosons to start with, so that the IR cutoff induced by the produced fermions and scalars may perhaps be larger than the thermal mass. Determining this IR cut-off goes beyond the purpose of this paper, where we just content ourselves with pointing out this potential limitation of using the naive thermal mass of Eq. (46) in the abelian case.
Other cutoffs. Below we discuss processes that may regulate the IR divergence of Eq. (38) at values µ > m c,s . We anticipate that, in the physics case of our interest, we find that the thermal mass and the phase-space saturation discussed above provide IR cutoffs that are stronger than the effects we discuss next.
3-to-2 boson scattering. In the IR limit k 2 ⊥ → 0, the occupation number of emitted vector bosons diverges. At some point, in the case of non-abelian gauge theories, we could expect the occupation number of emitted vector bosons to be large enough to trigger 3-to-2 scatterings. In that case, the population of emitted vector bosons is depleted and the pressure stops growing. In App. C, we compute the resulting value of the IR cut-off µ 3→2 at tree level, and find that µ 3→2 < m sat for typical values of the parameters. We then consider m sat as the effective IR cut-off of our emissions. 6 Energy-momentum conservation. Successive boson emissions lower the energy and momentum of the parent particle (E b , p b ). The impossibility to radiate more energy and momentum than what is initially available, therefore, tames the IR divergence even when the cut-offs discussed so far go to zero. To describe this effect, we rely on a numerical Monte-Carlo simulation which we describe in Sec. 5 and App. G. We find that the number of emissions and the resulting ∆p saturate in the IR, see Fig. 16 of App. G. This leads to an effective IR cut-off which depends on the initial energy E a . However we find that the corresponding IR cut-off is smaller than the thermal mass and backreaction can be neglected.
Prescription for calculations. In the remaining part of the text, we consider the IR cut-off µ in Eq. (41) as a free parameter. For numerical applications, we will consider two benchmark scenarios according to whether the IR cut-off is set, in the case of abelian gauge theory, by the thermal mass µ = α 1/2 T nuc in Eq. (46) or, in the case of non-abelian gauge theory, by the screening length µ = m sat in Eq. (52) resulting from phase space saturation µ = m c,s m th = α 1/2 T nuc (abelian gauge theory), Max [m th , m sat ] (non-abelian gauge theory).

Perturbativity breakdown
Integrated perturbative splitting probability. We integrate x in Eq. (38) over the range defined in Eq. (21) with m c (z) → m c,s . We obtain the probability to emit a vector boson with transverse momentum k ⊥ in the collinear limit The last integration over k 2 ⊥ , Eq. (22), is logarithmically dependent on the IR cut-off µ. Upon identifying µ = m c,s and in the limit µ m c,h E a , the integrated vector boson emission probability becomes We recover the standard double logarithm for integrated probability of vector boson radiation (see e.g. [91][92][93][94]). The soft-collinear divergence is here regularized by the IR cut-off µ, which is not the energy threshold of some particle detector like in a collider context, because here all emitted vector bosons contribute to the observable ∆p and so are physical. We take the IR cut-off in our context from processes due to the existence of (high) particle densities, see Sec. 3.4 for more details.
Numerical estimates in physical scenarios. We now evaluate the emission probability P E , for different benchmark values in Table 1. The incoming energy is given by the thermal energy 3T nuc [95] boosted in the wall frame We assume that the bubble wall Lorentz factor γ is set by its value evaluated at the time of bubble collision, see Eq.
Without much loss of generality, we choose We assume that the vector boson mass in the broken phase is set to We choose to characterize the amount of supercooling by the ratio T nuc /T start where T start is the temperature when the universe becomes vacuum-dominated 7 Here N e ≡ ln (T start /T nuc ) gives the number of e-folds of inflation generated during the supercooled phase transition. As shown in Table 1, at large supercooling or large coupling constant, e.g. α 0.3 and T nuc 10 −3 T start , the perturbative calculation in Eq. (59) cannot be trusted and we must account for the possibility to radiate multiple vector bosons. 4 The splitting probability at all leading-log orders

Multiple boson emission
Exchanged momentum. For multiple c emission, the exchange momentum along z is where we have denoted by X and K ⊥ the sum of energies and transverse momenta of the n emitted vector bosons The values of the z momenta of the emitted vector bosons p z c i depend on whether they are transmitted or reflected by the wall 7 In the low-temperature expansion of the thermal potential, i.e. when the abundance of the particles i thermally correcting the potential of φ is Boltzmann-suppressed, mi(φ) Tstart, one has where Tc is the critical temperature when the two minima of the potential are degenerate, cf. App. A in [12], and where g * ,i and g * ,f are the number of relativistic degrees of freedom before and after the phase transition. with Note that Eq. (66) assumes that successive emissions take place in the same (xz) plane and Eq. (68) neglects the transverse recoils of the successive emissions on the momentum of the parent particle. In App. B.3, we show that those two approximations are very good.
Soft-collinear limit. To first order in the soft x i 1 and collinear k ⊥,i 1 limit, we obtain

Sudakov resummation
Poisson distribution. At leading-log order, the many-boson emission distribution follows a Poisson distribution, see e.g. Eq. (2.30) of [96] or Eq. (6.86) of [92], and the mean value of an observable O can be computed from where The product factors in Eq. (71) account for the n indistinguishable leading-log real emissions while the exponential resums the leading-log virtual corrections. The matching between the matrix element of the virtual correction in the argument of the exponential and the perturbative splitting probability stems from unitarity The real emissions are assumed to be independent and we neglect correlations which are higher order effects [96].
Mean exchange momentum. We can write where we have made use of Eq. (73). This resembles the average for a Poisson distribution n × λ n e −λ /n! = λ λ n−1 e −λ /(n − 1)! = λ. We conclude that in the soft-collinear limit in which ∆p is an additive observable, see Eq. (70), the average exchanged momentum in which both real and virtual leading-log emissions are resummed is identical to the naive expectation using the perturbative splitting probability as in the original article [57].
Note however that in the steps leading to Eq. (74), we have neglected the dependence of the kinematics boundaries on earlier emissions, i.e. we have kept fixed integration boundaries in Eq. (72). Instead, the energy and momentum of the parent particle should get depleted with the number of emissions. The proper average exchanged momentum ∆p taking into account these backreaction effects is given in Eq. (226) of App. G. Since this analytical formula is not easily tractable, we propose to include the effect of backreaction with a Monte-Carlo simulation in the next section, Sec. 5.

Analytical estimate
From plugging the fixed kinematic boundaries in Eq. (22) inside Eq. (72) and Eq. (74), in the limit µ m c,h E a and m c,s m c,h , we obtain where ∆p R and ∆p T are the contributions from reflected and transmitted vector bosons and where in Eqs. (76) and (78) we have used the splitting probability of Eq. (58), and in Eqs. (77) and (79) we have expanded the results of the integrals for µ m c,h E a . We conclude that the contribution from soft emitted vector bosons which are reflected at the wall boundary is at least comparable to the contribution from the transmitted ones (as anticipated in the fixed-order calculation in [58]). For parameters of physical interest, ∆p R never never reaches its asymptotic behaviour of Eq. (77), but stays slightly smaller, as we will study next. We thus write a ready-to-use approximated resummed analytical result, sum of the reflected and transmitted contributions, as In Fig. 3, we show the analytical estimate for κ evaluated with the IR cut-off µ, either set by the thermal mass α 1/2 T nuc , or by the screening mass m sat of phase-space-saturated boson bath, see Eq. (57). In Fig. 4, we confront the analytical treatment with the numerical one, based on a MC algorithm, which is discussed in the next section, Sec. 5. In App. B, we compute different corrections to Eq. (80) due to properly computing the mode function of particle c, and to not expanding the various square root functions leading to Eqs. (76) and (78). We summarize the list of corrections in Table 2 and Fig. 12. We also anticipate that Eq.
The IR cut-off µ is set to (brown lines) the thermal mass α 1/2 T nuc and to (green lines) the mass m sat from the collective behaviour of the phase-space-saturated emitted vector bosons (valid for non-abelian gauge theories only), see Eq. (57). The presence of the plateau can be understood by the fact that m sat does not depend on T nuc /T start for terminal-velocity walls, cf. Eq. (52). Inside the yellow region, we compute γ 10, cf. Eq. (104), such that interactions between neighboring incoming particles can not be safely neglected during the time of wall crossing and our analysis may break down.

Fate of reflected vector bosons
Previously we have discussed the possibility for radiated c particles to be reflected by the wall boundary whenever their would-be momentum in the hidden phase p c,h , defined by Eq. (69), is negative. We now discuss possible corrections to ∆p R in Eq. (75) due to the presence of those particles right in front of the wall.
We expect the momentum of the reflected c particles to change due to scatterings with the incoming a particles in the thermal bath. At a distance of the wall set by their mean free path in the plasma frame, see App. D where p c,z m c,h is their typical momentum, the c particles are expected to come back in direction of the wall. If their new momentum p c,z is larger than m c,h , they are transmitted to the symmetric phase. Otherwise, if p c,z < m c,h they are reflected again and so on. Those cycles of multiple reflections ended up by transmission are expected to give corrections to ∆p in Eq. (80).
The presence of this population of reflected c particles is expected to give corrections to the dispersion relation and to the computation of m sat in Eq. (52).
Due to successive scatterings with the reflected c particles, the incoming a particles are expected to lose a small fraction of their momentum E a . The exchanged momentum ∆p in Eq. (75) is independent of E a in the limit E a m c,h but goes to zero in the limit E a ∼ m c,h . The latter only arises in the region of extreme supercooling (e.g. T nuc /T start ∼ 10 −7 if φ = 1 TeV), see purple vertical line in Fig. 17. We expect the depletion of E a due to scatterings with reflected c particles, to shift the position of this purple vertical line to the right.
The successive scatterings of the incoming a particles with the reflected c particles, which act as a medium by itself, is expected to induce further splitting radiations, in addition to the one induced by the bubble wall. The presence of reflected particles then induces further splitting radiations, which in turn induce further reflected particles and so on.
Contenting ourselves with qualitative comments, we leave the quantitative study of those new effects for future works.

Monte Carlo simulation
In this section, we numerically generate the shower of emitted bosons with a Monte-Carlo (MC) simulation and compute the resulting momentum exchanged with the wall.

The Sudakov form factor
The survival probability. A necessary ingredient to realize a MC simulation is the probability P NE (k 1 , k 2 ) of non-emission in the interval k 1 ≤ k ⊥ ≤ k 2 [91][92][93][94]. The latter can be derived from the probability of non-emission in the infinitesimal interval [k 2 as follows. dP E is given in Eq. (58). We consider the finite interval [k 2 where P E (k 1 , k 2 ), is the perturbative probability to emit a vector boson with transverse mo- with which, in the limit where k 1 m c,h k 2 and m c,s m c,h , reduces to Resummation of leading-log real and virtual corrections. Another way to recover the survival probability P NE in Eq. (83) is The first exponential factor includes real emissions outside the interval [k 1 , k 2 ], while the second includes virtual emissions (inside the interval [µ, E a ]), both resummed to all leading-log orders.

The algorithm
Motivation. In order to give support to our analytical estimate in Sec. 4.3, but also to include backreaction, see App. G, we simulate the shower of emitted vector bosons with a Monte Carlo (MC) algorithm [91,96,97], which we now describe.
Recipe. Starting from the hard scale k ⊥,0 = E a , we generate the transverse momentum of the first boson k ⊥,1 by solving the equation  Figure 5: Trajectory in the phase space (x, k perp ) (also called Lund plane, e.g. [94]) for two given Monte-Carlo simulations. We set T nuc /T start = 10 −4 , α 1/10 and φ = 1 TeV, which imply m c,h /E a 6 × 10 −7 and γ 2 × 10 10 . Left: We set the IR cut-off to the thermal mass µ = α 1/2 T nuc , cf. Eq. (46), which leads to µ/m c,h 7 × 10 −6 . Right: We set the IR cut-off to the self-energy of phase-space-saturated boson bath µ = m sat , cf. Eq. (52), which leads to µ/m c,h 0.01. The orange points linked by arrows are the kinematics variables of the successive emitted bosons generated by the MC algorithm. The ones which pass the threshold E c > m c,h will be able to enter the broken phase (blue region) while the others will be reflected at the wall boundary (green region). In the red region, the emission energy is smaller than the thermal mass µ, and therefore emission is kinematically forbidden.
where P NE is the Sudakov factor in Eq. (83) and R is a random number between 0 and 1. The energy of the boson x 1 E a is a pseudo-random number generated from the perturbative splitting probability in Eq. (41), i.e. by solving The kinematics of the second emitted vector boson (k ⊥,2 , x 2 ) are determined the same way with k ⊥,0 replaced by k ⊥,1 , and so on. We stop the shower whenever the transverse momentum becomes smaller than the IR cut-off, see Eq. (41) In App. B.3, we present a MC algorithm which takes into account azimuthal emission angles and transverse recoils of successive boson emissions on the momentum of the parent particle.
Backreaction. The depletion of the energy-momentum of the parent particle as the emission continues, discussed in App. G, is taken into account after replacing the x upper boundaries of Eq. (88) by and the x upper boundaries of Eq. (89) by and by stopping the cascade whenever Results. Fig. 5 shows the phase space trajectory of one given MC simulation. The resulting momentum exchanged with the wall is given by the master formula in Eq. (66). In Fig. 4, we show that numerical computations based on MC shower and the analytical estimates of Sec. 4.3 agree up to percent level. 6 The bubble wall velocity 6.1 The final retarding pressure Non-confining PTs. The goal of this paper was to compute the retarding pressure at all leading-log orders (LL) for non-confining PTs, meaning PTs where particles simply acquire a mass in the broken phase. We have obtained where P LO is given in Eq. (4), which we rewrite here and P LL follows from Eqs. (6) and (80), In Eq. (97) g a is the number of relativistic degrees of freedom of particle a and ν a = 1 (3/4) for bosons (fermions). We also have boosted the phase space volume d 3 p a to the wall frame by introducing the Lorentz factor of the wall in the plasma frame γ. In Eq. (98) C abc is the charge factor for the SM and can be found in [1], and κ ≈ 4. In our numerical analysis we use a,b,c ν a g a C abc = 100 for simplicity. The level of approximation leading to κ ≈ 4 is studied in App. B. The value of the IR cut-off µ is discussed in Sec. 3.4. In the case of abelian gauge theory, we set it to the thermal mass µ = α 1/2 T nuc and, in the case of non-abelian gauge theory, we set it to the screening mass µ = m sat resulting from phase space saturation, see Eq. (57).
The regime of validity of Eq. (98) is µ m c,h E a and thin walls. Instead, outside those regimes, we should use Eq. (97)  We conclude that after having performed a leading-log Sudakov resummation, we have recovered the linear γ increase found in [57], in contrast to [1] in which P LL ∝ γ 2 has been found (see App. F for more details on this discrepancy). In the left panel of Fig. 6 we display the resulting pressure Eq. (97) at the time of bubble wall collision, using the full analytical result for ∆p = ∆p R + ∆p T in Eqs. (76) and (78). In that figure, we also show the regime of extremely small T nuc /T start , where particles are too soft in the wall frame to enter the broken phase and so the pressure decreases faster, see App G.2 for more details.
Confining PTs. We review here the findings of Ref. [12] about the pressure in confining PTs. For confining PTs, the friction pressure is not due to particles getting a mass but instead due to particles becoming strongly-coupled. At small supercooling T nuc T start , one expects the pressure to be controlled by the formation of bound-states, and to conserve the scaling of Eq. (96) where the sum is operated over a spectrum of bound-states of mass m i which depends on the model. At large supercooling T nuc T start , particles entering the bubble are separated by d ∼ 1/T nuc and therefore are far from each other compared to the confining distance f −1 with f ≡ φ . Since the confining force grows linearly with the distance d, confinement with closest neighbors through flux-tube would cost too much energy. Instead, one expects particles entering the bubble to form flux-tube attached to the wall, where the confining scale f is at its weakest value. Those strings are expected to fragments into bound states. Besides, in order to conserve color charge, during this process, a particle must be ejected from the wall. The conversion of the momentum p a of the incoming particles into string confining energy plus the ejection of a quark is expected to reduce the momentum of the wall by an amount where E string is the center-of-mass energy of the string One obtains the retarding pressure where g TC = g g + 3gq 4 with g g (g q ) the relativistic number of techni-quarks (techni-gluons). We conclude that the retarding pressure on bubble walls of confining PTs has the same γ scaling as non-confining PTs.

The terminal Lorentz factor
Thanks to the dependence of the pressure on γ, bubble walls cannot be accelerated forever but instead they reach a terminal Lorentz factor, γ LL , when the driving pressure from the vacuum energy ∆V is compensated by the friction pressure where P LO and P LL are given by Eq. (96) and Eq. (98) for non-confining PTs. Introducing P LL ≡ γP LL , we obtain Right: Bubble wall Lorentz factor γ computed with the results of this paper, see Eq. (106). On the left of the peak, bubble walls collide before they reach their terminal velocity γ LL determined by the friction pressure. This is the so-called run-away regime. The associated γ at collision is proportional to T nuc , which is the inverse of the bubble size at nucleation time, see Eq. (105). On the right of the peak, bubble walls reach their terminal velocity before collision, which is where results from this study matter.
where in the second line we assumed ∆V P LO . In the non-abelian scenario, the value of the IR cut-off µ depends on γ, cf. Eq. (57), such that the computation of γ LL in Eq. (104) might require iterations.
Before pressures equilibrate, γ grows linearly with the bubble radius, see App. H. It can happen that bubble wall collisions occur before γ saturates to its terminal value γ LL . This is the so-called run-away regime. In that case, the bubble wall Lorentz factor at the time of collision is given by (see App. H) where L w,tot = c w /T nuc , and in the following we use c w = 1 for simplicity. In the general case, the Lorentz factor at collision time is given by We show the bubble wall Lorentz factor in Fig. 7. In the yellow region, bubble walls collide before they reach their terminal velocity γ LL . This is the so-called run-away regime in which GW are dominantly produced by the scalar field gradient. In contrast, in the green region, bubble wall reach their terminal velocity before collision and the GW signal is dominated by sound-waves and turbulence. In the right-hand panel, α ≡ ∆V /ρ rad = (T start /T nuc ) 4 is the latent heat of the phase transition in the Bag model [49], ∆V is the vacuum energy of the phase transition. For P LO we use Eq. (96) and account for the contributions from t, W ± , Z and h.

Source of the GW signal
Whether the walls run away or not changes the energy budget of the expanding bubbles drastically [49,80], and thus changes the dominant contribution to the GW production. If γ run γ LL , most of the vacuum energy is used for accelerating the bubble walls and the dominant source of GW is the anisotropic part of the stress-energy tensor of the scalar field kinetic term. In contrast, if γ run γ LL , then most of the vacuum energy is converted into thermal and kinetic energy of the thermal plasma through friction, leading to a GW spectrum dominated by the contribution from sound waves and turbulence, though the fate of the highly relativistic fluid must be investigated carefully [98,99]. 8 We do not report the formula for the GW spectrum here but we instead refer to the reviews [37,39]. The classification of strong 1stOPT according to their GW sources can be visualized in Fig. 7 in the case of minimal electroweak phase transition models with φ = 174 GeV, and in Fig. 8 in the case of strong 1stOPT of arbitrary scale φ . We assumed that the pressure P is given by Eq. (95). We can see that GW are sourced by scalar field in presence of large supercooling, or when the 1stOPT occurs at very high energies φ 10 10 GeV. On the right of the blue dashed line drawn for α = 0.1, effects due to multiple reflections of vector bosons in front of the wall must be considered, see Sec. 4.4 and App. D, but we leave them for further studies. On the left of the purple dashed lines the pressure is smaller than in Eq. (98), see App. G.2, so that bubble walls still run away.

Summary and outlook
Particles passing the bubble wall of a cosmological first-order phase transition undergo splitting radiation. This is analogous to the classical radiation emitted whenever a charged particle passes from one medium into another, see e.g. Chap. 13 of [108]. Radiation from particles in the cosmological bath exert a pressure on the bubble walls which affect their velocity, and in turn the physics of quantities that depend on it (gravity waves, dark matter, the baryon asymmetry, primordial black holes, topological defects, etc). In this paper we made progress in the computation of this pressure.
In Sec. 3 and App. A, we reviewed the perturbative splitting probability, which contains an IR logarithmic divergence. We improved over previous literature by discussing possible origins for the IR cutoff in Sec. 3.4, and by quantifying the effect of finite wall thickness and other approximations in App. B.
In the regime of large gauge coupling constant and/or large supercooling, the probability can exceed unity, which calls for resummation, see Table 1. In Sec. 4, we performed the resummation at all leading-log orders, at both real and virtual levels, using the master formula in Eq. (71). We found that the averaged momentum ∆p transferred to the wall is IR-dominated, more precisely it is dominantly due to radiated gauge bosons with energies in the ballpark of their mass in the broken phase (the order parameter), m c,h . We also found that the contribution from the reflected vector bosons is at least of the same order of the contribution from the transmitted ones. We pointed out additional novel effects due to the population of reflected bosons in Sec. 4.4, and we left their detailed study to future works. In Sec. 5, we confirmed our analytical result using a Monte-Carlo simulation of the splitting processes, see Fig. 4.
Based on these results, we deduced the friction pressure on the wall in Sec. 6. Our final result for the pressure, with leading logs resummed, is where γ is the Lorentz boost of the bubble wall and T nuc is the nucleation temperature, 9 g is the gauge coupling, m c,h is the mass of the gauge bosons in the broken phase and µ is an IR cutoff at most of the order of a fraction of m c,h . This result applies as long as the energy of the incoming particles in the wall frame well exceeds the particle masses. We provided more precise and ready-to-use expressions for P LL in Eq. (98) and for µ in Eq. (57). Our results are compatible with the friction pressure for confining phase transition, which also scales as P ∝ γ f T 3 nuc (f being the confining scale), as first computed in [12] with the formalism of the gluon flux tube, see Sec. 6.1. In Sec. 6.2 we finally discussed implications for the terminal Lorentz factor, which we display in Fig. 6 for different values of the energy scale of the PT. In Sec. 6.3, we deduced the source of the GW spectrum depending on the energy scale of the transition φ and on the amount of supercooling, see Fig. 8.
Our results constitute a step towards a better understanding of cosmological strong firstorder phase transitions. Future directions that may be worth a better understanding include i) effects from the fate of reflected vector bosons, like the impact on the pressure of multiple reflections, see Sec We compute the transition amplitude of the splitting radiation X(p a ) → V (p c ) X(p b ) with X being a fermion. See [88] for the pioneering paper and [57,[109][110][111] for more recent derivations.
The transition amplitude reads The squared amplitude averaged over fermion spins is 10 Note that since the momentum along z is not conserved, the Ward identity is not satisfied, see App. F.1, and we cannot use the standard replacement pol. µ * ν → −g µν + · · · . Instead, we must sum over the physical polarizations + and − [88,112]. In the basis used for writing Eq. (14), (15) and (16), the transverse polarizations of the vector boson c read where θ ≡ k ⊥ /xE a is the emission angle and µ c ≡ m c /xE a is the mass fraction. We can check that they satisfy (p c ) · p c = 0. At leading order in x 1, θ 1, we get where C abc is the charge factor of the gauge group, e.g. for SU (N ) we have C qqg = N 2 −1 2N .

A.2 Splitting φφV T
When a and b are scalars, the transition amplitude reads At leading order in x 1, θ 1, we compute One possible example is a dark U (1) D with C φφV = 1.

A.3 Splitting V T V T V T
When a and b are bosons, the transition amplitude reads where the parameters related to the gauge structure are implicit. In the basis used for writing Eq. (14), (15) and (16), the transverse polarizations of the vector boson a, b and c read 10 The quantity which we call |V | 2 is actually 1 where ± are the helicities of a, b and c, χ ≡ x/(1 − x), θ ≡ k ⊥ /xE a is the emission angle and µ a,b,c ≡ m a,b,c /xE a are the mass fractions. We compute where C abc encapsulates the parameters related to the gauge structure, e.g. for SU (N ) we have C ggg = N . Upon averaging over initial polarizations, we obtain where ± stands for the helicity of c. Inside the wall. Instead, inside the wall we first propose using the WKB approximation, which is valid in the limit p z /p z2 1. It consists of injecting χ(z) = e i (S 0 +S 1 +S 2 2 +··· ) in the differential equation and of matching the terms which are of the same order in . Then we solve the infinite set of equations order by order in . We obtain At first order in and upon neglecting the prefactor, the M-matrix defined in Eq. (30) becomes where we have used that V is z-independent, see Eq. (31). We assume a tanh wall profile where L w is the wall thickness. 11 We then obtain the following WKB phase 11 The tanh wall profile which we consider here, accounts for the initial rising part of the wall, which interpolates between mc,s and m c,h . However, it neglects the subsequent multiple wall oscillations around m c,h . The thickness of the rising part of the wall is expected to be Lw c −1/2 vac φ −1 which for supercooled phase transition can be much smaller than the total thickness of the wall accounting for the multiple oscillations, L tot w Tnuc, cf. App. A of [12]. We leave the study of the impact of multiple wall oscillations on particle splitting for future studies.
From injecting Eq. (126) into Eq. (124), we obtain [57] M 2πVδ(∆pL w , ∆ 2 pL w ), withδ As expected, in the limit of vanishing mass difference ∆ 2 pL w → 0, the functionδ(∆pL w , ∆ 2 pL w ) approaches the Dirac δ function, see Fig. 9 δ(∆pL w , and conservation of momentum along z is restored. Using that |Γ(iy)| 2 = π/[y sh(πy)] for real  y where sh x is the sinus hyperbolic function [113], we obtain The splitting probability in Eq. (29) reduces to where Π(k ⊥ ) is defined in Eq. (39) and with shc(x) ≡ shx/x. We compute the average exchanged momentum, cf. Eq. (74) with dP a→bc in Eq. (133) and ∆p i in Eq. (70) sums over both transmitted and reflected bosons. We did not find an analytical expression for ∆p in Eq. (135) so we report its numerical value in Fig. 10. We nonetheless find an analytical approximation (E i (x) is the exponential integral special function), that is valid up to O(50%). We identify three behaviors for ∆p , We conclude that ∆p decreases logarithmically as soon as L w m −1 c,h and exponentially around L w ∼ m c,h /m 2 c,s , see Fig. 10. We remind the reader that the scaling of the pressure with the Lorentz boost γ arises from the phase space integration and not from ∆p, see e.g. Eq. (97). We also comment that Eq. (137) includes both the reflected and the transmitted contributions, so it properly includes the case where an emitted particle is softer than m c,h and is thus reflected.
In concrete scenarios we expect the wall thickness L w to be of the order of the inverse mass m −1 φ of the scalar field driving the phase transition, for which we estimate m 2 φ φ 2 /2 ∆V ≡ c vac φ 4 where ∆V is the vacuum energy difference, which implies Therefore, we expect the finite wall thickness to bring only logarithmic corrections for m c,s m c,h , see Eq. (137). However, in the regime m c,s m c,h , we expect the friction pressure to receive an exponential suppression factor, see Eq. (137), in addition to the UV suppression factor contained in Π(k ⊥ ) in Eq. (39). In Fig. 3, with dotted lines we show the impact of the wall thickness on ∆p , assuming that Eq. (138) holds.
Thin-wall limit. In light of the preceding paragraph, in the regime m c,s m c,h , effects coming from the mode functions inside the wall bring at most logarithmic corrections and we can safely approximate the wall profile by a Heaviside function. This is the subject of the next section, Sec. B.1.2. WKB breaks down. In order to further motivate the use of a Heaviside function, we would like to comment about two difficulties which we have to deal with if we rely on the WKB method, and which arise when the particle c is reflected First, the hypothesis of the WKB approximation, breaks down at the turning point, p c (z) = 0, that is present when a particle is reflected, which implies that the WKB expansion χ(z) = e i (S 0 +S 1 +S 2 2 +··· ) cannot be used. Second, the prefactor p z s /p z (z) of χ c (z) in Eq. (123) becomes infinite at the turning point. Regularization of the WKB solution near turning points has a well-known solution based on the Airy equation [114]. In order to avoid those difficulties and thanks to fact the the thin wall approximation is a good one (see Eq. (137) and discussion around it), in this paper we decide to not use the WKB approximation for computing the mode functions, see next section for more details.

B.1.2 Step potential in 1D
Klein-Gordon equation in presence of a step potential. As motivated by the preceding paragraph, we choose to model the wall by a Heaviside function with j = a, b, c. From solving the Klein-Gordon equation on each sides, we obtain χ j,s (z) = B j,s e ip j,s z + C j,s e −ip j,s z when z < 0, with Note that in contrast to [1,11,57], we have included the plane wave e −ip j,s z moving in the symmetric phase direction. We normalize the incoming wave-functions to B j,s = 1. Then we impose continuity of the wave-functions and of their derivative in z = 0 and we get We recover that in the high-energy limit p j m j , the particle j gets transmitted C j,s = 0 and B j,h = 1, while in the low-energy limit p j m j , the particle j gets reflected C j,s = −1 and B j,h = 0.
M-matrix. For relativistic walls, the particles a and b always satisfy p a m a and p b m b , such that C j,s = 0, B j,h = 1, for j = 1, 2.
The M-matrix in Eq. (35) becomes with Validity of the approximation in the main text. In the main text, for the sake of simplicity we neglect the reflected wave-function C c,s = 0 and we approximate B c,h = 1 in Eq. (142), see χ c (z) in Eq. (32) and the corresponding M -matrix in Eq. (35). In Fig. 11, we show that the latter approximation underestimates the value of ∆p computed from Eq. (146), by ∼ 20 %. The validity of the approximation (C c,s = 0, B c,h = 1) was expected since among the three terms of Eq. (146), the last one dominates over the others. In Table 2 of the next section, we compare the error due to the simplification of the mode function to other sources of error, discussed in App. B.2 and in the main text.  bottom until we obtain the most refined one. The percents % show the relative differences between two consecutive lines, except for the ones on the line entitled 'total' which shows the relative difference between the simplest and the most refined estimates. The relative difference between the analytical and numerical treatment can be appreciated in Fig. 12. The values of the present table were evaluated for µ = 10 −4 m c,h and φ = TeV as well as the same values of the parameters as in Fig. (12). Corrections due to the presence of reflected bosons in front of the wall, see Sec. 6, due to multiple wall oscillations, see footnote 11, due to the longitudinal vector boson, see footnote 3, due to Bose enhancement, see footnote 2, and due to non-perturbative effects, see Sec. 3.4, are left for further works.

B.2 Beyond the relativistic-soft-collinear limit
In order to write the perturbative splitting probability in Eq. (38), in addition to the simplified mode-function discussed in the previous section, we have assumed the relativistic limit for the phase space factor 1 2p z the soft and collinear limit x 1, k ⊥ E a for the vertex function in Eq. (31), and the relativistic, soft and collinear limit for the phase of the mode function in Eq. (34) While the relativistic limit is a very good approximation for particles a and b, it becomes incorrect for particle c when E c is close to m c,h or m c,s depending on whether c is transmitted or reflected. In this appendix, we compute the error caused by all these approximations. Our most precise analytical calculation of the exchanged momentum follows from, cf. Eq. (74) ∆p R,T = dP a→bc ∆p Θ(±p 2 c,h ), with the perturbative splitting probability given by where with B h , C s given by Eq. (144) and The factor W wall , which is defined in Eq. (134), accounts for the finite wall thickness L w . Since its impact is already studied in the Sec. B.1.1, in this subsection we set it to W wall = 1. The differences between the full analytical result in Eq. (150) and the simplified one in Eq. (76), (78) are listed in Table 2 along with the respective errors. For µ = 10 −4 m c,h and φ = TeV, the simplified formula in Eqs (76) and (78) only underestimates the total ∆p = ∆p R + ∆p T by O(5%). The corrections for other values of µ/m c,h (or T nuc /T start ) are shown in Fig. 12.

B.3 Azimuthal angle and transverse recoil
Azimuthal angle. The formula for ∆p given in Eq. (66), which we rewrite here assumes that the successive emissions occur in the same (xz) plane. Instead, Eq. (157) should be replaced by with the azimuthal angle φ i generated randomly at each emission where R 2π is a random number between 0 and 2π.
Transverse recoil. Also the formula for p z c i in Eq. (68) does not account for the recoils of the successive emission on the transverse momentum of the parent particle. Indeed, k ⊥,i is the transverse momentum relative to the actual emitter and because of the successive recoils, it must differ from the absolute transverse momentumk ⊥,i relative to the initial incoming momentum p a in Eq. (14). Instead, upon taking into account successive transverse recoils, Eq. (68) becomes with The recoils on the energy of the parent particle, x j , were already discussed in Sec. 5.2 and App. G. In Fig. 12, we compare the exchanged momentum ∆p calculated with Monte-Carlo simulations, using Eq. (70) (green line) which assumes the soft-collinear limit, single plane emission and which neglects transverse recoil, with the MC using Eq. (156), (158) and (160) (gray line). We conclude that transverse recoil and azimuthal angle can be safely neglected.

C.1 Estimation of the scattering rate
In the case of a non-abelian gauge theory, when the density of the emitted IR vector bosons becomes large, one must account for the possibility of scatterings that decrease the number of initial particles, because they can deplete the number of vector bosons and in turn lower the pressure on the wall. As an example of such scatterings, here we consider 3 → 2 processes, see Fig. 13. They become important if the rate Γ 3→2 is larger than the time it takes for these particles to cross the wall L −1 The rate reads where M is the 5-vector-boson scattering amplitude [115]. From plugging Eq. (47) into Eq. (163), we obtain . (164) Assuming that the five gluons are soft and collinear with at tree-level order [115], we obtain Upon integrating Eq. (164) over the range in Eq. (22) with m c (z) = µ, we obtain the scattering rate where C.2 Impact on the IR cut-off As argued above, the value of µ = µ 3→2 for which the population of vector boson is depleted is found after requiring that 3 → 2 processes happen within a wall length L w c We obtain We find that µ 3→2 is always smaller than the scale m sat in Eq. (52) below which phase space is saturated and perturbation theory breaks down. Therefore, the IR cut-off in Eq. (170), which relies on perturbation theory and where additionally in Eq. (166) only tree-level has been included, is not trustable and in this paper for non-abelian gauge theories we instead rely on the more conservative IR cut-off m sat in Eq. (52).

D Fate of the reflected c particles
In this appendix, we compute the typical distance l from the wall beyond which the reflected c particles, cf. Sec. 4.4, have exchanged enough momentum with the incoming a particles in order to come back in the direction of the wall.

D.1 Mean free path
Elastic cross-section. We assume that the reflected c particles scatter elastically with the a particle with the differential cross section where t is the usual Mandelstram variable, with (−t) > 0. We approximate the c particles to be massless and denote by f their momentum in the wall frame, see Fig. 14-right. Plasma frame. The threshold value for the deflection angle θ is given by the condition that the velocity of the scattered c particle in the z direction becomes equal to the wall velocity, see Fig. 14-left. The Mandelstam variable t corresponding to this deflection angle reads The relevant cross section is given by integrating (−t) above this threshold value. Since the cross section is dominated by small (−t) values, it becomes σ ∼ Multiplying by the number density of the incoming particles n ∼ (ζ(3)/π 2 )g * T 3 , we obtain the typical length scale in the plasma frame, after which the c particles come back in the direction of the wall Wall frame. In the wall frame, the threshold value of the deflection angle θ beyond which a given scattered c particle comes back in the direction of the wall, is given by the condition that c moves perpendicular to the wall, see Fig. 14-right. Energy and momentum conservation for this case becomes γT − f = p cos θ th : z momentum, p sin θ th = q : x, y momentum.
The solution for γT f is The Mandelstam t for this deflection angle becomes which is consistent with the plasma frame calculation. The relevant cross section reads Multiplying by the number density of the incoming particles n ∼ ζ(3) π 2 g * γT 3 , we obtain the mean free path in the wall frame The extra 1/γ is because the thermal bath is Lorentz contracted.
Account for multiple scatterings. We pursue the discussion in the plasma frame. The average deflection angle squared θ 2 after N scatterings θ 1 , · · · , θ N , can be determined with a random walk approach. Assuming that scatterings are independent from each other, we can write θ 2 = (θ 1 + · · · + θ N ) 2 = θ 2 1 + · · · + θ 2 N . Thus, we can simply calculate the average deflection angle squared for each, and sum up. Noting that (−t) ∼ γ 2 f 2 θ 2 and dσ/d(−t) ∼ α 2 /(−t) 2 , the average deflection squared per unit distance in the z direction reads The cutoff angle θ cut is related to the IR cut-off µ discussed in Sec. 3.4 by θ cut = µ/E a . Since the particle gets caught by the wall once the deflection accumulates to θ th ∼ √ 1 − v ∼ 1/γ, the mean free path becomes We conclude that upon taking multiple scatterings into account, the mean free path l only receives a logarithmic correction such that the big picture remains unchanged.

D.2 Motivation for further studies
A given reflected c particle gets scattered by incoming a particles before traveling the mean bubble separation when which implies

E Massless vector boson scenario
In this section, we consider the emitted vector boson to remain massless in the broken phase, up to expected thermal effects. We find that the contributions to the friction pressure are equal to the LO one, up to O(ζ a ).
Perturbative splitting probability. We suppose that the vector boson mass is phaseindependent and given by the thermal contribution m c,h = m c,s = µ while m a,h = m a,s = 0 and m b,h = m b,s = 0. The WKB suppression factor introduced along Eq. (34), in the soft-collinear limit, reads so that the splitting probability in Eq. (41) becomes Exchanged momentum averaged over resummed distribution. Since emitted vector bosons do not acquire a mass in the broken phase, they are never reflected against the wall boundary. The exchanged momentum in the soft X, x i 1 and collinear K ⊥ , k ⊥,i 1 limit, reads with where we recall that K ⊥ ≡ n i=1 k ⊥,i . We average over the resummed splitting distribution in Eq. (71). At first, since ∆p 1 is independent of the splitting kinematics, we have This coincides with the LO order piece in Eq. (2). Next we compute which implies where n = dP E is the mean number of emitted bosons. We compute which in the limit n 1, implying Finally, we use Eq. (74) to compute Final results. Hence, we conclude that in the limit where the vector boson is massless and m a = m b , the exchanged momentum in Eq. (191) is equal, up to O(ζ a ), to the LO contribution This is in contrast to [1], which, in the massless vector boson limit, have found ∆p ζ a E a , see App. F.

E.2 Case m b,h > m a,h
We suppose m b,h > m a,h in Eq. (189).
Computations. As in Eq. (191), the exchange momentum can be written ∆p = ∆p 1 + ∆p 2 + ∆p 3 with 2E a if n even, The asymptotic state particles are a → bc if the number of splitting n is odd while they are a → ac if n is even. Therefore, Eq. (71) becomes where n is the mean number of emitted bosons with µ given by the vector boson thermal mass µ = m c,s or larger. Next, we use Eq. (195) and Eq. (74) to compute which in the limit n 1, implies Finally, we use Eq. (74) to compute with µ given by the vector boson thermal mass µ = m c,s or larger.
and where n is given by Eq. (202). Again, we conclude that the contributions to the exchanged momentum coming from the emission of a massless vector boson, in the case m b,h ≥ m a,h , are only O(ζ a ) corrections to the LO result.
In the case where m a,h > m b,h , the expression in Eq. (189) has a pole when ∆p z,h = A h /2E a = 0, corresponding to the possibility in the broken phase for a to decay to bc, on-shell and without the need of the presence of any wall. After subtracting that pole, we expect the NLO correction to the exchange momentum ∆p to be of the same order as Eq. (206).
possible problem with Eq. (220). We suggest that the root of the problem may lie in the effective assumptions discussed earlier, that lead to unexpectedly satisfying the Ward identities. In contrast, our splitting probability does vanish when m c,h k ⊥ or m b,h k ⊥ (but also when m i,h → m i,s ), for vector bosons c which are respectively either massive, see Eq. (41), or massless, see Eq. (189).

F.3 Splitting at all orders
The average exchange momentum. In [1], the average exchanged momentum is computed according to where the perturbative splitting probability dP E is given by Eq. ( In that case, Eq. (71) becomes  (80), which overestimates the exchanged momentum in the limit µ → 0, the MC simulation in solid accounts for the depletion of the incoming energy-momentum due to the successive boson emissions, see Sec. 5.2. This implies a maximal number of boson emission and a maximal value for the exchanged momentum compatible with energy-momentum conservation. The impact of energy depletion is higher and higher as we decrease the incoming energy E a /m c . For simplicity, m c = m c,h in this plot.

F.4 Finite wall thickness
Finally, let us add that if the result ∆p ∼ k 2 ⊥ /2xE a ∼ ζ a E a claimed by [1] was right, the infinitely-thin-wall assumption L w ∆p 1 would violently break down. Indeed, the perturbative splitting probability used by [1] and given in Eq. (220) should receive the additional factor W wall (L w , k ⊥ , x) defined in Eq. (134), which evaluated at k 2 ⊥ /xE a ∼ ζ a E a , is exponentially suppressed where X and K ⊥ are defined in Eq. (67). The Θ functions prevent the emission of more momentum and energy than what is available, i.e the momentum p z b and energy E b in Eq. (15) must remain positive. We account for the effect of backreaction with our Monte-Carlo simulation, cf. Sec. 5, in which the Theta functions in Eq. (226) enter through the boundaries of the integral in Eq. (92) and through the condition in Eq. (93).
In Fig. 16, we can see that in contrast to the analytical result (dotted line), in the numerical study in which the effect of energy depletion due to successive emissions is included (solid . Therefore, at very strong supercooling the initial particle energy E a in the wall frame becomes smaller than the particle masses m a,h , m b,h and/or m c,h in the broken phase, such that a, b and/or c are reflected. lines), the number of emissions and the resulting ∆p saturates to a plateau. This could by itself regulate the logarithmic divergence if the thermal mass as well as other possible IR cut-offs were zero, see Sec. 3.4. However, in presence of the thermally induced cut-off µ = α 1/2 T nuc , we find that backreaction can be safely neglected.

G.2 Snowplow region: E a m i,h
The expressions for the friction pressure, P LO and P LL , induced by 1 → 1 and 1 → n processes, cf. Eq. (96) and Eq. (98), are valid in the regime m a,h , m b,h E a and µ m c,h E a , respectively. We now derive their expression in the regime where E a m a,h , m b,h and E a m c,h , in which the particles a, b and c are reflected by the wall boundary. Since the particles then accumulate in front of the wall, we call this the 'snowplow region'. 14 Assuming for the sake of simplicity that a and b are identical, the exchanged momentum between a and the wall is given by ∆p = 2p a , and the friction pressure at LO in the relativistic limit is, cf.
The exchanged momentum due to the Eq. (228) is shown by the dashed blue line in Fig. 17. Note that the region E a m c,h ∼ m a,h falls in the region where interactions between reflected particles and newly incoming particles can be neglected, cf. dotted purple and blue lines in Fig. 8, such that the possible corrections to the friction pressure discussed in Sec. 4.4 are not expected. It is interesting that a region E a m i,h , where particles are too soft to penetrate inside the bubble wall, exists in the large supercooling limit (as already pointed out in [12] for confining PTs). We leave its study for further works.
H Lorentz factor in the run-away regime H.1 Generic potential Energy conservation. The vacuum energy gained upon formation of a bubble of radius R is E bubble 4 3 π R 3 P, where P is the expanding pressure. As the bubble grows, the energy stored in the bubble wall is given by E wall 4 π R 2 γ σ, where σ is the surface energy of the wall in the wall frame σ ≡ ∞ 0 dr 1 2 (φ (r)) 2 + V (φ(r)) − V (φ(0)) .
Insuring energy conservation by equating Eq. (230) and Eq. (231), we obtain that in the runaway regime, the wall Lorentz factor grows linearly with the bubble radius, e.g. [80,116] In the supercooled and run-away regime, the expanding pressure is given by the vacuum energy difference P = c vac φ 4 while the surface tension can be estimated as where L tot. w is the total thickness of the wall, accounting for the rising part plus the multiple oscillations around the true vacuum value, until the oscillations become negligible. This should not be mistaken with L w introduced earlier, which is the thickness of the rising part of the wall only, see footnote 11. The Lorentz factor of the running-away bubble wall in Eq. (233) becomes γ = R 3 L tot. w .
(235) 15 Note that we recover the scaling PLL ∝ αγ 2 T 4 nuc found in [1], but only in the limit Ea m c,h .
Total wall thickness. The wall profile can be found by sticking the space-like profile with the time-like profile, see e.g. [117]. The former is solution of the euclidean equation of motion with φ E (0) = 0, and lim r→∞ φ E (r) = 0.
where s = √ t 2 − r 2 is now the time-like light-cone coordinate. The total wall thickness L tot. w is set by the sum of the characteristic length scales of the two profiles which are themselves set by the damping terms of the two respective Eqs. (236) and (237). Due to the matching condition φ (0) = φ E (0), the two damping terms are of the same order and we conclude that where the last equality defines the bubble radius at nucleation. The Lorentz factor of the running-away bubble wall in Eq. (235) becomes 16

H.2 Shallow potential
The friction pressure in Eq. (95) is suppressed by powers of T nuc /f 1, where f ≡ φ . Therefore, a run-away regime at bubble collision time requires supercooled phase transitions. Typically, supercooled phase transitions are generated by shallow zero-temperature potentials, namely potentials where the curvature close to the false vacuum is much smaller than the curvature f close to the true vacuum. In that case, the tunneling exit point is very close to the false minimum and the bounce action is only sensitive to the scale f exp(−c/ ). We expect both temperature and bubble radius at nucleation to be related to that scale, and therefore 17 R nuc c w T −1 nuc , 16 In [80], the wall tension is expressed in term of the bubble radius at nucleation Rnuc, obtained from minimizing the total energy E bubble + E wall σ = RnucP/2, such that the Lorentz factor of the running-away bubble wall in Eq. (233) becomes γ = 2R 3Rnuc .
However, as the bubble expands, the scalar field undergoes damped oscillations toward the true vacuum, such that Eq. (240) under-estimates the surface tension σ. With Eq. (242), we claim that Eq. (240) only under-estimates σ by an O(1) factor. 17 Note that the authors of [39] have set Rnuc ∼ f −1 instead of Rnuc ∼ T −1 nuc . For typical potentials leading to large supercooling, e.g. Coleman-Weinberg [76] or light-dilation [74], the latter choice is the correct one, see App. A of [12].
where c w is a model-dependent numerical factor. At the time of collision, the bubble size is given by 18 R coll (8π) 1/3 v w β −1 , where v w is the bubble wall velocity, β is the expansion coefficient of the exponent of the nucleation rate around the typical transition time t * Γ ∝ e β(t−t * )+··· , which gives the inverse of the bubble propagation time. Therefore, in the run-away regime (v w = 1 and P > 0), the Lorentz factor at the time of collision in Eq. (242) becomes which we have use in the main text, cf. Eq. (105). 18 The factor 8π comes from the expected number of bubbles nucleating in a volume V . To derive this, one has to take into account the fact that the nucleating bubble must be in the false vacuum. Taking the wall velocity to be unity for simplicity, and taking the nucleation rate per unit time and volume to be Γ(t) = Γ * e βt , we sum up the differential nucleation probability for time interval [tn, tn + dtn] to get (expected # of bubbles) = V × (prob. for the nucleation point to be in the false vacuum) × (prob. for nucleation between [tn, tn + dtn])