# Particle shells from relativistic bubble walls

Iason Baldes, Maximilian Dichtl, Yann Gouttenoire, Filippo Sala

► To cite this version:

Iason Baldes, Maximilian Dichtl, Yann Gouttenoire, Filippo Sala. Particle shells from relativistic bubble walls. Journal of High Energy Physics, 2024, 2024 (07), pp.231. 10.1007/JHEP07(2024)231. hal-04525393

HAL Id: hal-04525393

https://hal.science/hal-04525393v1

Submitted on 20 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RECEIVED: May 22, 2024

ACCEPTED: June 11, 2024

PUBLISHED: July 24, 2024

# Particle shells from relativistic bubble walls

Iason Baldes<sup>D</sup>,<sup>a</sup> Maximilian Dichtl,<sup>b</sup> Yann Gouttenoire<sup>D c</sup> and Filippo Sala<sup>D</sup>

E-mail: iasonbaldes@gmail.com, maximilian.dichtl@lpthe.jussieu.fr, yann.gouttenoire@gmail.com, f.sala@unibo.it

**ABSTRACT:** Relativistic bubble walls from cosmological phase transitions (PT) necessarily accumulate expanding shells of particles. We systematically characterize shell properties, and identify and calculate the processes that prevent them from free streaming: phase-space saturation effects, out-of-equilibrium 2 → 2 and 3 → 2 shell-shell and shell-bath interactions, and shell interactions with bubble walls. We find that shells do not free stream in scenarios widely studied in the literature, where standard predictions will need to be reevaluated, including those of bubble wall velocities, gravitational waves (GW) and particle production. Our results support the use of bulk-flow GW predictions in all regions where shells free stream, irrespectively of whether or not the latent heat is mostly converted in the scalar field gradient.

**KEYWORDS:** Early Universe Particle Physics, Phase Transitions in the Early Universe

ARXIV ePRINT: 2403.05615

<footnote><sup>a</sup>Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France</footnote>

<footnote><sup>b</sup>Laboratoire de Physique Théorique et Hautes Énergies, CNRS, Sorbonne Université, F-75005 Paris, France</footnote>

<footnote><sup>c</sup>School of Physics and Astronomy, Tel-Aviv University, Tel-Aviv 69978, Israel</footnote>

<footnote><sup>d</sup>Dipartimento di Fisica e Astronomia, Università di Bologna and INFN sezione di Bologna, via Irnerio 46, 40126 Bologna, Italy</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Shell production</td><td>3</td></tr><tr><td>2.1</td><td>Phase transition parameters</td><td>3</td></tr><tr><td>2.2</td><td>Context: relativistic bubble walls</td><td>5</td></tr><tr><td>2.3</td><td>Shell production mechanisms and their properties</td><td>7</td></tr><tr><td>2.4</td><td>IR cut-off</td><td>15</td></tr><tr><td>2.5</td><td>Approach of the following calculations</td><td>17</td></tr><tr><td>3</td><td>Phase space saturation</td><td>17</td></tr><tr><td>3.1</td><td>Large plasma mass</td><td>19</td></tr><tr><td>3.2</td><td>Bose enhancement</td><td>19</td></tr><tr><td>3.3</td><td>Perturbativity break-down</td><td>20</td></tr><tr><td>4</td><td>Shell momentum loss</td><td>22</td></tr><tr><td>4.1</td><td>Possible processes inducing momentum losses</td><td>22</td></tr><tr><td>4.2</td><td>Total energy considerations</td><td>23</td></tr><tr><td>4.3</td><td>Shell reversal</td><td>26</td></tr><tr><td>4.4</td><td>Shell dissipation/metamorphosis</td><td>33</td></tr><tr><td>5</td><td>Thermalization</td><td>40</td></tr><tr><td>5.1</td><td>Processes changing individual but not total numbers</td><td>40</td></tr><tr><td>5.2</td><td>Number changing processes: setup</td><td>43</td></tr><tr><td>5.3</td><td>Computation of 3 → 2 rates</td><td>45</td></tr><tr><td>5.4</td><td>Results</td><td>48</td></tr><tr><td>6</td><td>Interaction with bubble walls</td><td>49</td></tr><tr><td>6.1</td><td>Wall-wall collision</td><td>49</td></tr><tr><td>6.2</td><td>Wall-shell collision</td><td>53</td></tr><tr><td>7</td><td>Gravitational waves from free-streaming shells</td><td>55</td></tr><tr><td>8</td><td>Conclusions</td><td>57</td></tr><tr><td>A</td><td>Shell production mechanisms</td><td>59</td></tr><tr><td>A.1</td><td>Radiated vector bosons</td><td>59</td></tr><tr><td>A.2</td><td>Large IR cut-off limit</td><td>62</td></tr><tr><td>A.3</td><td>Radiated scalar bosons</td><td>63</td></tr><tr><td>B</td><td>3 → 2 final phase space integration</td><td>64</td></tr><tr><td>B.1</td><td>Partial fraction decomposition of scalar products</td><td>65</td></tr><tr><td>B.2</td><td>Parametrization of the integration region</td><td>66</td></tr></table>
B.3 Basis integrals

67

C More on the computation of $\int d\Omega |\mathcal{M}|^2_{3\to 2}$

70

C.1 Detailed results of our computation

70

C.2 Simple upper estimates of the integrated amplitudes

70

# 1 Introduction

Cosmological phase transitions (PTs) constitute some of the most prominent changes that our universe underwent in its early age. The PTs predicted by the standard model (SM) of particle physics, associated with QCD confinement and electroweak (EW) symmetry breaking, are crossovers [1, 2]. Physics beyond the SM (BSM) can make the EW PT first order, a possibility that has been extensively studied to realise electroweak baryogenesis [3]. Other first-order PTs are predicted in extensions of the SM that explain, for example, the origin of the weak scale [4–7], neutrino oscillations [8], the strong CP [9, 10] and flavour [11] problems. First-order PTs have recently attracted an enormous amount of interest because they can source a spectrum of gravitational waves [12, 13] that is observable in foreseen experiments, such as LISA [14]. Speculatively, such a PT could already have been observed by pulsar timing arrays [15–17].

First-order PTs happen via the nucleation of bubbles of the broken phase of some symmetry, into the cosmological bath which is still sitting in the unbroken phase, for pedagogical introductions see [18, 19]. Phase transitions with relativistic bubble walls, in particular, have the brightest detection prospects in GWs and offer unique phenomenological possibilities, for example, as ultra-high energy colliders in the early universe [20], in the production of dark matter [21–27], the baryon asymmetry of the universe [28–32] and primordial black holes [17, 33–42].

The evolution of the relativistic bubble walls, and of the surrounding bath, has crucial implications in all the PT properties mentioned above. The evolution of relativistic walls has been the object of several studies, which identified and computed different sources of pressure, which in turn determine the wall velocities [43–49]. Much less attention has so far been given to the evolution of the shells of energetic particles, that necessarily accumulate around relativistic bubble walls due to their interactions with the cosmological bath, for example due to splitting radiation at the bubble walls [46]. To the best of our knowledge, both the origin and the evolution of such shells have been studied so far only partially, and in a few specific cases [23, 26, 46, 48]. Progress in this direction is a necessary ingredient to reliably compute many observational consequences of PTs, including GWs [50], and could lead to spectacular phenomenological consequences, such as heavy particle production [20].

In this paper we perform the first systematic study of the evolution of particle shells around relativistic walls of first-order PTs. We then identify the possible interactions that shells undergo until they collide with those from other bubbles, and compute the conditions that allow shells to free stream without interacting nor affecting the mechanisms that source them, thus conserving both their momenta and number densities. Our main results are summarised in figure 1.
The structure of the paper is as follows. In section 2 we provide a short review of first-order PTs, with a focus on wall velocities, and we list the different kinds of shells at bubble walls and we derive their properties. We then study the processes that could prevent shells from free streaming:

* In section 3 we study phase space saturation of the shells. This is associated with large finite-density corrections, that backreact on the shells' production mechanism and in turn affect the computation of all other processes considered in this paper;

* In section 4 we study momentum changes of the shell and bath particles due to $2 \to 2$ scatterings. They could for example push the shell inside the bubble and/or significantly change its properties. We also provide a general argument why such elastic scatterings cannot prevent bath particles from entering the bubbles;

* In section 5 we study thermalisation (i.e. efficient $3 \to 2$ number-changing interactions) of shells within themselves and with the bath, which would also affect the properties of both;

* In section 6 we study the shell interactions with the collided bubble walls.

In section 7 we estimate the GW spectrum in the free-streaming regime and in section 8 we conclude.

# 2 Shell production

## 2.1 Phase transition parameters

We consider a cosmological first-order phase transition taking place in the early universe with latent heat

$$ \alpha \equiv \frac{\Delta V}{\rho_{\text{rad}}} \Big|_{T_n} \equiv \left( \frac{T_{\text{eq}}}{T_n} \right)^4, \qquad (2.1) $$

where $T_n$ is the nucleation temperature and $T_{\text{eq}}$ is the temperature when the radiation energy density, $\rho_{\text{rad}}$, would drop below the vacuum energy difference $\Delta V$ of the potential at zero-temperature. That is to say,

$$ \Delta V \equiv \frac{\pi^2}{30} g_* T_{\text{eq}}^4, \qquad (2.2) $$

where $g_*$ is the effective number of degrees of freedom (dofs). Denoting $\Delta V = c_{\text{vac}}v_{\phi}^4$, where $v_{\phi}$ is the phase transition scale (i.e. the vacuum expectation value gained by the field driving the transition), we can write

$$ T_{\text{eq}} = \left( \frac{30c_{\text{vac}}}{g_*\pi^2} \right)^{1/4} v_{\phi}. \qquad (2.3) $$

Because of its connection to the latent heat parameter $\alpha$, we show all our plots as function of $T_{\text{eq}}/T_n$. It can be useful to note its relation with another PT strength parameter, $v_{\phi}/T_n$, used in the literature

$$ \frac{v_{\phi}}{T_n} = 4.3 \left( \frac{g_*}{106.75} \right)^{1/4} \left( \frac{0.1}{c_{\text{vac}}} \right)^{1/4} \frac{T_{\text{eq}}}{T_n}. \qquad (2.4) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** *Summary plots:* in the shaded colored areas the particle shells around relativistic bubble walls cannot be considered as free-streaming until wall collision, and thus their evolution will affect standard PT predictions like wall velocities, GW and particle production. We consider shell particles produced from LO interaction (top), Bremsstrahlung radiation, either transmitted (middle) or reflected (bottom), assuming an abelian (left) or non-abelian (right) gauge interaction, see the first two rows of table 1, for small ($T_n/T_{eq} = 1$ dashed) and large ($T_n/T_{eq} = 10^{-2}$, solid) supercooling. In the blue regions finite-density corrections become large and shells can not be described by perturbative theory anymore, see section 3 (eqs. (3.14) and (3.15)). In the orange region, momentum exchange between the shell and the thermal bath significantly affects the momentum of the shell particles in either the wall frame (eqs. (4.43) and (4.44), relevant only for reflected shells) or the bath frame ((eqs. (4.73) and (4.74))), see section 4. In the red region, 3-to-2 interactions within the shells become efficient and shells are expected to thermalise, see section 5 (eqs. (5.10) and (5.12)). For the readers interested in free streaming up until shell-shell collision (as opposed to up until wall-wall collision) see section 6 and the additional shaded regions in figure 17. In the brown region, relevant only to radiated shells in the non-abelian case, finite-density corrections drive the vector boson mass $m_{c,s} = \mu$ in the symmetric phase, cf. eq. (2.38), above its vacuum value $m_{c,h}^{\infty} = g_D v_{\phi} / \sqrt{2}$ in the Higgs (i.e. broken) phase. This cuts off vector boson production, see section 3 and appendix A.2, and our treatment breaks down. Similarly, our treatment breaks down in the pink region, where the IR cutoff due to the bubble size, $\beta$, is larger than $m_{c,h}$. In the grey region, the incoming bath particle energy, $E_a$, is insufficient to create the shell particle. In the black region relativistic LO shells do not exist (see section 2.3.2).</figcaption></figure>
Scenarios with $\alpha > 1$ are accompanied with a period of supercooling during which the universe inflates for a number of e-folds,

$$N_e \equiv \log \left( \frac{T_{\text{eq}}}{T_n} \right) = \frac{\log(\alpha)}{4}. \qquad (2.5)$$

The size of bubbles at collision, $R_c$, is related to the time derivative of the nucleation rate per unit volume, $\Gamma_V$, through [51]$^1$

$$R_c \approx \frac{\pi^{1/3}}{\beta}, \quad \text{where} \quad \beta \equiv \frac{1}{\Gamma_V} \frac{d\Gamma_V}{dt}. \qquad (2.6)$$

## 2.2 Context: relativistic bubble walls

Bubble walls reach relativistic velocities, $v_w \simeq 1$, if the driving pressure $\Delta V$ is larger than the leading order pressure, $P_{LO}$. This pressure comes from particles gaining a mass across the bubble wall, and reads [43]

$$P_{LO} = \frac{g_*}{24} \Delta m^2 T_n^2, \qquad (2.7)$$

where, for convenience, we have introduced the mass difference $\Delta m$ averaged over the population of the thermal bath,

$$\Delta m^2 \equiv \sum_i c_i \frac{g_i}{g_*} \Delta m_i^2, \quad \text{with} \quad g_* \equiv \sum_i g_i, \qquad (2.8)$$

and where $\Delta m_i$ is the mass acquired by species $i$ inside the bubble, $g_i$ is its number of relativistic dofs, and $c_i = 1 (1/2)$ for bosons (fermions). In this work, we suppose bubble walls to be relativistic, $\gamma_w \gg 1$, which implies that the following condition is satisfied

$$\alpha \gtrsim \alpha_{LO} \equiv 0.1 \left( \frac{\Delta m}{T_n} \right)^2. \qquad (2.9)$$

In fact, the Bodeker and Moore criterion [43] for bubble walls to be relativistic, defined by eq. (2.9), has a caveat. Friction pressure is non-monotonic and features a peak at the Jouguet velocity ($\sim$ speed of sound) [49, 52–55]. This is due to the presence of a compression wave heating the thermal bath in front of the wall, called hydrodynamic obstruction [52].$^2$ In principle, this peaked pressure can stop the wall from accelerating around the Jouguet velocity [56]

$$\xi_J = \frac{\sqrt{\frac{2}{3}\alpha + \alpha^2} + \sqrt{1/3}}{1 + \alpha}, \qquad (2.10)$$

even if eq. (2.9) is satisfied. In models that minimally extend the Standard Model, such as by adding a single singlet scalar, the enhancement of pressure has been found effective in

<footnote>$^1$The average distance between nucleation sites is $2R_c$.</footnote>

<footnote>$^2$A peak in the pressure as a function of velocity can also arise from a large reflection probability of the longitudinal component of vector bosons, in the particular case that most of the vector boson mass is already different from zero in the unbroken phase [48]. If, in addition, the parameters of the PT are such that this specific pressure is the one stopping the wall acceleration, then those reflected vector bosons also lead to a shell propagating ahead of the wall, which we leave for future study.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Bubble wall Lorentz factor $\gamma_w$ in eq. (2.14) as a function of the scalar field vacuum expectation value $v_\phi$ or confining scale for strongly-coupled theories (left) and the ratio of the nucleation temperature $T_n$ over the temperature $T_{eq}$ when the universe would become vacuum dominated (right). The supercooling ratio can be linked to the so-called latent heat fraction $\alpha = (T_{eq}/T_n)^4$, up to changes in dofs. In both panels we assume $\alpha \gtrsim \alpha_{LO}$ in eq. (2.9) such that walls reach relativistic velocities. The gauge coupling constant $g_D$ dictates the magnitude of quantum corrections to the frictional pressure, while the phase transition completion rate $\beta/H$ determines the bubble size upon collision.</figcaption></figure>

halting wall acceleration only for deflagration-type PTs with strength parameter $\alpha \lesssim 0.01$ [54]. So in PTs with $\alpha \gtrsim 0.01$ walls would start running away, until possible effects beyond the hydrodynamic regime (discussed later) become important. When instead the drop in number of relativistic dofs from the symmetric to broken phase is large (to be precise larger than about 15%, as it can happen e.g. in confining PTs), the recent study [49] suggests that an hydrodynamic obstruction to wall acceleration could arise also for $\alpha$ up to order 1, so reducing the parameter space where walls become relativistic. However, the results from [49] could be affected when going beyond the stationary regime which they presuppose. Indeed, upon accounting for transitory hydrodynamic effects, the even more recent ref. [57] finds that bubbles start to run away in much larger regions of parameter space than predicted in the stationary regime, i.e. at $\alpha$ way smaller than 1, thus raising a doubt about the solidity of the conclusions that one can draw in the stationary regime [49]. In this study, we assume that the pressure is monotonic, which essentially means applying the Bodeker and Moore criterion of eq. (2.9), while keeping in mind these interesting ongoing developments.

Assuming the Bodeker and Moore criterion is satisfied, $\alpha \gtrsim \alpha_{LO}$, bubble walls reach ultra-relativistic velocities $v_w \simeq 1$ and two distinct possibilities present themselves:

◇ Either the Lorentz factor $\gamma_w$ of the wall grows linearly with time until bubbles collide, with [44, 46]

$$ \gamma_{\text{run}} \approx \frac{R_c}{3R_n} \approx \frac{\pi^{\frac{1}{3}} T_n}{3c_w \beta} \approx 3 \times 10^{10} \left(\frac{0.5}{c_w}\right) \left(\frac{100}{g_*}\right)^{\frac{1}{4}} \left(\frac{T_n}{10^{-4} T_{\text{eq}}}\right) \left(\frac{10}{\beta/H_*}\right) \left(\frac{\text{TeV}}{v_\phi}\right) \left(\frac{0.1 v_\phi^4}{\Delta V}\right)^{\frac{1}{4}}, \quad (2.11) $$

where $R_c \simeq \pi^{1/3}/\beta$, $R_n = c_w/T_n$ and $c_w$ is a model-dependent order-one factor, see [23] for examples of its calculation.
<figure><img src="image_5.png" /><figcaption>Figure 3. In shaded regions, the friction pressure causes bubble walls to reach a constant Lorentz factor before collision, for different values of the supercooling ratio $T_n/T_{eq}$ (left) and the PT scale $v_\phi$ (right). $\gamma_{run}$ and $\gamma_{LL}$ are the Lorentz factors at collision in the run-away and friction-dominated regime respectively, cf. eqs. (2.11) and (2.13).</figcaption></figure>

◇ Alternatively, the driving pressure on the wall, $\Delta V$, can become balanced by quantum corrections to the friction pressure [44]. Resummed at leading-log order, the corrections read [46, 47]

$$ \mathcal{P}_{LL} \approx \gamma_w \alpha_D \Delta m_V T_n^3 \log \left( \frac{v_\phi}{T_n} \right), \qquad (2.12) $$

where $\Delta m_V$ is the mass difference of a vector boson across the bubble wall, and $\alpha_D \equiv g_D^2/4\pi$ is the associated fine structure constant. The balance occurs if $\gamma_w$ reaches

$$ \gamma_{LL} \approx \frac{3.3 \times 10^{13}}{\log(v_{\phi}/T_n)} \times \left(\frac{g_{*}}{100}\right)^{3/4} \left(\frac{1/30}{\alpha_{D}}\right)^{3/2} \left(\frac{10^{-4} T_{\text{start}}}{T_n}\right)^3 \left(\frac{\sqrt{2\pi\alpha_{D}}v_{\phi}}{\Delta m_V}\right) \left(\frac{\Delta V}{0.1 v_{\phi}^4}\right)^{1/4}, \qquad (2.13) $$

after which the Lorentz factor no longer increases.

The Lorentz factor at collision is therefore the minimum of eq. (2.11) and eq. (2.13),

$$ \gamma_{\text{coll}} \approx \frac{\gamma_{\text{LL}} \gamma_{\text{run}}}{\gamma_{\text{LL}} + \gamma_{\text{run}}} \qquad (2.14) $$

The bubble wall Lorentz factor $\gamma_w$ is shown in figure 2 as a function of the PT scale $v_\phi$ and the supercooling ratio $T_n/T_{eq}$. The regions in which bubble walls reach a terminal velocity are shaded in figure 3. In those regions, most of the latent heat $\Delta V$ of the PT is converted into the particle shells. Hence, they would become the main source for GW emission.

## 2.3 Shell production mechanisms and their properties

Diverse mechanisms can produce shells of primary particles X, or c, propagating either in front or behind bubble walls (to avoid clutter, in what follows we will dub kinematic quantities related to shell particles with either only X or only c). In section 2.3.1 we discuss the origin of the shells and report from the literature the typical values of the momenta and of the number of particles that constitute them. We then derive their number density in section 2.3.3 and their thickness in section 2.3.4. These will constitute necessary ingredients to determine the rates of the various interactions that govern the evolution of shells. These of course further depend on the interactions under consideration, that we will discuss in the rest of this paper.
<table><caption>**Table 1.** Estimates of quantities $\mathcal{N}$, $p_X$ and $\bar{L}_b$ that characterize the particle shells surrounding relativistic bubble walls, generated by a variety of selected channels among those described in section 2.3. All quantities are expressed in the bath frame and are valid for free-streaming shells. These quantities are necessary ingredients to determine whether shells free stream or not, e.g. the number density of shells is determined as a function of $\mathcal{N}$, $p_X$ and $L_b = \sqrt{\bar{L}_b^2 + 1/p_X^2}$ via eq. (2.24). The bubble radius at collision is $R_c$, and we define $L_m \equiv \log(m_{c,h}/m_{c,s})$, and $L_E \equiv \log(E_a/m_{c,h})$, where $E_a \simeq 3\gamma_w T_n$ is the energy of the incoming bath particle in the wall frame, and $m_{c,s}$, $m_{c,h}$, are the masses of the radiated particle in the symmetric and Higgs phases respectively. The order parameter of the PT is $v_\phi$, either the vacuum expectation value of the scalar field driving the PT or the confining scale. The effective shell thickness $\bar{L}_b$ and the momentum of the shell particles, denoted $p_c$ or $p_X$ depending on the shell under investigation, are given in the frame of the thermal bath, i.e. the one of the CMB. For derivation of these quantities see our discussion in section 2.3 plus the references/appendices reported in the table for each channel.</caption><thead><tr><th>Channel</th><th>Multiplicity <i>N</i> per incoming particle</th><th>Momentum of shell particles (<i>p</i><sub>c</sub> or <i>p</i><sub>X</sub>)</th><th><span style="text-decoration: overline;">L</span><sub>b</sub> = (<i>L</i><sub>b</sub><sup>2</sup> - <sup>1</sup>&frasl;<sub><i>p</i><sub>X</sub><sup>2</sup></sub>)<sup>1/2</sup><br>(<i>L</i><sub>b</sub> = effective shell thickness)</th></tr></thead><tbody><tr><td>Leading-order interaction (LO):<br><i>a</i> &rarr; <i>a</i><br>Particles acquiring a mass [43, 50]</td><td>1</td><td>&Delta;<i>m</i><sup>2</sup>/<i>T</i><sub>n</sub></td><td><i>R</i><sub>c</sub> / (2(&Delta;<i>m</i>/<i>T</i><sub>n</sub>)<sup>2</sup>)</td></tr><tr><td rowspan="2">Gauge interaction &alpha;<sub>D</sub> &lt;&lt; 4&pi;:<br>Bremsstrahlung radiation<br><i>a</i> &rarr; <i>bc</i><br>[44-47] and appendix A.1</td><td>transmitted</td><td>2<sup>&alpha;<sub>D</sub></sup>&frasl;<sub>&pi;</sub><i>L</i><sub>m</sub><i>L</i><sub>E</sub></td><td rowspan="2">&gamma;<sub>w</sub><i>m</i><sub>c,h</sub></td><td rowspan="2"><i>R</i><sub>c</sub> / (2&gamma;<sub>w</sub><sup>2</sup>)</td></tr><tr><td>reflected</td><td><sup>&alpha;<sub>D</sub></sup>&frasl;<sub>&pi;</sub><i>L</i><sub>m</sub><sup>2</sup></td></tr><tr><td>Gauge interaction &alpha;<sub>D</sub> &simeq; 4&pi;:<br>Hadronization<br>[23]</td><td>string<br>fragmentation<br>ejected<br>quarks</td><td><sup>&alpha;<sub>D</sub></sup>&frasl;<sub>&pi;</sub><i>L</i><sub>E</sub></td><td>&gamma;<sub>w</sub><i>v</i><sub>&phi;</sub></td><td><i>R</i><sub>c</sub> / (2&gamma;<sub>w</sub><sup>2</sup>)</td></tr><tr><td rowspan="2">Scalar interaction &lambda;&phi;<sup>4</sup>/4!:<br>Scalar Bremsstrahlung<br><i>a</i> &rarr; <i>bc</i><br>Appendix A.3</td><td>transmitted</td><td>&lambda;<sup>2</sup><i>v</i><sub>&phi;</sub><sup>2</sup> / 192&pi;<sup>2</sup><i>m</i><sub>c,h</sub><sup>2</sup></td><td>&gamma;<sub>w</sub><i>m</i><sub>c,h</sub><sup>2</sup> / <i>E</i><sub>a</sub></td><td rowspan="2"><i>R</i><sub>c</sub> / (2&gamma;<sub>w</sub><sup>2</sup>)</td></tr><tr><td>reflected</td><td>&lambda;<sup>2</sup><i>v</i><sub>&phi;</sub><sup>2</sup> / (32&pi;<sup>2</sup><i>E</i><sub>a</sub><sup>2</sup>)</td><td>&gamma;<sub>w</sub><i>m</i><sub>c,h</sub></td></tr><tr><td>Heavier particle production &lambda;&phi;<sup>2</sup>X<sup>2</sup>/4<br>(Azatov-Vanvlasselaer mechanism <i>&phi;</i> &rarr; XX)<br><i>M</i><sub>X</sub> &gt;&gt; <i>v</i><sub>&phi;</sub> [45]</td><td></td><td>&lambda;<sup>2</sup><i>v</i><sub>&phi;</sub><sup>2</sup> / (192&pi;<sup>2</sup><i>M</i><sub>X</sub><sup>2</sup> &times; &Theta;(<i>&gamma;</i><sub>w</sub> - <i>M</i><sub>X</sub><sup>2</sup>/<i>T</i><sub>n</sub><i>v</i><sub>&phi;</sub>))</td><td><i>M</i><sub>X</sub><sup>2</sup>/<i>T</i><sub>n</sub></td><td><i>R</i><sub>c</sub> / (2(<i>M</i><sub>X</sub>/<i>T</i><sub>n</sub>)<sup>2</sup>)</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** *Shell productions:* particles from the bath interact with bubble walls and radiate other particles which accumulate in shells, preceding and following the walls. See table 1 for the different shell production mechanisms at bubble walls.</figcaption></figure>

### 2.3.1 Shells: origin, momenta, number of particles

We describe the shell-production mechanisms, or channels, that are known to us below. For a subset of them we provide, in table 1, quantitative estimates for properties of the associated shells that will be useful in the rest of the discussion.

◇ **LO shell.** Particles that acquire a mass when passing through the bubble wall also lose part of their momentum in the wall frame, because of energy conservation. This implies that, in the bath frame, they accumulate in a shell following the wall. The mass gain happens via tree-level interactions, like a scalar potential for spin-0 bosons, Yukawa interactions for fermions, and the covariant derivative for spin-1 bosons. We dub these shells ‘leading-order’ (LO) in relation to the associated pressure on the wall [43].

◇ **Transmitted and reflected shells.** Particles that obtain a mass when passing through the wall, $c$, can also be emitted as quantum radiation by other particles that couple to them, $a$, via a process $a \to bc$ (more than two final state particles are possible depending on the interaction). The breaking of Lorentz boosts by the wall background makes these emissions possible even if they would be kinematically forbidden in vacuum, so that e.g. they happen also for $a = b$. Quantum radiation results in a shell of particles $c$ (and/or $b$) following the wall and in one preceding the wall: the first is generated by the particles radiated with enough energy to penetrate the wall, the second by those that do not have enough energy, $E_c \lesssim m_{c,h}$, and so are reflected, as illustrated in figure 4. If the PT is associated to the breaking of a gauge group then the radiation of gauge bosons $c$, by particles $a$ charged under that group, is enhanced for small energies of
the gauge bosons. The multiplicity and momentum of particles in these shells depend on the particular interaction under study. For concreteness we report in table 1 those generated at a PT associated with

1. the breaking of an Abelian gauge group, with Lagrangian

$$ \mathcal{L} = (D_{\mu}\phi)^{\dagger}D^{\mu}\phi + \bar{\Psi}i\slashed{D}\Psi, \qquad (2.15) $$

where $\phi$ is a scalar, $\Psi$ a fermion, and $D_{\mu} = \partial_{\mu} - ig_{D}Q_{D}A_{\mu}$ is the covariant derivative with $g_{D}$ and $Q_{D}$ the respective gauge coupling and charge.

2. the breaking of a non-Abelian gauge group, with Lagrangian

$$ \mathcal{L} = (D_{\mu}\phi)^{\dagger}D^{\mu}\phi + \bar{\Psi}i\slashed{D}\Psi - \frac{1}{4}F_{\mu\nu}^{a}F^{a\mu\nu}, \qquad (2.16) $$

where $F_{\mu\nu}^{a} = \partial_{\mu}A_{\nu}^{a} - \partial_{\nu}A_{\mu}^{a} + g_{D}f^{abc}A_{\mu}^{b}A_{\nu}^{c}$ with structure constants $f^{abc}$ now emphasises the possibility of radiating gauge bosons off gauge bosons upon wall crossing. Here the covariant derivative is given by the usual expression $D_{\mu} = \partial_{\mu} - ig_{D}A_{\mu}^{a}t^{a}$ where $t^{a}$ are the generators of the representation.

3. the breaking of a symmetry that gives mass to a scalar $\phi$, with self-coupling

$$ \mathcal{L} = \frac{\lambda}{4!} \phi^4. \qquad (2.17) $$

For both the Abelian and non-Abelian cases, in the table we use $\alpha_{D} = g_{D}^{2}/4\pi$, and normalize to $Q_{D} = 1$, but in general we could have species of different charges in the cosmological thermal bath. The gauge-radiation shells are associated with the radiation pressure first evaluated in [44] and first resummed in [46] (except for longitudinal gauge bosons, first included in [47]), reviewed in appendix A.1. Neither the $\phi^4$ interaction, studied in appendix A.3, nor a Yukawa-like coupling between fermions and scalars (not shown here), give rise to an IR-enhanced radiation. Therefore the induced pressure is subleading with respect to that from radiated gauge bosons and, importantly for our paper, the associated radiated reflected shells are way less dense than those from the radiation of gauge bosons.

◇ **Shells from a confining PT.** Particle production associated to confining PTs has first been modeled in [23]. The key aspect is that, if the bubble walls are fast enough, bath particles charged under the confining gauge group attach to the wall via a fluxtube when they are swiped by it. The fluxtube then breaks forming hadrons inside the bubble, that constitute a shell following the wall, and ejecting charged particles in the deconfined phase to conserve charge, that constitute a shell preceding the wall. We refer the reader to [23, 24] for more information and for the quantitative estimates about these shells, and to [58] for a later study supporting the picture proposed in [23].

◇ **Shells of heavy particles.** Particles $X$ much heavier than the scale of the PT, $M_{X} \gg v_{\phi}$, can be produced from other particles when they enter the wall, if their interaction with $X$ 'feels' the PT (i.e. if it contains the field $\phi$). The heavy $X$ particles
then constitute a shell following bubble walls. According to the Heisenberg uncertainty principle, the maximal momentum exchange $\Delta p \simeq M_X^2/\gamma_w T_n$ in the wall frame is given by the inverse of the wall thickness $L_{wall} \simeq v_\phi^{-1}$, which leads to the Heaviside function shown in table 1. Their emission and the associated pressure has first been computed in [45], for the computation of some of their other properties and their implications for dark matter or baryogenesis see e.g. [25, 26, 30, 59]. For concreteness, in table 1 we report estimates for shell quantities associated to the interaction

$$ \mathcal{L} = \frac{\lambda}{4} \phi^2 X^2. \qquad (2.18) $$

◇ **Shells from vector bosons acquiring a small part of their mass.** These shells arise in the particular case where the mass of a gauge boson is already non-zero in the unbroken phase, and it gets an extra subleading contribution upon the PT. The pressure on walls in this specific scenario has first been evaluated in [48], together with some of the quantities that characterize the associated shell.

◇ **Shells from decay of the wall.** The background field $\phi$ constituting the wall undergoes several oscillations inside the bubble, before relaxing to $v_\phi$ deep inside it. These oscillations experience friction both by Hubble and by emissions of the particles coupled to the background field. These particles are then produced by $\phi$ oscillations and constitute a shell following the wall.

◇ **Unruh and Casimir shells.** The motion of bubble walls acts as a time-dependent boundary condition for quantum fields. Hence, energy is not conserved and particles should be produced, analogously to what happens in quantum electrodynamics [60–62]. The production of those particles can be viewed as equivalent to Unruh radiation [63] or dynamical Casimir effects [64] according to whether the bubble wall is accelerating or expanding at constant velocity.

To our knowledge, the existence of the last two kinds of shells mentioned in this list have never been pointed out in the context of cosmological PTs and we leave them for future work.

### 2.3.2 Shell's boost factor

We will be interested in shells with significant bulk motion distinguishable from the thermal bath. That is, shells in which the particles move through the bath in a direction outward from the bubble nucleation site, and with significant Lorentz factors associated with the averaged velocity vector.

Consider leading-order (LO) shells, coming from the transition of bath particles from the symmetric to the Higgs phase, $a \to a$, which gives a corresponding change in mass $m_{a,s} \to m_{a,h}$. Take the wall moving in the z-direction. Then the change in momentum in the z-direction is given $p_x \simeq \Delta m_a^2/T_n = (m_{a,h}^2 - m_{a,s}^2)/T_n$. The averaged momentum of the shell particles is therefore relativistic in the bath frame provided

$$ p_x \gtrsim m_{a,h}. \qquad (2.19) $$
Assuming this is the case, this will also be the dominant momentum component of our shell particles. If $p_X \lesssim m_{a,h}$ the particle's momentum will still be dominated by the initial value from the thermal bath, $|\vec{p}| \sim T_n$, rather than the “kick” received by the passing wall. Taking the simple example with $m_{a,h} \gg m_{a,s}$, and a gauge boson shell with $m_{a,h} = g_D v_\phi$, eq. (2.19) translates into

$$g_D \gtrsim \frac{T_n}{v_\phi} \quad (\text{LO shell}). \qquad (2.20)$$

In all figures for the LO cases we will shade the regions where eq. (2.20) is not satisfied.

Going on to other types of shells, we read their momentum, $p_X$, off table 1 and see that for perturbative and non-perturbative gauge interactions we automatically have relativistic shells, provided the wall is relativistic $\gamma_W \gg 1$. This also holds for the scalar bremsstrahlung in case of reflected scalars. For transmitted scalars we again require the final state mass to be large compared to the temperature $m_{c,h}/T_n \gtrsim 1$ for the shell to be relativistic (similarly to the LO shell). Taking the mass to be $m_{c,h} \simeq \sqrt{\lambda} v_\phi$, one finds the condition

$$\lambda_{\phi} \gtrsim \left( \frac{T_n}{v_{\phi}} \right)^2 \quad (\text{Transmitted scalar bremsstrahlung shell}). \qquad (2.21)$$

Finally for the Azatov-Vanvlasselaer mechanism we also need $M_X/T_n \gtrsim 1$. But this is automatically satisfied even in the weakly supercooled limit in this scenario, $v_\phi \gtrsim T_n$, as the scenario operates in the limit $M_X \gg v_\phi$.

### 2.3.3 Shell particles' number density

The number density of the shell particles X in the wall frame, produced at a given time $t_X$ when the bubble radius is $r_X$, is given by

$$n_{X,w} = \mathcal{N}_X^i n_{i,w}(t_X) = \gamma_w(t_X) n_{i,b}, \qquad (2.22)$$

where $\mathcal{N}_X^i$ is the number of shell particles X produced per each particle of population i upon being swept by the bubble wall (see table 1), and where $n_{i,b} = g_i \zeta(3) T_n^3 / \pi^2$ is their number density in the bath frame. The number density of shell particles in the bath frame, $n_{X,b}$, then depends on the time elapsed since their production at $t_X$, or equivalently on how much the bubble has expanded after their production. Assuming that shells free stream until the bubble radius reaches the value $R > r_X$, their number density reads

$$n_{X,b} = 2\gamma_w(R)\left(\frac{r_X}{R}\right)^2 n_{X,w} = 2\gamma_w(R)\gamma_w(r_X)\left(\frac{r_X}{R}\right)^2 \mathcal{N}_X^i n_{i,b}, \qquad (2.23)$$

where the factor $2\gamma_w(R)$ arises from boosting the shell's density current from the wall to the bath frame, and the factor $(r_X/R)^2$ accounts for the spatial dilution of the number density due to the expansion of the bubble from $r_X$ to $R$.<fnref n="3" />

<footnote n="3">³Our derivation of the shells number density and effective thickness generalises the one for ejected quanta at confining PTs of [23] in two directions: i) it is valid for any shell, the shell dependence being encoded in $\mathcal{N}_X^i$; ii) it is valid irrespective of whether walls run away or reach a terminal velocity, while [23] assumed runaway walls.</footnote>
As we will prove in section 2.3.4, an effective shell's thickness $L_b$ exists such that, for any shell, $n_{X,b}$ is approximately constant at distances from the wall shorter than $L_b$, and goes to zero at larger distances. For simplicity, in this paper we will then use a constant value for the shells number density, obtained by dividing the total number of particles in that shell, $n_i \mathcal{N}_X^i 4\pi R_c^3/3$, by the effective volume of the shell, $4\pi R_c^2 L_b$, where $R_c$ is the bubble radius at collision:

$$n_{X,b} = \sum n_i \mathcal{N}_X^i \frac{R_c}{3L_b}, \qquad (2.24)$$

where all quantities are evaluated in the bath frame. We find it useful to also report the total number of particles $X$ produced per bubble

$$N_X = \sum n_i \mathcal{N}_X^i \frac{4\pi R_c^3}{3}. \qquad (2.25)$$

### 2.3.4 Shell's thickness

We start by noticing that the distance of shells' particles from the bubble wall cannot be localized to shorter values than their de Broglie wavelength $\sim 1/p_X$. This will constitute a lower bound to the effective thickness of any shell, that we write as

$$L_b^X = \sqrt{(\bar{L}_b^X)^2 + \frac{1}{p_X^2}}, \qquad (2.26)$$

where $\bar{L}_b^X$ is the contribution to the shell's thickness that depends on the shell's production dynamics. We will find that, throughout the parameter space considered in our study, the de Broglie contribution will dominate over $\bar{L}_b^X$ at very small $g_D$ and very large $v_\phi$. We now turn to the computation of $\bar{L}_b^X$ for the various shells.

Let us define $x$ as the radial distance of a given layer of the shell from the bubble wall, which is limited between $x=0$ (i.e. the shell's layer produced last before collision), and the position of the layer emitted first, which corresponds to $x = L_b^{\max}$. $x$ can be positive or negative depending on whether the shell of interest precedes or follows the walls. All the $x$ dependence in the number density eq. (2.23) is encoded in $r_X(x)$. While $L_b^{\max}$ defines the maximal extension of a shell, we anticipate that there always exists a value of $x$, which we define as $\bar{L}_b$, below which $n_{X,b}$ goes to a constant value and above which it gets very suppressed. We now derive $r_X(x)$ and $\bar{L}_b$.

The distance $x$ between a shell particle “X” and the wall is found by integrating the difference between the world lines of “X” and the bubble wall “w” from the moment of particle production $t_X$ until the time of bubble collision $t_c$,

$$x = \int_{t_X}^{t_c} dt [v_X(t_X) - v_w(t)] \approx \int_{t_X}^{t_c} dt \left( \frac{1}{2\gamma_w^2(t)} - \frac{1}{2\gamma_X^2(t_X)} \right) = -\frac{t_c - t_X}{2\gamma_X^2(t_X)} + \int_{t_X}^{t_c} dt \frac{1}{2\gamma_w^2(t)}. \qquad (2.27)$$

where we have used $v \simeq 1 - 1/(2\gamma^2)$ in the second equality, i.e. we have assumed $\gamma_{X,w} \gg 1$.

Let us start by considering shells that follow the wall, like LO shells, shells of heavy particles, and radiated transmitted shells. For them one has $\gamma_X \ll \gamma_w$, so that the second term in eq. (2.27) for $x$ is negligible and

$$x \approx \frac{r_X - R_c}{2\gamma_X^2(r_X)} \quad (\text{shells behind wall}) \qquad (2.28)$$
irrespective of whether walls runaway until collision or not. By inverting eq. (2.28) one obtains $r_X(x)$ and, via eq. (2.23), the dependence of $n_{X,b}$ in $x$. If $\gamma_X$ is constant in $r_X$, as it is for LO shells ($\gamma_X \simeq \Delta m^2/(T_n m)$, $\simeq \Delta m/T_n$ if $m \simeq \Delta m$) heavy particle shells ($\gamma_X \simeq M_X/T_n$) and radiated transmitted ones for terminal velocity walls ($\gamma_X \lesssim \gamma_w^{LL}$), one finds

$$n_{X,b}(x) = 2\gamma_w(R_c)\gamma_w(r_X(x)) \left(\frac{r_X(x)}{R_c}\right)^2 \mathcal{N}_X^i n_{i,b} \quad (2.29)$$

$$= 2\mathcal{N}_X^i n_{i,b} \left(1 + \frac{x}{\bar{L}_b^X}\right)^2 \times \begin{cases} (\gamma_w^{LL})^2, & \text{terminal velocity,} \\ (\gamma_w^{\text{run}}(R_c))^2 \left(1 + \frac{x}{\bar{L}_b^X}\right), & \text{runaway,} \end{cases} \quad (2.30)$$

where we used $\gamma_w^{\text{run}}(r_X) = \gamma_w^{\text{run}}(R_c)(r_X/R_c)$ given by eq. (2.11) and we have defined

$$\bar{L}_b^X = \frac{R_c}{2\gamma_X^2}, \quad X = \text{shells behind wall.} \quad (2.31)$$

$\bar{L}_b^X$ contributes to what we defined as the effective thickness of shell $X$ because for values of $x > -\bar{L}_b^X$ the shells number density is approximately constant, while it goes to zero as $x$ approaches $-\bar{L}_b^X$ (we remind that $x < 0$ for shells following the wall).

For shells made of radiated transmitted particles, this time in the case that bubbles collide while walls run away, we assume for simplicity that all radiated transmitted particles are produced in the regime of runaway walls, so that $\gamma_X(r_X) \propto r_X$.<fnref n="4" /> Eq. (2.28) then becomes $x = (R_c^2/r_X^2)(r_X/R_c - 1)R_c/(2\gamma_w^{\text{run}}(R_c)^2)$, its solution $r_X(x)$ satisfies $r_X(x) \to R_c$ for $x \ll R_c/(2\gamma_w^{\text{run}}(R_c)^2)$ and $r_X(x) \to R_c\sqrt{R_c/(2\gamma_w^{\text{run}}(R_c)^2 x)}$ for $x \ll R_c/(2\gamma_w^{\text{run}}(R_c)^2)$, so that

$$n_{\text{transm.,b}}^{\text{run}}(x) = 2\gamma_w^2(R_c)\mathcal{N}_{\text{transm.}}^i n_{ib} \times \begin{cases} 1, & x \ll \bar{L}_b^X \\ \left(\frac{\bar{L}_b^X}{x}\right)^{\frac{3}{2}}, & x \gg \bar{L}_b^X \end{cases} \quad (2.32)$$

where $\bar{L}_b^X$ has again the expression of eq. (2.31), this time of course using $\gamma_X = \gamma_w^{\text{run}}(R_c)$.

Let us now consider shells that precede the wall, like those formed by radiated reflected particles, or by ejected techniquanta in a confining PT. In this case $\gamma_X = p_X/m_X \sim \gamma_w v_\phi/\mu \gg \gamma_w$, where we have anticipated from section (2.4) that in the unbroken phase $m_X \sim \mu$ where $\mu$ is the IR cutoff. Eq. (2.27) for $x$ is then dominated by the second term and it becomes

$$x \approx \frac{(R_c - r_X)}{2\gamma_w^2 \text{coll}} \times \begin{cases} 1 & \text{terminal vel.} \\ \frac{\gamma_w \text{coll}}{\gamma_w^{\text{run}}(r_X)} & \text{runaway} \end{cases} \quad \text{(shells ahead of wall),} \quad (2.33)$$

where $\gamma_w$ is given in eq. (2.14). In the runaway case, we have again assumed for simplicity that all shell's particles are produced in the regime of runaway walls, so that $\gamma_w(r_X) \propto r_X$ (see again footnote 4). Using eq. (2.23), we then find

$$n_{X,b}(x) = 2\gamma_w^2(R_c)\mathcal{N}_X^i n_{ib} \times \begin{cases} \left(1 - \frac{x}{\bar{L}_b^X}\right)^2, & \text{terminal vel.} \\ \left(1 + \frac{x}{\bar{L}_b^X}\right)^{-3}, & \text{runaway,} \end{cases} \quad (2.34)$$

<footnote n="4">⁴Most of the shell is produced in the last stages of bubble expansion, which are the stages where walls swipe through most of the volume available to a bubble before collision. Therefore this approximation is accurate up to an $O(1)$ factor in the regions of parameter space where $\gamma_w^{LL}(R_c) \simeq \gamma_w^{\text{run}}(R_c)$, and to a much better precision as soon as $\gamma_w^{LL}(R_c) > \gamma_w^{\text{run}}(R_c)$.</footnote>
where we have defined, analogously to above,

$$ \bar{L}_b^X = \frac{R_c}{2\gamma_{w,\text{coll}}^2}, \quad X = \text{shells ahead of wall.} \qquad (2.35) $$

Finally, we find it useful to report the expressions for the contributions to the effective thickness of shells, $\bar{L}_b^X$ of eqs. (2.31) and (2.35), in the wall frame. They read

$$ \bar{L}_w^X = R_c \times \begin{cases} \frac{1}{2\gamma_{w,\text{coll}}}, & X = \text{shells ahead of wall,} \\ \frac{\gamma_{w,\text{coll}}}{\gamma_X^2 + \gamma_{w,\text{coll}}^2}, & X = \text{shells behind wall,} \end{cases} \qquad (2.36) $$

where the non-trivial factor, in the case of shells behind the wall, comes from the transformation law of velocities $v_X|_{\text{wall}} = (v_X - v_w)/(1 - v_w v_X) \approx \gamma_w^2/(\gamma_X^2 + \gamma_w^2)$. If $\gamma_X \lesssim \gamma_w$, then we have $\bar{L}_w \approx R_c/(2\gamma_{w,\text{coll}})$. The expressions for the shells thickness is reported in table 1, for convenience for specific shells.

## 2.4 IR cut-off

We anticipate that most squared amplitudes of the processes that we will study have IR singularities in vacuum QFT, once integrated over the final phase space. This is unphysical however, since interactions with the plasma (be it thermal, i.e. the bath, or out of equilibrium, i.e. the shell) screen long-range forces. This is equivalent to providing an IR cutoff, which takes the form of a Debye mass appearing in the propagators, which reads [65–67]

$$ m_D^2 \approx 2g_{\text{eff}} g_D^2 \int \frac{d^3 p}{(2\pi)^3} \frac{f(p, T)}{p^0} \approx 2g_D^2 \frac{n}{\langle E \rangle}, \qquad (2.37) $$

where $g_D$ is the coupling strength of the interaction of interest, $f$ is the phase-space distribution of the plasma particles that are charged under that interaction, and $g_{\text{eff}}$ is their number of degrees of freedom weighted by their charge. For example, in the case of an equilibrium thermal distribution $f(p, T)$ of gluons and fermions charged under a gauged SU(N), this reproduces the well-known result $m_D^2 \approx g_D^2 T^2(N/3 + N_f/6)$ [68], where $N_f$ is the number of Dirac fermions in the fundamental of SU(N). In case particles in the shells are charged under the gauge interaction of interest, e.g. as in the cases of LO shells and of shells formed by reflected or transmitted non-abelian gauge bosons, then eq. (2.37) implies that $m_D^2$ receives a contribution both from bath and from shell particles, and the latter can dominate over the former [23]. For PTs occuring in a gauge sector, once bubbles are nucleated and shells are created, the IR cutoff reads

$$ \mu^2 \approx \begin{cases} \frac{N_f}{6} g_D^2 T_n^2 & \text{U(1) gauge bosons,} \\ \left(\frac{N_f}{6} + \frac{N}{3}\right) g_D^2 T_n^2 + 2g_D^2 \frac{n_s}{p_X} & \text{SU(N) gauge bosons,} \\ \frac{N_f C_f}{8} g_D^2 T_n^2 + d_f g_D^2 \frac{n_s}{p_X} & \text{fermions or scalars charged under U(1) or SU(N),} \end{cases} \qquad (2.38) $$

where the shell densities $n_s$ depend on the shell of interest and can be obtained from section 2.3.3 plus table 1. The coefficient $C_f$ is the fermion color factor [69]. It is equal to
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Phase space saturation: a large occupation number of particles in the shell can lead to Bose-enhancement of the splitting radiation process (top) or to a break-down of perturbativity (bottom). See section 3 for the details.</figcaption></figure>

$C_f = 1$ for U(1) fermion and $C_f = 4/3$ for SU(N) fermion in the fundamental representation. The factor $d_f$ is an order one factor which we take equal to $d_f = 2$ for definiteness. In the Abelian case, the IR cutoff as given in eq. (2.38) is already a function of the fundamental parameters of the PT. This is not the case for non-Abelian gauge theories, for which we give below more explicit expressions for $\mu^2$. For LO shells, using $n_s$ and $p_X$ from table 1, we find that the non-Abelian contribution scales as the Abelian one,

$$ \text{Non-abelian, LO shells: } \mu^2 \sim g_D^2 T_n^2. \qquad (2.39) $$

For shells coming from radiated gauge bosons, either transmitted or reflected, particles in the shells are squeezed and the number density $n_X$ defined in eq. (2.24) reads

$$ n_X \approx \mathcal{N} \frac{\zeta(3)}{\pi^2} g_* T_n^3 \frac{R_c}{3L_b} \approx g_D \gamma_w^2 T_n^3 \left( \frac{g_*}{12\pi^4} \right) \left( \frac{\gamma_{\text{run}}}{\gamma_w} \right) \left( \frac{\mathcal{N}}{\alpha_D / \pi} \right) \left( \frac{\bar{L}_b}{L_b} \right), \qquad (2.40) $$

where $\bar{L}_b$ and $L_b$ are related by eq. (2.26). For non-abelian gauge sector, it leads to the IR cut-off, cf. eq. (2.38)

$$ \text{Non-abelian, radiated shells: } \mu^2 \approx g_D^4 R_c \frac{T_n^4}{m_{c,h}} \left( \frac{g_*}{6\pi^4} \right) \left( \frac{\mathcal{N}}{\alpha_D / \pi} \right) \left( \frac{\gamma_w m_{c,h}}{p_X} \right) \left( \frac{\bar{L}_b}{L_b} \right). \qquad (2.41) $$

Another IR cutoff is provided by the size of the bubbles at collision,

$$ \text{All shells: } \mu^2 \approx \beta^2. \qquad (2.42) $$
If the interaction of a shell particle is mediated by frequencies larger than the inverse of the collision time, $\beta$, then that interaction cannot be effective before bubbles collide, and so it cannot affect the propagation of shells. In other words, wavelengths larger than the spatial dimension of bubbles at collision should affect the physics of bubbles nor vice-versa. The IR cutoff due to bubble size, eq. (2.42), is subleading with respect to the IR cutoffs due to plasma effects, eq. (2.38), in most of the regions of our parameter space. For example it is smaller than $g_D T_n$ if

$$g_D \frac{M_{Pl} T_n}{v_\phi v_\phi \beta/H} \gtrsim 1. \qquad (2.43)$$

We anticipate that this condition will turn out to be satisfied at the border of all free-streaming regions, so that using the cutoff from the plasma-mass to determine that border, as we will do, is self-consistent. Deep in the regions of free streaming, far from that border, there will be parts of parameter space where the inequality of eq. (2.43) is violated: this will only mean that amplitudes have a larger cutoff than the one we employ, so they are smaller, so a fortiori free-streaming in those regions is confirmed. The IR cutoff in eq. (2.42) will however play a role in identifying a region, deep in the free-streaming area, where our predictions cannot be trusted due to $\mu \simeq m_{c,h}$, see section 3.1 and the regions shaded in pink in figures 1 and 6.

## 2.5 Approach of the following calculations

In the calculation of all effects that could prevent shells from free streaming (i.e. in sections 3, 4, 5 and 6), we use the determinations of the PT and shell properties within the free-streaming assumption. This allows one to self-consistently determine the regions of parameter space of a PT where shells are guaranteed to free-stream, which is the purpose of this paper. Also, to our knowledge, no one ever determined (outside of the hydrodynamic limit) the PT and shell properties in regions where shells are relativistic and do not free stream. Our study is then not only self-consistent, but also a needed first step to determine the behaviour of walls and shells in regions of no free-streaming, which we leave for future work. For example, in those regions, the Lorentz factor $\gamma_{coll}$ in eq. (2.14), the shell number density and typical momentum, $n_X$ and $p_X$, and the shell thickness $L_w$ (see table 1) will be affected. These will in turn induce a backreaction on all calculations that depend on them.

## 3 Phase space saturation

In this section, we identify the regions where finite density corrections become significant, either affecting the plasma mass discussed in section 2.4 or the particle occupation number in the shells. These corrections can lead to Bose enhancement or perturbativity breakdown, as illustrated in figure 5. The regions where these effects occur are individually shaded in figure 6, and their combined area is shaded in figure 1.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** *Phase space saturation:* in the blue regions, the phase space is saturated and the shells - their production and propagation - can not be described by perturbative theory, see conditions in eq. (3.14) and (3.15). In the orange region, the splitting radiation probability producing the shell should be corrected by Bose-enhancement factors, see condition in eq. (3.8). In the brown region, finite-density corrections drive the vector boson mass $m_{c,s} = \mu$, cf. eq. (2.38), above its broken vacuum value $m_{c,h}^{\infty} = g_D v_{\phi}/\sqrt{2}$, cutting off vector boson production, see appendix A.2. Our treatment based on $m_{c,s} \ll m_{c,h}$, can not be trusted in the brown region and is left to future works. Similarly, our treatment breaks down in the pink region, where the IR cutoff to shell processes due to the bubble size, $\beta$, is larger than $m_{c,h}$. We consider shell particles produced from LO interaction (top), Bremsstrahlung radiation, either transmitted (middle) or reflected (bottom), assuming abelian (left) or non-abelian (right) gauge interaction, see the first two rows of table 1, assuming the amounts of supercooling $T_n/T_{eq} = 1$ (dashed) and $T_n/T_{eq} = 10^{-2}$ (solid). We used $g_b = g_{SM} + g_{emit}$ with $g_{SM} = 106.75$ and $g_{emit} = 10 = N_F$. Pink, grey and black shaded areas as in figure 1.</figcaption></figure>
## 3.1 Large plasma mass

Due to the high density of shells, the plasma mass $m_{c,s} \simeq \mu$ in the symmetric phase can become larger than the vector boson mass $m_{c,h}^\infty \simeq g_D v_\phi / \sqrt{2}$ in the broken phase far away from the wall

$$ \mu^2 \gtrsim (m_{c,h}^\infty)^2 \quad \Rightarrow \quad g_D T_n R_c \left(\frac{T_n}{v_\phi}\right)^3 \left(\frac{g_*}{200}\right) \left(\frac{N}{\alpha_D/\pi}\right) \left(\frac{\gamma_w m_{c,h}}{p_X}\right) \frac{\bar{L}_b}{L_b} \gtrsim 1. \quad (3.1) $$

Note that the vector boson mass in the broken phase close to the wall, within the shell, is also supplemented by the plasma correction

$$ m_{c,h}^2 \simeq m_{c,s}^2 + (g_D v_\phi / \sqrt{2})^2, \quad \text{with} \quad m_{c,s} \simeq \mu. \quad (3.2) $$

As a result, for non-abelian gauge group in which the vector boson mass receives finite-density correction and in the ultra-relativistic limit in which this correction is large, the mass difference $\Delta m_c = m_{c,h} - m_{c,s}$ becomes small $\Delta m_c \ll m_{c,h}$. In appendix A.2, we show that in this regime the number of particles radiated becomes suppressed by $N_R \propto (\Delta m_c / m_{c,h})^3$ and $N_T \propto (\Delta m_c / m_{c,h})^2$, according to whether they are reflected on or transmitted though the wall. This suggests that the number of particles in the shell saturates. At the same time, the momentum exchange between an incoming particle and the wall is suppressed by $\langle \Delta p_R \rangle \propto (\Delta m_c / m_{c,h})^4$ and $\langle \Delta p_T \rangle \propto (\Delta m_c / m_{c,h})^{5/2}$, see appendix A.2 again, suggesting that the friction pressure drops, leading to an enhancement of the bubble wall Lorentz factor $\gamma_w$. We shade the regions where this occurs in brown in figures 1, 6, 10, 11 and 13. Our results can not be trusted inside those regions since the number of particles in the shells and their kinematics should be modified accordingly. We leave its exploration for future works. In an entirely analogous manner, we shade in pink the regions where the IR cutoff due to finite bubble-size, eq. (2.42) is larger than $m_{c,h}$, where all the limitations mentioned above also apply. Note that this region is independent of $T_n / T_{eq}$.

## 3.2 Bose enhancement

In the wall frame, consider particle a crossing the wall, splitting to produce particle b of roughly similar momentum as a, and soft boson x which will form the shell. The interaction Hamiltonian for the splitting process can be expressed as

$$ H_{\text{int}} = \mathcal{M}_0 a_c^\dagger a_b^\dagger a_a + \text{h.c.}, \quad (3.3) $$

where $a_x$ are the creation operators in Fock space. Then the transition amplitudes for emission and absorption read, respectively

$$ \mathcal{M}_{a \to bc} = \langle f_a - 1, f_b + 1, f_c + 1 | H_{\text{int}} | f_a, f_b, f_c \rangle = \mathcal{M}_0 \sqrt{f_a} \sqrt{1 \pm f_b} \sqrt{1 + f_c}, \quad (3.4) $$

$$ \mathcal{M}_{bc \to a} = \langle f_a + 1, f_b - 1, f_c - 1 | H_{\text{int}} | f_a, f_b, f_c \rangle = \mathcal{M}_0 \sqrt{1 \pm f_a} \sqrt{f_b} \sqrt{f_c}, \quad (3.5) $$

where $+/-$ refers to boson/fermion statistic. We deduce the interaction rate accounting for both emission and absorption

$$ |\mathcal{M}_{a \to bc}|^2 - |\mathcal{M}_{bc \to a}|^2 = |\mathcal{M}_0|^2 [f_a(1 \pm f_b) + f_c(f_a - f_b)], \quad (3.6) $$
where we have used that particle c is bosonic so a and b necessarily are of the same statistical type.<fnref n="5" /> We can see that Bose-enhancement can be neglected as long as $(f_a - f_b)f_c \ll f_a$, in which case we get

$$ |\mathcal{M}_{a \to bc}|^2 - |\mathcal{M}_{bc \to a}|^2 \approx |\mathcal{M}_0|^2 f_a. \qquad (3.7) $$

We now investigate the condition for neglecting Bose-enhancement

$$ f_c(p_c)(f_a(p_a) - f_b(p_b)) = f_c(p_c)(p_a - p_b) \frac{\partial f_a}{\partial p_a} \lesssim f_a(p_a). \qquad (3.8) $$

We now add a subscript 'p' to indicate quantities evaluated in the bath frame. In the absence of 'p', quantities are evaluated per default in the wall frame. Assuming a Maxwell-Boltzmann momentum distribution in the bath frame, $f_a(p_{a,p}) \propto \exp(-p_{a,p}/T_n)$, we can write

$$ \frac{\partial f_a}{\partial p_a} = \frac{\partial p_{a,p}}{\partial p_a} \frac{\partial f_a}{\partial p_{a,p}} = - \frac{f_a(p_a)}{\gamma_w T_n}, \qquad (3.9) $$

where we used that $p_{a,p} = \gamma_w p_a - \beta \gamma_w \sqrt{m_a^2 + p_a^2} \approx p_a/\gamma_w$. Plugging eq. (3.9) into eq. (3.8) leads to

$$ f_c(p_c)(f_a(p_a) - f_b(p_b)) = f_c(p_c)f_a(p_a) \frac{\Delta m_a^2}{2p_a\gamma_w T_n}, \qquad (3.10) $$

Replacing $p_a \approx 3\gamma_w T_n$ and $f_a(p_a) \approx 1$, and using eq. (2.11) we calculate

$$ f_c(p_c)(f_a(p_a) - f_b(p_b)) \approx \frac{1}{6} f_c(p_c) \left( \frac{\Delta m_a}{R_c T_n^2} \right)^2 \left( \frac{\gamma_{\text{run}}}{\gamma_w} \right)^2 \qquad (3.11) $$

Note that eq. (3.11) is general, in the sense that it applies to particles c belonging to any shell considered in this paper. We show the parameter space where Bose-enhancement matters in orange in figure 6. In most of the parameter space, we find the condition for neglecting Bose-enhancement $f_c(p_c)(f_a(p_a) - f_b(p_b)) \ll 1$ less stringent than the condition for non-perturbativity $g_D^2 f_c(p_c) \ll 1$ or $g_D^4 f_c(p_c) \ll (4\pi)^2$, shown in blue in the same figure and discussed in the next section.

### 3.3 Perturbativity break-down

Let us consider first shells produced by splitting radiation $a \to bc$. The occupation number $f_c(p_c)$ of radiated vector bosons is related to the vector boson wave function $A_\mu$ by

$$ (\partial A)^2 \sim p^2 A^2 \sim \int_{p_c} d^3 p p f_c(p) \sim p_c^4 f_c(p_c), \qquad (3.12) $$

where, in the second wiggle, we have used that the boson kinetic term 'corresponds' to the energy of the free gas of bosons. For a non-abelian theory, the hierarchy between the Lagrangian terms

$$ \mathcal{L} \supset \partial A \partial A + g A A \partial A + g_D^2 A A A, \qquad (3.13) $$

<footnote n="5">⁵The phase space factors, $f_a$, $f_b$, $f_c$ which appear in the forward, eq. (3.4), and reverse processes, eq. (3.5), are the same ones. There is no switch in the momentum direction.</footnote>
which is essential for perturbation theory to apply, breaks down as soon as

$$ \text{Non-abelian: } g_D^2 f_c(p_c) > 1. \qquad (3.14) $$

Instead, for an abelian theory, gauge bosons self-interaction terms are loop-suppressed so that perturbativity breaks down for

$$ \text{Abelian: } \frac{g_D^4}{(4\pi)^2} f_c(p_c) > 1. \qquad (3.15) $$

Similarly, for emission of self-interacting scalars $\phi$ with $\mathcal{L} \supset \lambda\phi^4$, perturbativity breaks down as soon as

$$ \text{Self-interacting scalar: } \lambda f_c(p_c) > 1. \qquad (3.16) $$

We now proceed in calculating the phase space occupation number of radiated particles in the shells. The $c$ particles get accumulated within a thin shell whose thickness $L_b$ is computed using kinematics arguments in section 2.3. The associated number density of particles in the bath frame is

$$ n_{c,p} \approx \mathcal{N} \frac{\zeta(3)}{\pi^2} g_* T_n^3 \frac{R_c}{3L_b}, \qquad (3.17) $$

We deduce the occupation number $f_c$ (number of particles per de Broglie wavelength cell)

$$ f_c \approx \frac{n_{c,p}}{k_{\perp}^2 p_{c,b}} \approx \frac{\zeta(3)}{3\pi^2} g_* \mathcal{N} \frac{R_c T_n^3}{k_{\perp}^2 p_{c,b} L_b} \qquad (3.18) $$

In figure 6, we show in blue the regions where perturbativity breaks-down, using eqs. (3.14) and (3.15). In this plot we accounted for the full expressions for $\mathcal{N}$, $p_{c,b}$ and $k_{\perp}$ given in appendix A.1. For each quantity $X$, we included both the mean $\langle X \rangle$ and width $\sigma_X$ by performing the simple quadratic sum $\sqrt{\langle X \rangle^2 + \sigma_X^2}$.

We now provide analytical estimates of the occupation numbers in the simplified cases where $L_b \simeq \bar{L}_b$. For shells of particles generated at leading order (LO), we can use $L_b \simeq \bar{L}_b \simeq R_c(T_n/m_{c,h})^2$, $\mathcal{N} = 1$, $p_{c,b} \simeq m_{c,h}^2/T_n$ and $k_{\perp} \simeq T_n$ as shown in the first row of table 1 (with $X \equiv c$), leading to

$$ f_c(p_c) = \zeta(3)g_s/\pi^2, \qquad (\text{LO Shell}). \qquad (3.19) $$

Replacing the quantities $L_b \simeq \bar{L}_b \simeq R_c/(2\gamma_w^2)$, $\mathcal{N} \simeq (\alpha_D/\pi)L_m^2$, $p_{c,b} \simeq 1.6\gamma_w m_{c,h}/L_m^{1/2}$ and $k_{\perp} \simeq T_n$ in eqs. (2.31) and (2.35), (A.7), (A.11) and (A.13) for shells of reflected radiated particles and $\mathcal{N} \simeq (2\alpha_D/\pi)L_m L_E$, $p_{c,b} \simeq 0.4\gamma_w m_{c,h}/L_E^{1/2}$ and $k_{\perp} \simeq T_n$ for transmitted radiated particles, cf. second row of table 1, the occupation number of particles $c$ reads

$$ f_c(p_c) = g_D^2 \left( \frac{g_*}{100} \right) \left( \frac{\gamma_w T_n}{m_{c,h}} \right) \times \begin{cases} L_m^{5/2}/8, & (\text{Reflected Shell}), \\ L_E^{3/2}, & (\text{Transmitted Shell}), \end{cases} \qquad (3.20) $$

with $L_m \equiv \log(m_{c,h}/m_{c,s})$, and $L_E \equiv \log(E_a/m_{c,h})$. For shells of particles produced by Azatov-Vanvlasselaer mechanism [25, 26, 30, 59], in eq. (3.18) we plug $L_b \simeq \bar{L}_b \simeq R_c(T_n/M_X)^2$
to be conservative, $\mathcal{N} = \lambda^2 v_{\phi}^2 / 192\pi^2 M_X^2$, $p_{c,b} \simeq M_X^2 / T_n$ and $k_{\perp} \simeq T_n$ as shown in the last row of table 1, leading to

$$f_c(p_c) = \frac{\zeta(3)}{\pi^2} g_* \left(y \frac{v_{\phi}}{M_X}\right)^2, \quad (\text{Azatov-Vanvlasselaer mechanism}), \quad (3.21)$$

which is much smaller than 1, implying that phase-space saturation effects (both the perturbativity breakdown given by eqs. (3.16) and the Bose enhancement given by eq. (3.8)) are negligible.

# 4 Shell momentum loss

In this section we consider elastic 2 → 2 interactions between shell and bath particles and determine the conditions necessary for shells to free stream without the momentum losses induced by these interactions. We identify the possible impediments listed in section 4.1, and shade the associated regions of no free-streaming individually in figures 10, 11 and their envelope in figure 1.

## 4.1 Possible processes inducing momentum losses

Consider a shell produced on either side of the wall and propagating outwards from the nucleation site in the same direction of the wall. Particles in the shell will see an incoming flux of particles from the thermal bath. Interactions with these bath particles may change the momentum of the shell particles and therefore modify the overall properties or propagation of the shell. Also the bath particles may be affected by the shell and perhaps the flux of particles reaching the wall be suppressed or otherwise modified. We describe below the different scenarios arising from shell-bath momentum exchange and picture them in figure 7.

* **Reversal of the shell:** consider a shell travelling in front of the wall. In the bath frame, both wall and shell travel at close to the speed of light, but the shell is slightly faster than the wall. Interactions with the bath particles may lead to a change in momentum of the shell particles, insufficient for dissipation, but sufficient to slow the shell particles so that the latter are caught by the wall. The picture in the wall frame is the following: shell particles travel outward with momentum $p_{X,w}$. Typically one has $p_{X,w} \approx m_X$ in the wall frame for Bremstrahlung type production. Incoming bath particles have momentum $\approx \gamma_w T_n$. These interact with the shell particles and if the change of momentum of the latter in the wall frame is $\Delta p_{X,w} \approx p_{X,w}$ before the shells collide, then the shell particles are typically caught by the wall. More concretely, the condition to avoid the shell reversal is given by

$$\frac{1}{p_{X,w}} \left. \frac{dp_{X,w}}{dt} \right|_{\text{wall}} < \frac{\gamma_w}{R_c}, \quad (4.1)$$

where the $\gamma_w$ factor takes into account the shorter propagation distance before wall collision in the wall frame. Note in the wall frame the density of bath particles is also Lorentz boosted $\sim \gamma_w T_n^3$. We calculate the conditions for shell reversal in section 4.3.
◇ **Reversal of the bath:** if a relativistic shell is travelling in front of the wall, its interactions with the bath could in principle reverse the bath particles, and prevent them from reaching the wall. This can occur only if the shell carries more energy than the bath in the wall frame. We provide a general argument in section 4.2.2, however, that this is energetically impossible because the shell itself comes from bath particles.

◇ **Shell dissipation/metamorphosis:** in the extreme case the shell may be partially or fully dissipated. Consider particles travelling outward with momentum $p_X$ in the bath frame. Depending on the nature of the shell, we can have $p_X \approx \gamma_w v_\phi$ where $\gamma_w$ is the Lorentz factor of the wall, or $p_X \approx m_X^2/T_n$, see table 1. Then if interactions with the bath particles leads to changes in the shell particle momentum $\Delta p_X \approx p_X$ before the shells collide, then the shell will become dissipated back into the thermal bath, provided the bath carries sufficient energy to do so. To ensure free streaming, we thus require the momentum exchange to small over the distance of propagation, here taken to be the bubble size at collision,

$$ \frac{1}{p_X} \left. \frac{dp_X}{dt} \right|_{\text{bath}} < \frac{1}{R_c}. \qquad (4.2) $$

However, if the total energy in the shell particles is larger than the total energy in the bath particles, $E_{\text{shell total}} > E_{\text{bath total}}$, the latter cannot fully dissipate the former. Then a short path length for $\Delta p_X \approx p_X$ interactions, indicates the shell will transfer an $O(1)$ fraction of its momentum to bath particles, which themselves form part of an outward propagating shell, but now with modified properties. The number density, species, and momenta of shell particle is changed, so we dub this metamorphosis. The energy condition to have either dissipation or metamorphosis is determined in section 4.2.1, and the associated free-streaming conditions are calculated in section 4.4.

The rest of this section is structured as follows. To set the stage, we first derive the required energy conditions, for shell dissipation/metamorphosis and for shell/bath reversal (showing the bath cannot be reversed in the wall frame). We then go on to evaluate the reversal and dissipation/metamorphosis path lengths for the shell, for different choices of underlying interactions, in greater detail.

## 4.2 Total energy considerations

To better understand the end state of each the shell/bath interactions, it is instructive to consider the total energy/momentum of the shell and bath. The free energy difference between the true and false vacua leads to a pressure gradient which works on the bubble wall. The bath carries an associated pressure but no net momentum in the bath frame. Through the microphysical processes of bath-wall interactions, momentum of the wall is transferred to the resulting shell. This gives a net outward going momentum to the shell which in turn encounters the bath particles (either inside or outside the bubble, depending on the nature of the shell). Interactions between the shell and bath lead to further transfers of energy and momentum between the two when the above free streaming conditions are violated.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Shell-bath momentum loss mechanisms: interactions between particles of the shells and particles of the bath generate momentum transfer between the two systems. Clockwise starting from top-left, different scenarios occur according to whether in the wall frame the shell typical momentum is reversed (shell reversal), the bath momentum is reversed (bath reversal), in the bath frame the shell momentum is changed by $O(1)$ (shell metamorphosis), or is fully dissipated (shell dissipation), see section 4 for the details. The last condition, shell dissipation, requires the bath to have more energy than the shell in the bath frame, see eq. (4.5). The parameter space where shell reversal occurs is shown in figure 10. The one where shell metamorphosis/dissipation occurs is shown in figure 11. Bath reversal in the wall frame is never energetically feasible, see eq. (4.8).</figcaption></figure>

### 4.2.1 Shell dissipation or metamorphosis?

As previously mentioned, our shell dissipation/metamorphosis condition takes as a threshold an $O(1)$ change in the momenta of the shell particles, but this does not mean the complete erasure of the shell, only a significant change from its character in the limit of free streaming. Here we derive the required total energy of the bath in order for it to dissipate the shell completely back into a thermal state. We work in the bath frame. The total energy of bath particles in a Hubble volume, as seen in the bath frame, is given by

$$E_{\text{bath total,b}} = \rho_{\text{bath}} V_{\text{Hubble}} = \frac{g_* \pi^2 T_n^4}{30} \frac{4\pi}{3H^3}. \quad (4.3)$$

Similarly, the total energy of relativistic shell particles in a Hubble volume, in the bath frame, is given by

$$E_{\text{shell total,b}} = n_{\text{bath}} \mathcal{N} p_X V_{\text{Hubble}} = \frac{g_* \zeta(3) T_n^3}{\pi^2} \frac{4\pi \mathcal{N} p_X}{3H^3}. \quad (4.4)$$
The ratio of the two is thus

$$R_{E \text{ shell diss.}} \equiv \frac{E_{\text{shell total,b}}}{E_{\text{bath total,b}}} \approx \frac{N p_X}{3T_n} \qquad (4.5)$$

The simplified prefactor 3 in the denominator of the above expression, in reality ranges between 2.7 and 3.15, where then former corresponds to a bath made purely of bosons and the latter to a bath made purely of fermions. (This subtlety is not captured by our simplified assignment of the same $g_*$ for $\rho_{\text{bath}}$ and $n_{\text{bath}}$.) Clearly, if $R_{E \text{ shell diss.}} > 1$, then the shell cannot be completely dissipated by the bath back into a purely thermal state before shell/bubble collision. The above condition should therefore be checked if the path length for dissipation/metamorphosis is short compared to $R_c$, in order to gain a better handle on the end state of the process.

As a simple example, we can take the reflected radiated vector bosons, for which $\mathcal{N} \approx \alpha_D L_m^2 / 4\pi$, $p_X \approx \gamma_w \sqrt{\alpha_D} v_\phi$, and for $\gamma_w$ assume the vacuum dominated runaway regime. Then we find

$$R_{E \text{ shell diss.}} \approx \frac{\alpha_D^{3/2} M_{Pl} \log^2(v_\phi/T_n)}{\sqrt{48 c_{\text{vac}} \pi^{2/3} v_\phi \beta / H}} \qquad (4.6)$$

which for typical values implies $R_{E \text{ shell diss.}} \gg 1$. Thus if the shell dissipation free streaming condition is violated in this scenario, the shell will not be completely dissipated back into a thermal state, but rather be metamorphosed into a mélange of original shell and accelerated bath particles, before colliding with opposing shells/bubbles. Other scenarios of interest can easily be checked using table 1 and the relevant $\gamma_w$ factor. The more precise energy condition, eq. (4.5), will be plotted together with the path lengths in our summary plots below, distinguishing regimes of dissipation/metamorphosis regarding the bath frame momentum of the shell.

### 4.2.2 Shell or bath reversal?

We can also ask a similar question for shell or bath reversal in the wall frame. We are interested in whether the shell or bath particles change direction in the wall frame and so perform the calculation in said frame for convenience.<fnref n="6" /> The total energy of bath particles in the wall frame is

$$E_{\text{bath total w}} = \gamma_w E_{\text{bath total,b}} \qquad (4.7)$$

The total energy of the shell particles in the wall frame is

$$E_{\text{shell total w}} = \frac{E_{\text{shell total,b}}}{\gamma_w} \qquad (4.8)$$

If the ratio of shell to bath energies in the wall frame is smaller than one, then the shell particles do not have carry enough energy to reverse the bath particles in the wall frame. In the wall

<footnote n="6">⁶Working in the bath frame leads to additional complications, e.g. because the bath particles gain a net outward momentum during the reversal process, so it is not as easy as simply taking into account the energy change required for shell reversal in the bath frame, $\delta E_X = -\hat{t}/(4T) = p_X^2/(\gamma_w^2 T)$.</footnote>
frame the bath particles therefore reach the bubble without significant changes in momentum. Moreover, it is energetically allowed for the same bath particles to reverse the shell particles so the latter are pushed into the bubble. Actually, for consistency, we always have

$$R_{E \text{ shell/bath rev.}} \equiv \frac{E_{\text{shell total. w}}}{E_{\text{bath total w}}} = \frac{N p_X}{3 \gamma_w^2 T_n} < 1. \qquad (4.9)$$

This can be understood as follows. The wall can be treated as an inertial frame over the timescale of the particle transitions through it. Thus we can boost into the effectively time independent wall frame, for which energy conservation holds between the incoming particles (the bath), and their products, of which the shell is a subset. Indeed, shell particles produced earlier in the bubble expansion are produced at smaller $\gamma_w$, so their energy contribution is actually smaller than this estimate. Thus, in the wall frame, the total energy of the shell is smaller than the total energy of the bath. In other words $R_{E \text{ shell/bath rev.}} < 1$: out of the bath and the shell, it is only the latter which can be reversed in the wall frame. In other words, eq. (4.9) is satisfied in the entire parameter space of our interest.

Having derived the energy requirements, we now turn to calculating the relevant path lengths for shell reversal and shell dissipation/metamorphosis.

## 4.3 Shell reversal

### 4.3.1 Basic picture

We consider the particles created at the wall which are reflected back into the unbroken phase and calculate their reversal path length. Gauge bosons are of particular interest because they obtain an enhanced production rate compared to fermions and scalars. The gauge bosons have a typical momentum in the wall frame $p_{X,w} \approx m_{c,h} \approx g_D v_\phi$ and zero vacuum mass (thus their momentum in the bath frame is $p_X = \gamma_w p_{X,w} \approx \gamma_w g_D v_\phi$, as in table 1).

We now derive the conditions for shell particles (either gauge bosons or others) to be sent back to the wall via interactions with the bath particles. We also consider Compton scattering, which is of relevance when U(1) gauge bosons interact with charged fermions or scalars in the bath, and which gives a longer path length. The limits from Compton scattering are therefore the relevant ones in the absence of t-channel gauge boson exchange processes. Finally we consider a non-gauged PT in which a scalar boson is reflected back into the symmetric phase.

### 4.3.2 Simple estimates

We work in the wall frame. The basic picture is illustrated in figure 8. (The same situation, but limited to Møller scattering, has also been previously been considered in [46].) The shell particle momentum before and after scattering is denoted $p_1$ and $p_3$ respectively. The bath particle momentum is denoted $p_2$ or $p_4$. The four momenta are

$$p_1 = (p_{X,w}, \ 0, \ 0, \ p_{X,w}) \qquad (4.10a)$$

$$p_2 = (\gamma_w T_n, \ 0, \ 0, \ -\gamma_w T_n) \qquad (4.10b)$$

$$p_3 = (q, \ 0, \ q, \ 0) \qquad (4.10c)$$

$$p_4 = (p, \ 0, \ -p \sin \theta_{\text{th}}, \ -p \cos \theta_{\text{th}}), \qquad (4.10d)$$
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The scattering — as seen in the wall frame — leading to the threshold at which point the reflected particle will head back to the wall.</figcaption></figure>

where $\theta_{\text{th}}$ is the minimal deflection angle of bath particles which lead to the reversal of the momentum of reflected shell particles, i.e. the deflection angle for which $\vec{p}_3$ becomes parallel to the wall. Here we are ignoring particle masses and $\gamma_w$ is the Lorentz factor of the wall in the bath frame. There are three unknowns, $q, p, \theta_{\text{th}}$, and luckily for us three equations from energy and momentum conservation:

$$ \gamma_w T_n + p_{X,w} = q + p \qquad (4.11a) $$

$$ \gamma_w T_n - p_{X,w} = p \cos \theta_{\text{th}} \qquad (4.11b) $$

$$ q = p \sin \theta_{\text{th}}. \qquad (4.11c) $$

The solution to these equations is

$$ q = \frac{2p_{X,w}\gamma_w T_n}{p_{X,w} + \gamma_w T_n} \approx 2p_{X,w} \qquad (4.12a) $$

$$ p = \frac{p_{X,w}^2 + (\gamma_w T_n)^2}{p_{X,w} + \gamma_w T_n} \approx \gamma_w T_n \qquad (4.12b) $$

$$ \sin \theta_{\text{th}} = \frac{2p_{X,w}\gamma_w T_n}{p_{X,w}^2 + (\gamma_w T_n)^2} \approx \frac{2p_{X,w}}{\gamma_w T_n}, \qquad (4.12c) $$

where we have also given the approximate solution in the limit $\gamma_w T_n \gg p_{X,w}$. Clearly then a small deflection of the incoming bath particle is sufficient to reverse the reflected particle. The minimum momentum transfer squared needed to reverse the momentum of reflected shell particles is then

$$ -\hat{t} = -(p_1 - p_3)^2 = 2p_1 \cdot p_3 = 4p_{X,w}^2. \qquad (4.13) $$

Note this is far below the COM energy squared $\hat{s} = 2p_1 \cdot p_2 = 4\gamma_w T_n p_{X,w}$. The COM momentum is likewise found in the massless limit $p_{CM}^2 = \hat{s}/4$.

**Møller scattering.** We assume a spin-averaged matrix element inspired by t-channel gauge boson exchange, as in Møller scattering in regular fermionic (or scalar) QED, see figure 9-left, and consider the leading contribution in the $(\hat{t} \to 0, \hat{u} \to -\hat{s})$ limit

$$ |M|^2 \approx \frac{128\pi^2 \alpha_D^2 \hat{s}^2}{\hat{t}^2}, \qquad (4.14) $$
<figure><img src="image_11.png" /><figcaption>Figure 9. Møller (left) and Compton (right) scattering between shell and bath particles.</figcaption></figure>

where we have added a factor 2 to effectively also account for the ($\hat{u} \to 0$, $\hat{t} \to -\hat{s}$) limit. This gives us

$$ \frac{d\sigma}{d\hat{t}} = \frac{|\mathcal{M}|^2}{16\pi\hat{s}^2} \approx \frac{8\pi\alpha_D^2}{\hat{t}^2}. \qquad (4.15) $$

For other examples, such as Bhaba scattering, or gluon scattering in an SU(3), we find $\mathcal{O}(1)$ differences, so for practical purposes the same as above. Other choices would not give an IR enhancement, so the eventual path length before reversal is minimized by an assumption of $t$-channel gauge boson exchange. All relevant cross sections for the $2 \to 2$ momentum-loss processes, computed in this section 4, are reported synthetically in table 2. The effective cross section to have one scattering impart the necessary momentum exchange is then

$$ \sigma_{\text{eff}} = \int_{-4p_{\text{CM}}^2}^{-4p_{X,w}^2} d\hat{t} \frac{d\sigma}{d\hat{t}} \approx \frac{2\pi\alpha_D^2}{p_{X,w}^2}. \qquad (4.16) $$

The upper boundary is the minimal momentum exchange in order to reverse the momentum of reflected shell particles, eq. (4.13), and the lower boundary is the maximal momentum exchange available. We use this to find the rate of such scatterings in the wall frame,

$$ \Gamma_{\text{wall}} = n_{\text{bath,w}}\sigma_{\text{eff}}v_{\text{Møl}} \approx \frac{4\gamma_w\alpha_D^2 g_* \zeta(3) T_n^3}{\pi p_{X,w}^2}, \qquad (4.17) $$

where we have used the bath particle density in the wall frame

$$ n_{\text{bath,w}} \approx \frac{\gamma_w g_* \zeta(3) T_n^3}{\pi^2}, \qquad (4.18) $$

and $v_{\text{Møl}} \approx 2$. This result leads to an effective path length before reversal

$$ l_{\text{wall}} \approx \frac{\pi p_{X,w}^2}{4\gamma_w \alpha_D^2 g_* \zeta(3) T_n^3}. \qquad (4.19) $$

The path length in the bath frame is larger by a factor of $l_{\text{bath}} = \gamma_w l_{\text{wall}}$. We will eventually also study the effects of multiple soft scatters not taken into account in the above. But we first consider some other possibilities for the matrix element.

**Compton scattering in fermion QED.** In the massless limit, see figure 9-right, the polarization and spin averaged matrix element squared reads

$$ |\mathcal{M}|^2 \approx -32\pi^2 \alpha_D^2 \left( \frac{\hat{s}}{\hat{u}} + \frac{\hat{u}}{\hat{s}} \right). \qquad (4.20) $$
Note the divergence in the limit $\hat{u} = (p_1 - p_4)^2 \to 0$. In terms of the scattering angle, this divergence occurs when the gauge boson scatters back directly toward the wall. The divergence is cut off by the finite fermion mass in the propagator. Using $\hat{u} = 2m_f^2 - \hat{s} - \hat{t}$, we write the effective cross section as

$$ \sigma_{\text{eff}} \approx \frac{1}{16\pi\hat{s}^2} \int_{-\hat{s}}^{-m_f^2} d\hat{u} |\mathcal{M}|^2 \approx \frac{2\pi\alpha_D^2}{\hat{s}} \log\left(\frac{\hat{s}}{m_f^2}\right) \approx \frac{\pi\alpha_D^2}{2\gamma_w T_n p_{X,w}} \log\left(\frac{4\gamma_w p_{X,w} T_n}{\mu^2}\right), \quad (4.21) $$

where we used $\hat{s} = 4\gamma_w T_n p_{X,w}$. The quantity $m_f^2 = \mu^2$ is the IR cut-off of the fermion propagator due to thermal and finite-density correction, see last line of eq. (2.38). Note the suppression compared with eq. (4.16). The effective path length in the wall frame is then

$$ l_{\text{wall}} \approx \frac{\pi p_{X,w}}{\alpha_D^2 g_* \zeta(3) T_n^2 \log\left(\frac{4\gamma_w p_{X,w} T_n}{\mu^2}\right)}. \quad (4.22) $$

The path length in the bath frame is therefore $l_{\text{bath}} = \gamma_w l_{\text{wall}}$ which is larger by a factor $\sim \gamma_w T_n / p_{X,w}$ than the equivalent estimate for Møller scattering in eq. (4.19).

**Compton scattering in scalar QED.** In the massless limit, the polarization averaged matrix element squared reads

$$ |\mathcal{M}|^2 \approx 64\pi^2 \alpha_D^2, \quad (4.23) $$

and therefore

$$ \sigma_{\text{eff}} = \frac{|\mathcal{M}|^2}{16\pi\hat{s}} \approx \frac{4\pi\alpha_D^2}{\hat{s}}. \quad (4.24) $$

The effective path length in the wall frame is then

$$ l_{\text{wall}} \approx \frac{\pi p_{X,w}}{2\alpha_D^2 g_* \zeta(3) T_n^2}. \quad (4.25) $$

The path length in the bath frame is therefore $l_{\text{bath}} = \gamma_w l_{\text{bath}}$.

**Non-gauged phase transition.** We now consider a case with no gauge bosons present. It is instructive to consider renormalizable scalar self interactions and a yukawa interaction between the scalar and some fermion

$$ \mathcal{L} \supset \frac{\mu_{\phi}^2}{2} \phi^2 + \frac{\mu_3}{3!} \phi^3 + \frac{\lambda}{4!} \phi^4 + \frac{y}{\sqrt{2}} \phi \bar{f} f + m_f \bar{f} f + \text{H.c.} \quad (4.26) $$

In principle, the scalar can also interact with other scalars in the bath, but these effects can easily be taken into account in a similar way as what we do here for the self interactions. The momentum of the reflected radiated scalar is set by $p_{X,w} \approx m_{\phi,t}$, its mass in the broken phase.
For the $\phi\phi \to \phi\phi$ scattering, and after suppressing $O(1)$ factors, we have

$$|\mathcal{M}|^2 \approx \lambda^2 + \lambda\mu_3^2 \left[ \frac{1}{(\hat{t} - m_{\phi,f}^2)} + \frac{1}{(\hat{s} - m_{\phi,f}^2)} \right] + \mu_3^4 \left[ \frac{1}{(\hat{t} - m_{\phi,f}^2)} + \frac{1}{(\hat{s} - m_{\phi,f}^2)} \right]^2. \quad (4.27)$$

Here the first term comes from the quartic scalar interaction, the third term from the $t$-channel and $s$-channel scalar exchange amplitude squared, and the middle term is the interference between the first and third. The effective scalar mass in the false vacuum, including thermal effects, is denoted $m_{\phi,f}$. The leading three terms of the effective cross section reads

$$\sigma_{\text{eff}} \approx \frac{\lambda^2}{32\pi\hat{s}} + \frac{\lambda\mu_3^2}{32\pi\hat{s}^2} \log\left(\frac{\hat{s}}{4p_{X,w}^2}\right) + \frac{\mu_3^4}{32\pi\hat{s}^2 p_{X,w}^2}, \quad (4.28)$$

where we have assumed $p_{X,w} \gtrsim m_{\phi,f}$ (otherwise the final two terms are further suppressed by the finite propagator mass). A natural expectation is $\mu_3 \sim m_{\phi,t}$ (or smaller in case of an approximate $Z_2$ symmetry), from which it immediately follows that the $\lambda^2$ term typically dominates, unless $\lambda \lesssim m_{\phi,t}/(\gamma_w T_n)$. As the quartic coupling is also unavoidably present in front of the bubble (the cubic can be suppressed by symmetry reasons in the unbroken phase), we focus on this coupling for simplicity. Then the shell reversal path length in the wall frame is

$$l_{\text{wall}} \approx \frac{64\pi^3 p_{X,w}}{\zeta(3)g_*\lambda^2 T_n^2}, \quad (4.29)$$

which in the bath frame is $l_{\text{bath}} = \gamma_w l_{\text{wall}}$.

For $\phi f \to \phi f$ scattering we have $s$-channel $f$ mediation and $t$-channel $\phi$ mediation, which after dropping $O(1)$ factors, gives

$$|\mathcal{M}|^2 \approx y^4 + \frac{\mu_3 y^3 m_f}{(\hat{t} - m_{\phi,f}^2)} + \frac{\mu_3^2 y^2 \hat{t}}{(\hat{t} - m_{\phi,f}^2)^2}. \quad (4.30)$$

Here the first term comes from the $s$-channel amplitude squared, the second from the interference, and the third from the $t$-channel squared. We thus find

$$\sigma_{\text{eff}} \approx \frac{y^4}{16\pi\hat{s}} + \frac{\mu_3 m_f y^3}{16\pi\hat{s}^2} \log\left(\frac{\hat{s}}{4p_{X,w}^2}\right) + \frac{\mu_3^2 y^2}{16\pi\hat{s}^2} \log\left(\frac{\hat{s}}{4p_{X,w}^2}\right), \quad (4.31)$$

where $m_f$ is the fermion mass in the false vacuum. The first term will typically dominate, unless the Yukawa coupling is very small, $y < \mu^3 m_f / (\gamma_w T_n m_{\phi,t})$, or $y \lesssim \mu_3 / \sqrt{\gamma_w T_n m_{\phi,t}}$, in which case the second and/or third terms become important. For brevity, focusing on the large $y$ case, we have a reversal path length in the wall frame

$$l_{\text{wall}} \approx \frac{32\pi^3 p_{X,w}}{\zeta(3)g_* y^4 T_n^2}. \quad (4.32)$$

Translated to the bath frame the length is instead $l_{\text{bath}} = \gamma_w l_{\text{wall}}$, which is of course qualitatively similar to the scalar self interaction and the scalar Compton scattering path lengths.
### 4.3.3 Integral method

In section 4.3.2, we have only considered scattering processes which are “hard” enough to reverse of momentum of reflected shells. We now consider the possibility for the momentum of shell particles to be reversed under the effect of a large number of “soft” scattering processes. We again begin by working in the wall frame. From the above discussion, for reversal, we need to change the initial momentum (and energy) $p_1$ by an $O(1)$ factor in the wall frame. Using energy conservation, the change in the gauge boson momentum magnitude is

$$ \delta p_{X,w} = \delta E_X = -(E_2 - E_4). \qquad (4.33) $$

Our first task is to find $\delta p_{X,w}$ as a function of $\hat{t}$. In the COM frame the momenta are

$$ p'_{1} = (p_{\text{CM}}, \quad 0, \quad 0, \quad p_{\text{CM}}), \qquad (4.34a) $$

$$ p'_{2} = (p_{\text{CM}}, \quad 0, \quad 0, \quad -p_{\text{CM}}), \qquad (4.34b) $$

$$ p'_{3} = (p_{\text{CM}}, \quad 0, \quad p_{\text{CM}} \sin \theta_{\text{CM}}, \quad p_{\text{CM}} \cos \theta_{\text{CM}}) \qquad (4.34c) $$

$$ p'_{4} = (p_{\text{CM}}, \quad 0, \quad -p_{\text{CM}} \sin \theta_{\text{CM}}, \quad -p_{\text{CM}} \cos \theta_{\text{CM}}). \qquad (4.34d) $$

Here $\hat{s} = 4p_{\text{CM}}^2 = 4\gamma_{\text{w}}T_n p_{X,w}$. To go from the wall frame to the COM frame requires a relativistic boost, $v_{\text{boost}} \simeq 1$, in the direction of the wall with Lorentz factor

$$ \gamma_{\text{boost}} = \frac{p_{\text{CM}}}{2E_1} = \frac{E_2}{2p_{\text{CM}}} = \frac{1}{2} \sqrt{\frac{\gamma_{\text{w}}T_n}{p_{X,w}}}. \qquad (4.35) $$

Note we have the relation between the scattering angle in the COM frame and the momentum exchange

$$ \cos \theta_{\text{CM}} = 1 + \frac{\hat{t}}{2p_{\text{CM}}^2} \qquad (4.36) $$

Boosting from the COM frame back into the wall frame we have

$$ E_2 = \gamma_{\text{boost}}(E'_2 - v_{\text{boost}}p'_{2z}) = 2\gamma_{\text{boost}}p_{\text{CM}}, \qquad (4.37a) $$

$$ E_4 = \gamma_{\text{boost}}(E'_4 - v_{\text{boost}}p'_{4z}) = \gamma_{\text{boost}}p_{\text{CM}}(1 + \cos \theta_{\text{CM}}) = 2\gamma_{\text{boost}}p_{\text{CM}}\left(1 + \frac{\hat{t}}{4p_{\text{CM}}^2}\right). \qquad (4.37b) $$

Using these we find

$$ \delta p_{X,w} \simeq -\frac{\hat{t}}{4p_{X,w}}. \qquad (4.38) $$

We have to be careful with the above formula, as it captures the change in the magnitude of the momentum, but we can only reverse a particle once. Thus we should cut-off the weighting for $-\hat{t} \gtrsim 4p_{X,w}^2$, by making the replacement

$$ \delta p_{X,w} \rightarrow \text{Min} \left[ -\frac{\hat{t}}{4p_{X,w}}, \quad p_{X,w} \right], \qquad (4.39) $$

in regions of phase space where such hard scatterings occur. Then the rate to lose an $O(1)$ fraction of momentum in the wall frame is given by

$$ \frac{1}{p_{X,w}} \frac{dp_{X,w}}{dt} \simeq \frac{n_{\text{bath,w}} v_{\text{M0l}}}{p_{X,w}} \int_{-4p_{\text{CM}}^2}^{-\mu^2} d\hat{t} \frac{d\sigma}{d\hat{t}} \delta p_{X,w}. \qquad (4.40) $$
<figure><img src="image_12.png" /><figcaption>**Figure 10. Shell reversal:** regions in which the momentum of shells of reflected radiated bosons is reversed in the wall frame, due to Compton-type (orange, eq. (4.44)) or Møller-type (blue, eq. (4.43)) scatterings with bath particles. The symbol “V” denotes vector bosons while “ψ” denotes fermions and scalars. The subscripts “s” and “b” denote particles from the shell and from the bath, respectively. In the Abelian case, Møller-type scattering can only happen if there are fermions or scalars in the shells, which itself depends on whether $2 \to 2$ pair production processes are active or not, see eq. (5.5). The region for the scalar self interaction is similar to the Compton scattering region for $\lambda = g_D^2$. We consider reflected shell particles produced from Bremsstrahlung radiation, assuming abelian (left) or non-abelian (right) gauge interaction, see first two rows of table 1, assuming the amounts of supercooling $T_n/T_{eq} = 1$ (dashed) and $T_n/T_{eq} = 10^{-2}$ (solid). We used $g_b = g_{SM} + g_{emit}$ with $g_{SM} = 106.75$ and $g_{emit} = 10 = N_F$. Grey shaded areas as in figure 1.</figcaption></figure>

**Møller scattering.** We apply this formula to our Møller scattering cross section to find

$$ \frac{1}{p_{X,w}} \frac{dp_{X,w}}{dt} \approx n_{\text{bath,w}} v_{\text{Møl}} \left[ \int_{-4p_{X,w}^2}^{-\mu^2} d\hat{t} \frac{-2\pi\alpha_D^2}{p_{X,w}^2 \hat{t}} + \int_{-4p_{CM}^2}^{-4p_{X,w}^2} \frac{8\pi\alpha_D^2}{\hat{t}^2} d\hat{t} \right] \quad (4.41a) $$

$$ \approx \frac{4\gamma_w\alpha_D^2\zeta(3)g_*T_n^3}{\pi p_{X,w}^2} \left[ 1 + \log\left(\frac{4p_{X,w}^2}{\mu^2}\right) \right]. \quad (4.41b) $$

The streaming length before reversal in the wall frame is therefore

$$ l_{\text{wall}} \approx \frac{\pi p_{X,w}^2}{4\gamma_w \alpha_D^2 g_* \zeta(3) T_n^3 \left[ 1 + \log\left(\frac{4p_{X,w}^2}{\mu^2}\right) \right]} \quad (4.42) $$

This finally leads us to the streaming length in the bath frame $l_{\text{bath}} = \gamma_w l_{\text{wall}}$. Note up to the logarithmic suppression factor appearing in the denominator, this is the same estimate as using the simple $\sigma_{\text{eff}}$ method. Assuming a thermal mass cut-off $\mu^2$ as in eq. (2.38), the suppression is $\sim \log(v_{\phi}^2/T_n^2)$ in the Abelian case.

### Compton scattering in fermion and scalar QED.

Here the scatterings are dominantly hard, so the use of the integral method will not change the estimates of reversal path length compared to the effective cross section approach.

### Non-gauged PT.
Here the leading terms are also hard scattering dominated, unless $\lambda$ and/or $y$ are sufficiently small to suppress the first terms in the effective cross sections of
eqs. (4.28) and (4.31). So the integral method will not change the results derived above for the reversal path length. The (typically subdominant) $\propto 1/\bar{t}^2$ term in the $\phi\phi \to \phi\phi$ scattering receives a modest $\sim \log(m_{\phi,t}/m_{\phi,f})$ correction.

### 4.3.4 Summary of reversal path lengths

For Møller type scattering, i.e. t-channel gauge boson exchange, we can use the path length for reversal in the bath frame

$$l_{\text{bath}} \approx \frac{\pi p_{X,w}^2}{4\alpha_D^2 g_* \zeta(3) T_n^3 \log\left(\frac{4p_{X,w}^2}{\mu^2}\right)} = \frac{\pi p_X^2}{4\gamma_w^2 \alpha_D^2 g_* \zeta(3) T_n^3 \log\left(\frac{4p_X^2}{\gamma_w^2 \mu^2}\right)}. \quad (4.43)$$

For Compton scattering involving fermions

$$l_{\text{bath}} \approx \frac{\gamma_w \pi p_{X,w}}{\alpha_D^2 g_* \zeta(3) T_n^2 \log\left(\frac{4\gamma_w p_{X,w} T_n}{\mu^2}\right)} = \frac{\pi p_X}{\alpha_D^2 g_* \zeta(3) T_n^2 \log\left(\frac{4p_X T_n}{\mu^2}\right)}. \quad (4.44)$$

And for Compton scattering involving scalars

$$l_{\text{bath}} \approx \frac{\gamma_w \pi p_{X,w}}{2\alpha_D^2 g_* \zeta(3) T_n^2} = \frac{\pi p_X}{2\alpha_D^2 g_* \zeta(3) T_n^2}. \quad (4.45)$$

For non-gauged PTs the quartic interaction typically gives the shortest reversal path length although some care must be taken in this case with the precise parameter space and field content, using the methods outlined above. For fermionic interactions with sizable Yukawa couplings, the path length is similar to the quartic interaction, with the replacement $\lambda^2 \to y^4$. We thus write

$$l_{\text{bath}} \approx \frac{64\gamma_w\pi^3 p_{X,w}}{\zeta(3)g_*(\lambda^2 + 2y^4)T_n^2} = \frac{64\pi^3 p_X}{\zeta(3)g_*(\lambda^2 + 2y^4)T_n^2}, \quad (4.46)$$

In the above, $\mu^2$ is the gauge boson thermal mass from eq. (2.38), $p_{X,w}$ is the initial gauge boson momentum in the wall frame, and we reintroduce the momentum in the bath frame, $p_X = \gamma_w p_{X,w}$, to aid later comparison using table 1. Finally in order, say, to check whether the shells meet before particle reversal, one simply compares the above path lengths to the required propagation distances, typically $\sim R_c$ in the bath frame. In figure 10 we compare path lengths, for selected shells, to the bubble radius for shell reversal for a shell of radiated and reflected gauge bosons.

## 4.4 Shell dissipation/metamorphosis

### 4.4.1 Basic picture

We now consider the shell dissipation/metamorphosis path length. By this we mean any process which changes the momentum of shell particles in the bath frame by an $O(1)$ factor. Note this could still leave an expanding shell with significant, albeit altered, mean momenta and particle types and number densities. We shall make additional comments, in the context of specific examples, clarifying the two possibilities when relevant below. (Similar calculations to those below, in the case of non-gauged PTs, have been used in [26].)
<table><caption>**Table 2.** List of $2 \to 2$ processes accounted in this work for the estimation of the momentum exchange between the shell and the bath (shell reversal in section 4.3 and dissipation/metamorphosis in section 4.4) and the $2 \to 2$ inelastic between shell-shell and shell-bath entering the thermalization conditions in section 5.1. For $d\sigma/d\hat{t}$, we give the leading divergent behaviour, in the limit of massless particles. The greater the magnitude of the divergence as $\hat{t} \to 0$ or $\hat{u} \to 0$, the higher the efficiency of the corresponding process (depending also on the required momentum exchange). The symbol “V” denotes vector bosons, “ψ” denotes fermions, “φ” denotes scalars, and “X” denotes any charged particles. The subscripts “s” and “b” denote particles from the shell and from the bath, respectively.</caption><thead><tr><th colspan="2">Channel</th><th>dσ/dt̂</th><th colspan="2">Application</th></tr></thead><tbody><tr><td colspan="2">Møller scattering<br>X_s X_b → X_s X_b</td><td>8πα<sub>D</sub><sup>2</sup>/t̂<sup>2</sup><br>cf. eq. (4.15)</td><td colspan="2"></td></tr><tr><td rowspan="3">Compton<br>scattering</td><td>Fermion QED<br>V_s ψ<sub>b</sub> → V_s ψ<sub>b</sub></td><td>2πα<sub>D</sub><sup>2</sup>/ŝu<br>cf. eq. (4.20)</td><td rowspan="2">Shell reversal</td><td rowspan="2">Shell dissipation<br>and metamorphosis</td></tr><tr><td>Scalar QED<br>V_s φ<sub>b</sub> → V_s φ<sub>b</sub></td><td>4πα<sub>D</sub><sup>2</sup>/ŝ<sup>2</sup><br>cf. eq. (4.24)</td></tr><tr><td colspan="2">Non-gauged scalar<br>φ_s φ<sub>b</sub> → φ_s φ<sub>b</sub></td><td>λ<sup>2</sup>/32πŝ<sup>2</sup><br>cf. eq. (4.28)</td><td colspan="2"></td></tr><tr><td rowspan="2">Pair<br>production</td><td rowspan="2">Fermion QED<br>V_s V_s → ψ_s ψ_s<br>V_s V_b → ψ_b ψ_s</td><td>8πα<sub>D</sub><sup>2</sup>/t̂ŝ<br>cf. eq. (5.3)</td><td rowspan="2">Thermalization</td><td rowspan="2"></td></tr><tr><td>2πα<sub>D</sub><sup>2</sup>/ŝ<sup>2</sup><br>cf. eq. (5.1)</td></tr><tr><td colspan="2">Scalar QED<br>V_s V_s → φ_s φ_s<br>V_s V_b → φ_s φ_b</td><td></td><td colspan="2"></td></tr></tbody></table>

In our calculation of the shell reversal, we were interested whether a shell propagating in front of the bubble wall would remain there, or be sent back into the bubble. Thus the primary interest was for shells of the reflected radiated bosons of transition radiation. The shell dissipation/metamorphosis, however, is relevant not only for the shells considered in the reversal, but also for shells formed behind the bubble wall through the processes summarized in section 2.3.

### 4.4.2 Simple method

We work in the bath frame and consider $2 \to 2$ scattering between a shell and bath particle. We assume the shell particles are relativistic in the bath frame so that $E_1 \equiv E_X \simeq p_X$. The initial momenta are approximately

$$p_1 = (p_X, \ 0, \ 0, \ p_X) \qquad (4.47a)$$

$$p_2 = (T_n, \ 0, \ 0, \ -T_n). \qquad (4.47b)$$

The COM energy squared is $\hat{s} \simeq m_X^2 + 4p_X T_n$ and the COM momentum squared is $p_{CM}^2 \simeq 4p_X^2 T_n^2/\hat{s}$. To bring the bath particle energy from the bath to the COM frame requires a relativistic boost in the positive $z$-direction with $v_{boost} \simeq 1$ and Lorentz factor $\gamma_{boost} \simeq p_{CM}/2T_n \gg 1$. Consider now a scattering between the two particles in the COM frame with a scattering angle $\theta_{CM}$, corresponding to Mandelstam variable $\hat{t} = -2p_{CM}^2(1 - \cos \theta_{CM})$. In
the COM frame the four-momenta are

$$p'_{1} = (p_{\text{CM}}, 0, 0, p_{\text{CM}}) \qquad (4.48)$$

$$p'_{2} = (p_{\text{CM}}, 0, 0, -p_{\text{CM}}) \qquad (4.49)$$

$$p'_{3} = (p_{\text{CM}}, 0, p_{\text{CM}}s_{\theta_{\text{CM}}}, p_{\text{CM}}c_{\theta_{\text{CM}}}), \qquad (4.50)$$

$$p'_{4} = (p_{\text{CM}}, 0, -p_{\text{CM}}s_{\theta_{\text{CM}}}, -p_{\text{CM}}c_{\theta_{\text{CM}}}), \qquad (4.51)$$

The change in momentum for the shell particle can be found by boosting back into the bath frame, and is given by

$$\delta p_X = E_3 - E_1 = -\gamma_{\text{boost}} v_{\text{boost}} p_{\text{CM}} (1 - \cos \theta_{\text{CM}}) = \frac{\hat{t}}{4T_n}. \qquad (4.52)$$

Thus to achieve $\delta p_X \approx -p_X$ requires rather hard scattering $\hat{t} \approx -4T_n p_X$, of the order of $\hat{s}$.

**Møller scattering.** For t-channel gauge boson exchange processes, such as Møller scattering, we have as above,

$$|\mathcal{M}|^2 \approx 128\pi^2 \alpha_D^2 \frac{\hat{s}^2}{\hat{t}^2}. \qquad (4.53)$$

and thus

$$\frac{d\sigma}{d\hat{t}} = \frac{|\mathcal{M}|^2}{16\pi\hat{s}^2} \approx \frac{8\pi\alpha_D^2}{\hat{t}^2}. \qquad (4.54)$$

For the hard scattering processes we are interested in for shell dissipation, we use $\hat{t} \approx -4T_n p_X$, so we take an effective cross section

$$\sigma_{\text{eff}} \approx \frac{2\pi\alpha_D^2}{p_X T_n}. \qquad (4.55)$$

Note the above ignores a suppression factor because the required momentum exchange is close to $\hat{s}$, nevertheless, it should give a suitable estimate up to $\mathcal{O}(1)$ factors. Remembering that we are working in the bath frame, so that the bath number density is

$$n_{\text{bath}} \approx \frac{g_* \zeta(3) T_n^3}{\pi^2} \qquad (4.56)$$

we find the scattering rate

$$\Gamma = n_{\text{bath}}\sigma_{\text{eff}}v_{\text{Møl}} = \frac{4\alpha_D^2 g_* \zeta(3) T_n^2}{\pi p_X}. \qquad (4.57)$$

Hence, the dissipation length in the bath frame is

$$l_{\text{bath}} = \frac{\pi p_X}{4\alpha_D^2 g_* \zeta(3) T_n^2}. \qquad (4.58)$$

It is instructive to compare this to the equivalent calculation for the shell reversal path length, eq. (4.19), in the case of reflected radiated gauge bosons. Noting $p_X = \gamma_W p_{X,w} \approx \sqrt{\alpha_D} \gamma_W v_\phi$, we find the path length for dissipation is a factor $\gamma_W T_n / (\sqrt{\alpha_D} v_\phi)$ larger than for reversal. Note the
path length is indeed longer as kinematics requires $\gamma_W > \sqrt{\alpha_D}v_\phi/T_n$ for the gauge bosons to be produced. As the large momentum change required for dissipation/metamorphosis, $\delta p_X \approx p_X$ given in eq. (4.52), exceeds the momentum change required for reversal, $\delta p_X \approx p_X^2/(\gamma_W^2 T_n)$ given eqs. (4.38) and (4.13), the former would necessarily imply also the latter. So it is self consistent to have a dissipation/metamorphosis path length longer than or approximately coinciding with reversal.

**Compton scattering in fermion QED.** The matrix element is given by

$$ |\mathcal{M}|^2 \approx -32\pi^2\alpha_D^2 \left(\frac{\hat{s}}{\hat{u}} + \frac{\hat{u}}{\hat{s}}\right) \qquad (4.59) $$

Thus we have an effective cross-section

$$ \sigma_{\text{eff}} \approx \frac{1}{64\pi\hat{s}p_{\text{CM}}^2} \int_{-4p_{\text{CM}}^2}^{-m_f^2} d\hat{u} |\mathcal{M}|^2 \approx \frac{\pi\alpha_D^2}{2p_{\text{CM}}^2} \log\left(\frac{4p_{\text{CM}}^2}{m_f^2}\right) \approx \frac{\pi\alpha_D^2\hat{s}}{8p_X^2 T_n^2} \log\left(\frac{16p_X^2 T_n^2}{\mu^2\hat{s}}\right), \qquad (4.60) $$

where we have assumed an IR cutoff from the fermion thermal mass $m_f^2 \approx \mu^2$ with $\mu$ given by the second line of eq. (2.38), and used $p_{\text{CM}}^2 \approx 4p_X^2 T_n^2/\hat{s}$. Note that scattering precisely at the u-channel singularity would correspond to replacing a shell particle of one type (say a gauge boson) with a shell particle of another type (say a fermion). This changes the nature of the shell, but one may hesitate to label it as dissipation. Nevertheless, at somewhat more moderate $\hat{t}$, the size of the final state momenta are also significantly altered, so this remains a valid estimate up to logarithmic factors. The above leads to a dissipation length in the bath frame of

$$ l_{\text{bath}} \approx \frac{\pi p_X}{\alpha_D^2 g_* \zeta(3) T_n^2 \log\left(\frac{4p_X T_n}{\mu^2}\right)}. \qquad (4.61) $$

Again we compare with the reversal path length in the case of reflected radiated gauge bosons, eq. (4.22), and now find it coincides with the dissipation path length. This is due to the u-channel singularity, which means hard scatterings dominate up to the IR cutoff in the effective particle mass. Physically this makes sense provided we treat the derived lengths as approximate up to $O(1)$ factors.

**Compton scattering in scalar QED.** As before, for Compton scattering in massless scalar QED we have

$$ |\mathcal{M}|^2 \approx 64\pi^2\alpha_D^2, \qquad (4.62) $$

and therefore

$$ \sigma_{\text{eff}} \approx \frac{4\pi\alpha_D^2}{\hat{s}}. \qquad (4.63) $$

As in the fermion QED case, the scatterings of interest are dominantly hard, but there is now no $\hat{u} \to 0$ singularity which makes interpretation easier. Accordingly, the dissipation path length in the bath frame is

$$ l_{\text{bath}} \approx \frac{\pi p_X}{2\alpha_D^2 g_* \zeta(3) T_n^2}. \qquad (4.64) $$
Making the comparison to the reversal path length in the case of reflected radiated gauge bosons, eq. (4.25), we find the path lengths coincide. Thus the scatterings which reverse the shell in the wall frame, also change its constituent particle momenta by $O(1)$ factors in the scalar QED case.

**Non-gauged PT.** We again study this case by assuming interactions of the form in eq. (4.26). For the $\phi\phi \to \phi\phi$ scattering, using the approximate matrix element squared, eq. (4.27), we find

$$ \sigma_{\text{eff}} \approx \frac{1}{32\pi\hat{s}} \left( \lambda^2 + \frac{\lambda\mu_3^2}{\hat{s}} + \frac{\mu_3^4}{\hat{s}^2} \right). \qquad (4.65) $$

This gives a dissipation path length in the bath frame of

$$ l_{\text{bath}} \approx \frac{64\pi^3 p_X}{g_* \zeta(3) T_n^2 \left( \lambda^2 + \frac{\lambda\mu_3^2}{4T_n p_X} + \frac{\mu_3^4}{16T_n^2 p_X^2} \right)} \qquad (4.66) $$

In the case of a hard scattering dominated process (the leading term given by $\lambda^2$), the dissipation and reversal path lengths coincide.

Turning now to scattering with fermions, $\phi f \to \phi f$, using eq. (4.30), we have

$$ \sigma_{\text{eff}} \approx \frac{1}{16\pi\hat{s}} \left( y^4 + \frac{\mu_3 y^3 m_f}{\hat{s}} + \frac{\mu_3^2 y^2}{\hat{s}} \right). \qquad (4.67) $$

Note we do not integrate to small $-\hat{t}$, as we require large momentum shifts, and so do not find any logarithmic enhancements for the $t$-channel mediated process, as we did for the shell reversal). The effective cross section gives a bath frame dissipation path length

$$ l_{\text{bath}} \approx \frac{32\pi^3 p_X}{g_* \zeta(3) T_n^2 \left( y^4 + \frac{\mu_3 y^3 m_f}{4T_n p_X} + \frac{\mu_3^2 y^2}{4T_n p_X} \right)} \qquad (4.68) $$

Again the dissipation and reversal path lengths coincide for the hard scattering dominated process (governed by the $y^4$ term).

### 4.4.3 Integral method

For our $2 \to 2$ scattering, we see from eq. (4.52) that the momentum change of the shell particle in the bath frame is given by $\delta p_X = -\hat{t}/4T_n$. To take into account the possibility of a large number of soft scatterings adding up to give a momentum change of order $p_X$, we can obtain an estimate of the path length by using an integral, as in eq. (4.40). The only difference is that here we are working in the bath frame. Furthermore, we do not need to cut-off the change in momentum weighting, as hard scatterings correspond to $\delta p_X \approx p_X$, and not more. Accordingly, the momentum loss rate is given by

$$ \frac{1}{p_X} \frac{dp_X}{dt} \approx \frac{n_{\text{bath}} v_{\text{Møl}}}{p_X} \int_{-4p_{\text{CM}}^2}^{-\mu^2} d\hat{t} \frac{d\sigma}{d\hat{t}} \delta p_X \approx -\frac{n_{\text{bath}} v_{\text{Møl}}}{4p_X T_n} \int_{-4p_{\text{CM}}^2}^{-\mu^2} d\hat{t} \frac{d\sigma}{d\hat{t}} \hat{t}. \qquad (4.69) $$

**Møller scattering.** We apply the above to $t$-channel gauge boson exchange, using eq. (4.54), and obtain

$$ \frac{d \log p_X}{dt} \approx \frac{4\alpha_D^2 \zeta(3) g_* T_n^2}{\pi p_X} \log \left( \frac{4p_X T_n}{\mu^2} \right). \qquad (4.70) $$
Thus the dissipation path length in the bath frame is given by

$$l_{\text{bath}} \approx \frac{\pi p_X}{4\alpha_D^2 \zeta(3) g_* T_n^2 \log\left(\frac{4p_X T_n}{\mu^2}\right)}, \quad (4.71)$$

which is logarithmically suppressed compared to the simple estimate, eq. (4.42), and replaces it as our preferred approximation. In the case of reflected radiated gauge boson shells, the ratio of dissipation to reversal path lengths remains the same when using the simple of integrated estimates, up to differences in the logarithmic factor. That is, up to logarithmic corrections, the dissipation path length is a factor of $\gamma_W T_n / (\sqrt{\alpha_D} v_\phi)$ longer than the reversal path length.

**Compton scattering in scalar and fermion QED.** As for the shell reversal, these dissipation path lengths are determined by hard scatterings, so the simple estimates are unchanged by use of the integral method.

**Non-gauged PT.** In the $\phi\phi \to \phi\bar{\phi}$ process we have a $1/\hat{t}^2$ term which receives a logarithmic correction. We find

$$l_{\text{bath}} \approx \frac{64\pi^3 p_X}{g_* \zeta(3) T_n^2 \left[ \lambda^2 + \frac{\lambda \mu_3^2}{4T_n p_X} + \frac{\mu_3^4}{16T_n^2 p_X^2} \log\left(\frac{4T_n p_X}{\mu^2}\right) \right]}, \quad (4.72)$$

where $\mu^2 = m_{\phi,f}^2$ or $\mu^2 = m_{\phi,t}^2$ depending on whether the scattering is occurring in the symmetric or broken phases respectively. The path length for $\phi f \to \phi f$ does not receive any such correction.

### 4.4.4 Summary of dissipation/metamorphosis lengths

We now summarize the dissipation path lengths due to $2 \to 2$ interactions between shell and bath particles. When the shell interacts with the bath via t-channel gauge boson exchange the dissipation path length in the bath frame is approximately

$$l_{\text{bath}} \approx \frac{\pi p_X}{4\alpha_D^2 \zeta(3) g_* T_n^2 \log\left(\frac{4p_X T_n}{\mu^2}\right)}, \quad (4.73)$$

as in Møller scattering. For Compton scattering with fermions we obtain

$$l_{\text{bath}} \approx \frac{\pi p_X}{\alpha_D^2 g_* \zeta(3) T_n^2 \log\left(\frac{4p_X T_n}{\mu^2}\right)}, \quad (4.74)$$

and for Compton scattering with scalars we obtain

$$l_{\text{bath}} \approx \frac{\pi p_X}{2\alpha_D^2 g_* \zeta(3) T_n^2}. \quad (4.75)$$

For non-gauged phase transitions, the leading effect will typically be given by

$$l_{\text{bath}} \approx \frac{64\pi^3 p_X}{g_* \zeta(3) T_n^2 (\lambda^2 + 2y^4)}, \quad (4.76)$$

coming from the quartic scalar interactions or Yukawa interactions. Once a more precise field content and parameter space is defined, the techniques developed above can be of course
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Shell dissipation and metamorphosis: in the blue and orange regions, the shell particles undergo $O(1)$ changes to their momentum in the bath frame, due to Møller (see eq. (4.73)) and Compton (see eq. (4.74)) scattering interactions with the bath, respectively. (The region for the scalar self interaction is very similar to the Compton scattering one for $\lambda = g_D^2$.) On one side of the gray lines, in the bath frame there is more energy in the shell than in the bath suggesting that it is energetically possible for the momentum of bath particles to be reversed in the bath frame (shell metamorphosis). On the other side of the gray lines, the bath contains enough energy to dissipate the momentum of shell particles in the bath frame (shell dissipation). It is important to emphasize that the gray lines solely represent energetic considerations, and that they are relevant only when they fall within the blue and orange regions. Only within the shaded regions momentum loss is efficient, and momentum is lost via metamorphosis or dissipation depending on whether one is on one side or the other of the gray line. The symbol “*V*” denotes vector bosons while “*ψ*” denotes fermions and scalars. The subscripts “*s*” and “*b*” denote particles from the shell and from the bath, respectively. In the Abelian case, Møller-type scattering can only happen if there are fermions or scalars in the shells, which itself depends on whether 2 → 2 pair production processes are active or not, see eq. (5.5). We consider shell particles produced from LO interaction (top), Bremsstrahlung radiation, either transmitted (middle) or reflected (bottom), assuming abelian (left) or non-abelian (right) gauge interaction, see first two rows of table 1, assuming the amounts of supercooling $T_n/T_{eq} = 1$ (dashed) and $T_n/T_{eq} = 10^{-2}$ (solid). We used $g_b = g_{SM} + g_{emission}$ with $g_{SM} = 106.75$ and $g_{emission} = 10 = N_F$. We fixed the PT completion rate $\beta/H = 20$ to maximise the GW signal even though the results of this paper are only weakly dependent on this value. Grey and black shaded areas as in figure 1.</figcaption></figure>
used to check precisely the dominant contribution to the dissipation path length. We remind the reader that in these expressions, $p_X$ is the shell particle momentum as measured in the bath frame, as found in table 1, and $\hat{s} \simeq 4p_X T_n$. Parameter space in which the shell is dissipated/metamorphosised for LO shells and radiated gauge bosons are shown in figure 11. Also shown is the energy condition, eq. (4.5), delineating the regions in which there is or is not sufficient energy in the bath for complete shell dissipation.

# 5 Thermalization

While the thermal bath on its own is a system in thermal equilibrium, the particles in the shell are strongly boosted and strongly compressed. Interactions between the thermal bath and the shell, as well as interactions within the shell itself, give rise to out-of-equilibrium processes, which — if happening sufficiently fast — can lead to a new state of equilibrium, with the particles’ energies and densities different with respect to the free-streaming case.

In section 5.1, we will first consider inelastic $2 \to 2$ processes, or equivalently $2 \to 2$ annihilations, which change the type of particles present but not the overall number. This can occur via either shell-bath or internal shell-shell interactions. In sections 5.2–5.4, we then go on to the more challenging problem and consider $2 \to 3$ and $3 \to 2$ type interactions, which also change the total number of shell particles. The processes leading to thermalization are summarized in figure 12.

## 5.1 Processes changing individual but not total numbers

### 5.1.1 Inelastic $2 \to 2$ processes between the shell and bath

Let us first consider $2 \to 2$ annihilations where one initial particle belongs to the shell and the other one belongs to the bath. We anticipate that their inclusion does not result in stronger constraints on free-streaming, than the momentum-loss processes considered in section 4 (where we took into account $2 \to 2$ processes between the shell and the bath). We give reasoning for this below.
The momentum loss could, in principle, proceed in two qualitatively different manners:

◇ First of all, the addition of inelastic scatterings between the shell and the bath simply increases the number and type of interactions and therefore reduces the free streaming length. This effect can be captured by an $\mathcal{O}(1)$ enhancement of the total effective cross section for momentum loss — unless there is also an additional soft enhancement — which could lead to qualitatively more efficient scattering (e.g. giving Møller type scattering $|M|^2 \propto \hat{s}^2/\hat{t}^2$ when such a divergence was previously absent).

◇ Secondly, a two step process: efficient conversion via inelastic scattering to charged particles, followed by Møller type scattering. In particular, we saw in section 4.3, that for shell reversal, elastic scattering at $-\hat{t} \ll \hat{s}$ played an important role in the case of $t$-channel gauge boson exchange (due to the singular behaviour $|M|^2 \propto \hat{s}^2/\hat{t}^2$). Thus, for shell reversal, there is a strong sensitivity to the composition of the shell, i.e. whether it contains charged particles. We therefore need to check whether the consideration of inelastic processes between the shell and the bath (which could perhaps be enhanced at $-\hat{t} \ll \hat{s}$) could change the shell composition, and therefore also the type of scatterings leading to shell reversal, i.e. Møller instead of Compton scattering.

We now check these two possibilities. In both cases, the interactions between the shell and bath occur at high $\sqrt{s}$ compared to broken phase particle masses, so we can work in the massless limit. Let us consider a shell of gauge bosons encountering gauge bosons of the bath. The polarization averaged matrix element in scalar QED for $V + V \to \phi + \phi^*$ is given by

$$ |M|^2 \approx 32\pi^2\alpha_D^2, \qquad (5.1) $$

Resulting in a cross section

$$ \sigma v_{\text{Møl}} = \frac{|M|^2}{8\pi\hat{s}} = \frac{4\pi\alpha_D^2}{\hat{s}}, \qquad (5.2) $$

where we used that $v_{\text{Møl}} = 2$. The polarization averaged matrix element in fermion QED for $V + V \to f + \bar{f}$ is given by

$$ |M|^2 = 64\pi^2\alpha_D^2 \left( \frac{\hat{t}}{\hat{u}} + \frac{\hat{u}}{\hat{t}} \right). \qquad (5.3) $$

Giving an effective cross section

$$ \sigma v_{\text{Møl}} \approx \frac{1}{8\pi\hat{s}^2} \int_{-\hat{s}}^{-\mu^2} d\hat{u} |M|^2 = \frac{16\pi\alpha_D^2}{\hat{s}} \log \left( \frac{\hat{s}}{\text{Max}[\mu^2, -\hat{t}_{\text{min}}]} \right), \qquad (5.4) $$

where $\mu^2$ is the IR cutoff from the effective plasma mass, and $-\hat{t}_{\text{min}}$ is the minimum momentum exchange relevant for the process under consideration. We now comment on what enters the log factor depending on the process under consideration. We have $-\hat{t}_{\text{min}} \sim \hat{s} \approx 4T_n p_X$ for shell dissipation/metamorphosis from eq. (4.52), $-\hat{t}_{\text{min}} \sim 4p_{X,w}^2 \ll \hat{s}$ for shell reversal from eq. (4.13), and $-\hat{t}_{\text{min}} \sim \mu^2$ if we are simply interested in changing the composition of the shell. All relevant cross sections for the $2 \to 2$ annihilation processes, computed in this section 5.1, are reported synthetically in table 2.
Having found the cross sections, now let us consider the effect of the inelastic scatterings on shell momentum loss. For the scalar QED case, there is obviously no singular behaviour, and the inelastic cross section behaves the same way as the elastic one we previously saw in eqs. (4.24) and (4.63), and used in section 4. Also the shell composition is only altered at the same rate as momentum loss, so there is no change, up to $O(1)$ factors, due to the two-step process.

For the fermion QED case, we have $-\hat{t}, -\hat{u} \to 0$ singularities for $VV \to f\bar{f}$. But the same type of singularity has already been taken into account in eqs. (4.21) and (4.60) for elastic scattering in section 4. Thus there is no enhanced momentum loss for shells of abelian gauge bosons in either the scalar or fermion QED cases. Moreover, all the constrained regions in section 4 remain unchanged, up to $O(1)$ factors, when including inelastic processes between the shell and the bath. We therefore do not display these regions in our plots, to avoid clutter.

The above discussion has focused on shells of abelian gauge bosons. Let us argue that the same conclusions hold for other type of shells. First of all, for shells of non-abelian gauge bosons, charged fermions, or charged scalars — which all scatter via gauge boson exchange — we have already taken into account the soft enhancement in the matrix elements leading to momentum loss, $|M|^2 \propto \hat{s}^2/\hat{t}^2$, in our previous discussion in section 4 (in general both Møller and Compton scatterings take place for such shells). Secondly, if we have a non-gauged PT, there exist divergences of the type $|M|^2 \propto m^4/\hat{t}^2$ from eqs. (4.27), and $|M|^2 \propto m^2/\hat{t}$ from (4.30) ($m$ is some mass dimension one term), for which we have already accounted (amongst others, these are the most singular). But the addition of inelastic interactions will not give anything further enhanced, such as $|M|^2 \propto \hat{s}^2/\hat{t}^2$ (where one has $\hat{s} \gg m^2$), simply because by definition no gauge bosons are involved. Thus our conclusions of section 4 remain also for other types of shells.

### 5.1.2 Inelastic $2 \to 2$ processes within the shell

Let us now consider inelastic $2 \to 2$ processes within the shell, i.e. both initial annihilating particles belong to the same shell. If efficient, these processes alter the identities of particles making up the shells. For example shells made of radiated reflected gauge bosons would become made of gauge bosons plus symmetry-breaking scalars plus, if they exist in the spectrum, light fermions charged under the gauge group, with relative abundances set by the respective numbers of internal degrees of freedom.

Even if efficient, however, $2 \to 2$ processes starting with two shell particles do not alter the typical momenta of shells particles, $p_X$, nor their overall number densities, reported in table 1. Therefore, we will eventually include such $2 \to 2$ equilibrated shells within our definition of free-streaming shells, so that we will not display the areas where $2 \to 2$ processes are efficient in the summary figure 1. This is in contrast with the $3 \to 2$ processes studied below which alter both $p_X$ and the overall shell's number density.

The caveat regarding the $2 \to 2$ processes and momentum loss discussed in section 5.1.1, however, should be kept in mind. The two-step momentum loss process. In the case of radiated and reflected abelian gauge boson shells, we must also check whether $2 \to 2$ processes within the shell allow for shell momentum-loss via: (i) conversion to charged scalars/fermions, (ii) which then undergo more efficient Møller scattering with incoming charged bath particles.
Having made this specification, we calculate the rate of $2 \to 2$ processes within the shells, e.g. annihilation of vector bosons into scalar bosons or fermions. If the probability that such processes occur before collision is larger than unity, then we should expect the particles composition of the shell (but not the shell's overall number density) to reach relative chemical equilibrium. In the opposite case,

$$P_{2 \to 2} = n_{X,b} \sigma_{2 \to 2} v_{M\o l} R_c \ll 1, \qquad (5.5)$$

the composition of shell remains conserved during its propagation. The quantity $n_{X,b}$ is the shell number density given in eq. (2.24) and the bubble radius, $R_c$, is the effective distance over which the shell particles propagate before collision. For definiteness, we consider the process of vector bosons annihilation into scalar bosons or fermions. Its cross-section $\sigma$ multiplied by the Møller velocity $v_{M\o l}$, in the bath frame, can be approximated with naive-dimensional analysis as

$$\sigma_{2 \to 2} v_{M\o l} \approx \frac{g_D^4}{16\pi} \frac{1}{p_X^2}. \qquad (5.6)$$

For example, eq. (5.6) is indeed the result one obtains by considering the annihilation of vector bosons into a pair of scalars driving the PT.

In figure 13, we show in blue the regions where inelastic $2 \to 2$ shell-shell (labelled ‘ss’) interactions become efficient. As shells still free stream even if such processes are efficient, as we have argued in the beginning of this section, we leave their more complete study to future work.

We can now also address the two-step shell momentum-loss process with the initial annihilation taking place within the shell. Let us first consider shell reversal. Comparing the ‘ss’ vector boson annihilation region of figure 13, with the Compton scattering region of figure 10, it turns out there are regions in which annihilation is efficient, but where reversal via Compton scattering is not. Therefore — provided the annihilation into charged scalars or fermions is kinematically allowed — there exist regions of parameter space in which radiated and reflected shells of gauge bosons can convert into charged particles, which are then efficiently reversed via Møller scattering back to the wall. Thus, when such annihilations are efficient, then processes of shell reversal are important also for abelian reflected shells, and we shade the associated region both in figure 10 and in figure 1. We obtain those regions by requiring that both the annihilation $VV \to \phi\phi^*$ (with the cross section given in eq. (5.6)) and the momentum loss induced by $\phi\phi \to \phi\phi$ (where one initial $\phi$ belongs to the shell and one to the bath, and where we use $d\sigma/d\hat{t}$ from eq. (4.15)) are efficient at the same time. Coming finally to shell dissipation and metamorphosis, the one implied by the two-step process leads to no-free-streaming regions which overlap with those derived in section 4.4, as we show in figure 11. This is because the shell dissipation/metamorphosis requires significant momentum exchange, $-\hat{t}_{min} \sim \hat{s}$, so differences between Møller and Compton scattering are negligible.

## 5.2 Number changing processes: setup

In the rest of this section we consider those processes that involve number changing-interactions, see figure 12. To leading order in the couplings the relevant number-changing
<figure><img src="image_14.png" /><figcaption>**Figure 12.** *Thermalization.* LEFT: 3-body interactions between particles of the bath and of the shells can be efficient enough to thermalize the shell with the bath. Red and blue circles represent incoming bath and shell particles respectively, while empty legs represent outgoing particles. We denote the incoming particles originating from the shell and the bath by ‘s’ and ‘b’, respectively. The notations ‘ssb’ and ‘bss’ describe the same 3-body process — two shell particles interacting with one bath particle — given by the same cross-section. However, they are distinguished by their interaction probabilities: ‘ssb’ signifies the interaction probability from the perspective of a shell particle, whereas ‘bss’ reflects this probability from the standpoint of a bath particle. The ratio of interaction probabilities is given by $P_{ssb}/P_{bss} \approx n_b L_{scat}^s/(n_s L_{scat}^b)$, where $n_i$ denotes the number density, and $L_{scat}^i$ represents the effective interaction distance before bubble collision, as detailed in table 3. A similar distinction applies to the ‘sbb’ and ‘bbs’ configurations with $P_{sbb}/P_{bbs} \approx n_b L_{scat}^s/(n_s L_{scat}^b)$. RIGHT: inelastic 2-body and 3-body interactions between particles within the shells can be efficient enough to thermalize the shell. See section 5 for the details.</figcaption></figure>

processes are $2 \to 3$ and $3 \to 2$. At the initial stages of the evolution of a shell, it is enough to take into account only $3 \to 2$ processes, because $2 \to 3$ ones are suppressed by the final phase space. In the case of dense shells, the rate of $3 \to 2$ processes is further enhanced over the one for $2 \to 3$ processes due to an extra factor of the large shell’s number densities. Therefore, for the purpose of determining the regions where shells free stream, it is enough to compute the effect of $3 \to 2$ interactions. $2 \to 3$ interactions would have to be taken into account when shells approach some equilibration.

In case phase-space saturation is negligible (i.e. $1 \pm f \simeq 1$ with $f$ the occupation number of the particle of interest, see section 3), we can then write the effect of number changing interactions on the number density of particles $i$ as

$$ \frac{dn_i}{dt} \approx \Gamma_{3 \to 2}^i n_i, \qquad (5.7) $$

where the interaction rate is given by

$$ \Gamma_{3 \to 2}^i \approx \frac{1}{256\pi^2} \frac{n_j n_k}{E_1 E_2 E_5} \times \int d\Omega |M|^2, \qquad (5.8) $$

where $E_\ell$ is the energy of particle $\ell$ and $n_\ell$ is its number density, and where we have used the labeling of particles and momenta

$$ 125 \to 34 \qquad (5.9) $$

so $i, j, k \in \{1, 2, 5\}$ and $i \neq j \neq k \neq i$. We further specify that we label with 1 and 2 the initial state particles which belong to the same population. More explicitly 1 and 2 can be either
<table><caption>**Table 3.** Typical values of a particle's energy, number density, and effective distance over which it undergoes 3 → 2 interactions, depending on whether the particle belongs to a shell or to the bath. All quantities are given in the frame of the thermal bath. $g_b$ counts the bath degrees of freedom that participate in the scatterings of interest, $g_{emit}$ those that give rise to the shell of interest. See table 1 for the values of $p_x$, $N$ and $L_b$ for the various shells.</caption><thead><tr><th><i>i</i></th><th><i>E</i><sub><i>i</i></sub></th><th><i>n</i><sub><i>i</i></sub></th><th><i>L</i><sup><i>i</i></sup><sub>scat</sub></th></tr></thead><tbody><tr><td>bath particle</td><td>3<i>T</i><sub><i>n</i></sub></td><td><i>g</i><sub><i>b</i></sub> <sup>ζ(<i>3</i>)</sup>/<sub>π<sup>2</sup></sub><i>T</i><sub><i>n</i></sub><sup>3</sup></td><td><i>L</i><sub><i>b</i></sub></td></tr><tr><td>shell particle</td><td><i>p</i><sub><i>x</i></sub></td><td><i>N</i> <sup><i>R</i><sub><i>c</i></sub></sup>/<sub>3<i>L</i><sub><i>b</i></sub></sub> × <i>g</i><sub>emit</sub> <sup>ζ(<i>3</i>)</sup>/<sub>π<sup>2</sup></sub><i>T</i><sub><i>n</i></sub><sup>3</sup></td><td><i>R</i><sub><i>c</i></sub> ≃ (<i>π</i>)<sup>1/3</sup>β<sup>−1</sup></td></tr></tbody></table>

two shell particles, in which case 5 can be either shell or bath, or two bath particles, in which case 5 can only be a shell particle because usual number changing interactions within the bath alone are not the subject of our study. For example, if we are interested in the evolution of bath particles due to scatterings with two shell particles, then in eqs. (5.7) and (5.8) $i = 5, j = 1, k = 2$. $\int d\Omega |M|^2$ is the standard expression of the final two-particle phase-space integration of the spin-averaged squared matrix element $|M|^2$ in the center-of-mass frame. The total probability of a particle $i$ undergoing a 3 → 2 interaction before walls collide is

$$P_{3 \to 2}^i \approx \Gamma_{3 \to 2}^i \times L_{\text{scat}}^i, \qquad (5.10)$$

where $\Gamma_{3 \to 2}$ is the interaction rate and $L_{scat}^i$ is the effective distance (or equivalently time, because particles are ultra-relativistic) over which 3 → 2 interactions are possible for the particle $i$ under consideration. The condition

$$P_{3 \to 2} \ll 1, \qquad (5.11)$$

will then imply that number changing interactions, between the bath and the shell or within the shell itself, do not affect the propagation of shells nor the evolution of the bath.

## 5.3 Computation of 3 → 2 rates

### 5.3.1 Energies, densities and scattering lengths

A particle $i$ from the bath can only interact with shell particles when it traverses the shell. Thus, $L_{scat}^i$ represents the effective thickness of the shell, corresponding to $L_b$ listed in table 1. On the other hand, a particle $i$ originating from the shell has the potential to interact with either two other shell particles, one shell and one bath particle, or two bath particles throughout its entire journey. Hence, in this context, $L_{scat}^i$ is equivalent to the bubble radius.

Having quantified $L_{scat}^i$ in eq. (5.10), we now turn to the quantities that enter $\Gamma_{3 \to 2}^i$ in eq. (5.8). The energy $E_i$ and the number density $n_i$ of initial bath particles are of course the thermal ones. The energies and densities of initial shell particles can both be read off table 1, where the density is obtained by multiplying the bath one by $N$. A summary of these results is reported in table 3. We are left with the computation of $\int d\Omega |M|^2$, which we address next.
<table><caption>**Table 4.** U(1) and SU(N) columns. Parametric dependence of $\int d\Omega|M|^2$ in the scalar products, the gauge coupling $g_D$ and N. Only the leading terms in $g_D$, N and the scalar products are kept, where in the ssb and bbs cases $\mu^2 \ll p_1 \cdot p_2 \ll p_1 \cdot p_5 \simeq p_2 \cdot p_5$ (we remind we defined 1 and 2 to belong to the same population, i.e. to a shell 's' or to the bath 'b'), while in the sss case there is no hierarchy among the different scalar products and $\mu^2 \ll p_i \cdot p_j$. The reported values assume that shells are made of vector bosons, and that they couple to fermions. $\int d\Omega|M|^2_{125\to34}$ have the same or a more suppressed scaling in case i) there are no fermions charged under the gauge group, or ii) shells are made of fermions and/or scalars but not vectors (as e.g. reflected scalar shells in a global PT or shells from the Azatov-Vanvlasselaer mechanism). See appendix C for results that keep O(1) factors, subleading terms and the identities of initial and final states. **p<sub>1</sub>·p<sub>2</sub> and p<sub>1,2</sub>·p<sub>5</sub> columns.** Values of the scalar products of interest, where $k_\perp^2$ and $p_x$ depend on the shell of interest. $p_x$ can be read off table 1, while $k_\perp^2 = O(g_D^2 v_\phi^2 + T_n^2)$ for both reflected and transmitted radiated vectors, $k_\perp^2 = O(T_n^2)$ for shells of particles getting a mass, see eq. (5.12) and text for more details. We add the IR cutoff $\mu^2$ to all scalar products.</caption><thead><tr><th>∫ dΩ|M|<sup>2</sup><sub>125→34</sub></th><th>U(1)</th><th>SU(N)</th><th>p<sub>1</sub>·p<sub>2</sub></th><th>p<sub>1,2</sub>·p<sub>5</sub></th></tr></thead><tbody><tr><td>125 = sss</td><td>4πg<sup>6</sup><sub>D</sub> <sup>1</sup>/<sub>p<sub>i</sub>·p<sub>j</sub></sub> ln <sup>p<sub>i</sub>·p<sub>j</sub></sup>/<sub>μ<sup>2</sup></sub></td><td>4πg<sup>6</sup><sub>D</sub> <sup>1</sup>/<sub>N<sup>2</sup></sub> ln <sup>p<sub>i</sub>·p<sub>j</sub></sup>/<sub>μ<sup>2</sup></sub></td><td>k<sub>⊥</sub><sup>2</sup> + μ<sup>2</sup></td><td></td></tr><tr><td>ssb</td><td>4πg<sup>6</sup><sub>D</sub> <sup>1</sup>/<sub>p<sub>1</sub>·p<sub>2</sub></sub> ln <sup>p<sub>1</sub>·p<sub>2</sub></sup>/<sub>μ<sup>2</sup></sub></td><td>4πg<sup>6</sup><sub>D</sub> <sup>p<sub>1</sub>·p<sub>5</sub></sup>/<sub>(p<sub>1</sub>·p<sub>2</sub>)<sup>2</sup></sub> ln <sup>p<sub>1</sub>·p<sub>2</sub></sup>/<sub>μ<sup>2</sup></sub></td><td>k<sub>⊥</sub><sup>2</sup> + μ<sup>2</sup></td><td></td></tr><tr><td>bbs</td><td>4πg<sup>6</sup><sub>D</sub> <sup>1</sup>/<sub>p<sub>1</sub>·p<sub>2</sub></sub> ln <sup>p<sub>1</sub>·p<sub>5</sub></sup>/<sub>μ<sup>2</sup></sub></td><td>4πg<sup>6</sup><sub>D</sub> N<sup>2</sup> <sup>p<sub>1</sub>·p<sub>5</sub></sup>/<sub>(p<sub>1</sub>·p<sub>2</sub>)<sup>2</sup></sub> ln <sup>p<sub>1</sub>·p<sub>2</sub></sup>/<sub>μ<sup>2</sup></sub></td><td>T<sub>n</sub>p<sub>x</sub> + μ<sup>2</sup></td><td></td></tr></tbody></table>

### 5.3.2 Integrated amplitudes

We compute $\int d\Omega|M|^2$ for all possible $3 \to 2$ processes involving in the initial state at least one gauge boson V, scattering with other V's and/or fermions and/or scalars charged under the gauge group.

We perform these computations for the cases of both an Abelian and a non-Abelian gauge symmetry, which for concreteness we take as U(1) and SU(N), because the self-interactions of the vectors eventually lead to important differences between the two.

We compute the matrix elements with the help of FeynRules [70], FeynArts [71], and FeynCalc [72–74]. We refer the reader to appendix B for a detailed explanation of our parametrisation and of our integration of the matrix elements over the final phase space phase, and to appendix C for detailed results of our computation of $\int d\Omega|M|^2$ in terms of scalar products of the 4-momenta in the center-of-mass frame $p_1+p_2+p_5 = p_3+p_4$, including for definiteness only the terms that are leading order in $p_i \cdot p_j/\mu^2$. Indeed, we anticipate from section 5.3.3 that one has the hierarchy $\mu^2 \ll p_1 \cdot p_2 \ll p_1 \cdot p_5 \simeq p_2 \cdot p_5$ in case at least one initial particle belongs to the bath, while all scalar products are of the same order in case all initial particles belong to the same shell. In the same appendix C, we also report simple estimates of the upper limits one can expect on $\int d\Omega|M|^2$, as a check of our results and as an orientation to obtain analogous ones in the future. For convenience of the reader, we report in table 4 the parametric dependence, of the leading terms in $\int d\Omega|M|^2$, on the scalar products, the gauge coupling $g_D$ and, for SU(N), on N.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Shell thermalization: in the coloured regions (on the left of the dashed lines) the $2 \to 2$ and $3 \to 2$ interactions specified in the insets (with “s” and “b” denoting incoming particles from shells and bath respectively), see eqs. (5.5) and (5.11), become efficient and shells do not free stream, for $T_n/T_{eq} = 1(10^{-2})$. We show shell particles produced from LO interaction (top), Bremsstrahlung radiation, either transmitted (middle) or reflected (bottom), assuming abelian (left) or non-abelian (right) gauge interaction, see first two rows of table 1. We used $g_b = g_{SM} + g_{emit}$ with $g_{SM} = 106.75$ and $g_{emit} = 10 = N_F$. Grey and black shaded areas as in figure 1.</figcaption></figure>
### 5.3.3 Scalar products

Depending on the identity of the initial-state scatterers, the scalar products have a different dependence on the parameters of the PT and of the theory, implying different results for the probability that the associated 3 → 2 interaction happens before collision, which is the question of our interest. We therefore now turn to derive expressions for the scalar products in terms of the parameters of the theory.

**Shell-shell-shell.** Let us start by the simple case where the 3 initial particles all belong to the shell. Despite in the bath frame the energies of the three particles are very large, their scalar products do not track of course their frame-dependent energies, but rather the typical spread of their momenta, $p_i \cdot p_j \simeq k_\perp^2$. This spread depends on the shell of interest. We use

$$k_{\perp}^2 = \begin{cases} T_n^2 & \text{particles acquiring mass} \\ \langle k_{\perp, \mathcal{R}}^2 \rangle + \sigma_{k_{\perp, \mathcal{R}}}^2 & \text{reflected radiated vectors} \\ \langle k_{\perp, \mathcal{T}}^2 \rangle + \sigma_{k_{\perp, \mathcal{T}}}^2 & \text{transmitted radiated vectors} \end{cases} \quad (5.12)$$

where $\langle k_{\perp, \mathcal{R}, \mathcal{T}} \rangle$ and $\sigma_{k_{\perp, \mathcal{R}, \mathcal{T}}}$ are derived in appendix A.1 and can be read off, respectively, in eqs. (A.12) and (A.13).

**Shell-shell-bath and bath-bath-shell.**

Here $p_1$ and $p_2$ refer to particles that either both belong to the shell or both belong to the bath. In particular for shell-shell-bath interactions we have $p_1$ and $p_2$ for the shell particles and $p_5$ for the bath particle; for bath-bath-shell interactions we have $p_1$ and $p_2$ for the bath particles and $p_5$ for the shell particle. Again, in some frame $p_1$ and $p_2$ are nearly parallel (e.g. if they denote bath particles and one works in the wall frame), but their scalar product does not feel the large energies that they can have in that frame. One then has the usual $p_1 \cdot p_2 \simeq T_n^2$ for bath particles, and $p_1 \cdot p_2 \simeq k_\perp^2$ for shell particles, whose origin has already been explained in the shell-shell-shell case. One instead has $p_{1,2} \cdot p_5 \simeq p_\chi T_n$, where $p_\chi$ is the momentum of the given shell in the bath and can be read off table 1. We report the resulting values for the various scalar products in table 4.

In this paper we are not interested in the impact of 3 → 2 interactions on the evolution of the bath and the shell, but only in determining the region where they are relevant or not. We still find it worth to note that, by a careful investigation (see appendix B) of just the scattering kinematics, it is straightforward to show that the two final state particles are of one type shell and one type bath, i.e. one final particle has the typical energy (up to $\mathcal{O}(1)$-factors) of a bath one and is traveling inside the bubble, and the other particle has the typical energy of a shell one and is traveling along with other shell particles.

## 5.4 Results

Using eq. (5.10) for the probability $P_{3 \to 2}^i$ that a particle $i$ undergoes a 3 → 2 interaction, with the rate from eq. (5.8) and the other inputs from tables 3, 4, the condition $P_{3 \to 2}^i < 1$ is not satisfied in the colored regions in figure 13, where we shade in different colors regions where (in parenthesis the momentum to assign to each particle, according to our definitions)

◇ $i(p_1), j(p_2), k(p_5)$ are all shell particles;
◇ $i(p_1), j(p_2)$ are shell particles and $k(p_5)$ is a bath one;

◇ $i(p_5)$ is a shell particles and $j(p_1), k(p_2)$ are bath ones;

◇ $i(p_5)$ is a bath particles and $j(p_1), k(p_2)$ are shell ones;

◇ $i(p_1), j(p_2)$ are bath particles and $k(p_5)$ is a shell one.

Among the 5 cases above, the one where the three initial particles all belong to the shell dominates the shaded regions for radiated shells (both transmitted and reflected) in the abelian case. This can be understood with the fact that, in those cases, $n_{\text{bath}}/n_{\text{ej}} L_{\text{eff}}/R_c \ll 1$. The shaded regions for radiated shells, in the non-abelian case, instead all lie within the parameter space where our computation cannot be trusted because the plasma mass is larger than the vacuum one, see section 3.1 and figure 6.

For completeness, in figure 13 we also shade the regions where $2 \to 2$ processes within shells are efficient, discussed in section 5.1: while they do not alter the total shell's number density and typical particle momentum and so do not make shells depart from free-streaming, they could affect the identities of the particles making up the shell.

Overall, thermalization processes are the strongest ones in determining the regions where radiated shells from an abelian gauge theory do not free stream, as well for LO shells from non-abelian theories at large $v_\phi$ values, see figure 1. An example consequence of these findings is that the first processes to consider, in order to determine the evolution of radiated and transmitted shells at PTs associated with an abelian gauge group, would be $3 \to 2$ ones all involving shell particles in the initial state. We leave this interesting direction for future work.

# 6 Interaction with bubble walls

Field dynamics after wall-wall collision has been investigated by a number of authors, starting from [75]. The study of particle production from wall-wall collisions was first conducted in [76], followed by subsequent applications and refinements in [21, 27, 28, 77–83]. In section 6.1 we characterize wall-wall collisions following [21] for the analytical expressions and [81] for the numerical treatment. In section 6.2 we use these existing modelings to derive novel conclusions about the propagation of shells.

## 6.1 Wall-wall collision

We model two incoming bubble walls as two Heaviside functions propagating in opposing direction

$$ \phi(z, t < 0) = v_{\phi} + \frac{v_{\phi}}{2} \left[ \tanh \left( \gamma_w \frac{z+t}{L_0} \right) - \tanh \left( \gamma_w \frac{z-t}{L_0} \right) \right], \quad (6.1) $$

where $L_0$ is the wall thickness in the wall frame. The evolution of the scalar field profile during bubble collision can be determined analytically by integrating the Klein-Gordon equation

$$ (\partial_t^2 - \partial_z^2)\phi(z, t) + V'(\phi) = 0, \quad (6.2) $$
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Shell-wall interactions: after passing through one another, walls 1 and 2 progressively transform their gradient energy into an oscillating scalar field condensate (labelled as “wall collision remnants”) that develops between them as they continue to move apart. Shell 1, associated with wall 1, can interact both with the gradient energy of wall 2 and with the wall collision remnants. Similarly, shell 2, accompanying wall 2, can interact both with the gradient energy of wall 1 and the wall collision remnants. See section 6 for details of shell-wall interactions. Of course, shells 1 and 2 can also interact with each other, see [20].</figcaption></figure>

with eq. (6.1) as boundary condition. Taking $V(\phi) = m_{\phi}^2(\phi - v_{\phi})^2/2$, one obtains [21]

$$ \phi(z, t > 0) = v_{\phi} \left[ 1 + \frac{L_0}{\gamma_w} \int_0^{\infty} dp_z \frac{p_z}{\sqrt{p_z^2 + m_{\phi}^2}} \frac{\cos(p_z z)}{\sinh\left(\frac{\pi L_0 p_z}{2\gamma_w}\right)} \sin\left(\sqrt{p_z^2 + m_{\phi}^2} t\right) \right]. \quad (6.3) $$

We show the analytical solution for the scalar profile in figure 15, which we find to be in good agreement with the profile obtained from numerically integrating the equation of motion in eq. (6.2) following the recipe of [81], and shown in figure 16.

After walls pass through each other, their energy gets released into oscillating modes of the scalar field as the walls continue to propagate. The energy stored in the bubble wall in the bath frame at the time of collision $t_c$ is

$$ \rho_{\text{wall}}(t_c) \sim \gamma_w^2 m_{\phi}^2 v_{\phi}^2, \quad (6.4) $$

over a thickness of $L_{\text{wall}} = L_0/\gamma_w \sim (\gamma_w m_{\phi})^{-1}$ in the bath frame in the case of a non-supercooled PT. With supercooling, $T_n \ll m_{\phi}, v_{\phi}$, the wall surface tension is larger and the above energy density instead exists over a thickness

$$ L_{\text{wall}} \sim \frac{1}{\gamma_w T_n}, \quad (6.5) $$

in the bath frame, which one finds taking into account regions of energy density both from the leading edge of the bubble in which the potential term dominates and oscillations of the
<figure><img src="image_17.png" /><figcaption>Figure 15. Evolution of the scalar field profile after wall collision with Lorentz factor $\gamma_w = 10$ (left) and $\gamma_w = 100$ (right). We can distinguish the peaked energy distribution stored in the walls even after collision, and the more dilute energy fraction of the scalar field oscillating around the true vacuum. We used the analytical solution in eq. (6.3) of the scalar field equation of motion. This closely resembles the numerical solution shown in figure 16.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16. Time evolution (from top to bottom) of the scalar field profile, its field value on the left and its energy density on the right, after wall collision. We numerically integrated the scalar field equation of motion in eq. (6.2). The wall Lorentz factor, measured as the ratio of the profile periodicity in the wall front to the periodicity in the oscillating condensate is $\gamma_w \simeq 50$.</figcaption></figure>

field about the true vacuum behind the wall [23, 46]. (The temperature does not appear in the above modelling, eqs. (6.1)–(6.3), but we include such a correction to better match the supercooled case.) The energy density of a scalar field oscillating in a quadratic potential is

$$ \rho_{\text{osc}} \sim m_{\phi}^2 v_{\phi}^2 / 2. \qquad (6.6) $$

Conservation of the scalar field energy after collision, i.e. at $t > t_c$, gives

$$ E_{\text{wall}}(t) + E_{\text{osc}}(t) = E_{\text{wall}}(t_c), \qquad (6.7) $$

where $E_{\text{wall}}(t) \simeq \rho_{\text{wall}}L_{\text{wall}}4\pi R^2$ and $E_{\text{osc}}(t) \simeq \rho_{\text{osc}}L_{\text{osc}}4\pi R^2$ are the energies of the bubble wall and oscillating condensate, respectively, and $L_{\text{osc}} = R - R_c$ is the growing distance between the wall and the collision point, with $R_c$ the bubble radius at collision. We obtain

$$ \rho_{\text{wall}}(t) \simeq \rho_{\text{wall}}(t_c) \left(\frac{R_c}{R}\right)^2 - \rho_{\text{osc}} \frac{R - R_c}{L_{\text{wall}}}. \qquad (6.8) $$
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Wall-shell collision: in the blue region, shells interact efficiently with the oscillating scalar condensate left behind collided walls. In the orange region, shells interact efficiently with the peak of scalar field gradient in bubble walls. We show shell particles produced from LO interaction (top), Bremsstrahlung radiation, either transmitted (middle) or reflected (bottom), assuming abelian (left) or non-abelian (right) gauge interaction, see first two rows of table 1. We used $g_b = g_{SM} + g_{emit}$ with $g_{SM} = 106.75$ and $g_{emit} = 10 = N_F$. Grey and black shaded areas as in figure 1.</figcaption></figure>
where we have used that $L_{wall}$ is constant in $t$ after collision. The wall energy is transferred into the oscillating condensate after a time $t_{damp}$ when $\rho_{wall}(t_{damp}) \simeq 0$ in eq. (6.8). We conclude that after collision, the walls continue propagating during a distance

$$L_{osc}(t_{damp}) \equiv R(t_{damp}) - R_c \approx \frac{\rho_{wall}(t_c)}{\rho_{osc}} L_{wall} \sim \gamma_w^2 L_{wall} \sim \frac{\gamma_w}{T_n}, \quad (6.9)$$

before being completely dissolved in the center-of-mass frame of the two colliding walls, which for simplicity we approximate to coincide with the bath frame in this work. In the second equality of eq. (6.9) we have assumed that $R_c/R \simeq 1$, indicating that $L_{osc}(t_{damp}) \ll R_c$. This condition holds in the terminal velocity regime characterized by $\gamma_w/T_n \ll R_c$, and remains acceptable to an order of magnitude in the runaway regime where $\gamma_w/T_n \lesssim R_c$.

## 6.2 Wall-shell collision

Shells of particles produced during the initial propagation of the walls can either interact with the collided wall (the peak of stress-energy tensor) or the oscillating left-over (the diluted part), see figure 14. We model the two scalar field configurations as Bose-Einstein condensate with number densities

$$n_{wall} \sim \gamma_w m_{\phi} v_{\phi}^2, \quad n_{osc} \sim m_{\phi} v_{\phi}^2 \quad (6.10)$$

and energies

$$E_{wall} \sim \gamma_w m_{\phi}, \quad E_{osc} \sim m_{\phi}, \quad (6.11)$$

respectively. The interaction probability of a particle $X$ with the two aforementioned scalar field configurations reads

$$P_{\text{shell-wall}} = n_{\text{wall}} L_{\text{wall}} \sigma_{\text{eff}}(E_{\text{X}}, E_{\text{wall}}) = \frac{m_{\phi} v_{\phi}^2}{T_n} \sigma_{\text{eff}}(E_{\text{X}}, E_{\text{wall}}), \quad (6.12)$$

$$P_{\text{shell-osc}} = n_{\text{osc}} L_{\text{osc}} \sigma_{\text{eff}}(E_{\text{X}}, E_{\text{osc}}) = \frac{\gamma_w m_{\phi} v_{\phi}^2}{T_n} \sigma_{\text{eff}}(E_{\text{X}}, E_{\text{osc}}), \quad (6.13)$$

where the effective cross section, $\sigma_{eff}$, depends on the underlying interaction and shell particle energy $E_X \simeq p_X$.

As an example, consider a gauge interaction between the shell particle and quantum of the wall peak or oscillating condensate. Consider a change in momentum of order $E_X$. This may either be induced through t-channel gauge boson exchange, like Møller scattering, or Feynman diagrams without t-channel gauge bosons, like Compton scattering. The former features a soft enhancement, as we have seen previously in section 4, leading to

$$\sigma_{\text{eff}} \sim \frac{\alpha_D^2}{(-\hat{t}_{\text{min}})} \quad (\text{Møller scattering}) \quad (6.14)$$

where $-\hat{t}_{min} = \text{Min}[E_X^2, E_X E_{wall}]$, giving the required momentum exchange.<fnref n="7" /> In the case of Compton scattering we instead just have

$$\sigma_{\text{eff}} \sim \frac{\alpha_D^2}{E_X E_{\text{wall}}} \quad (\text{Compton scattering}), \quad (6.15)$$

as the cross section is hard scattering dominated.

<footnote n="7">⁷The real scalar quanta of the condensate couple off-diagonally with the imaginary component of the original complex field to the gauge bosons. However, the momentum exchange is large compared to the masses, so effective massless gauge boson exchange is valid to consider.</footnote>
### 6.2.1 Leading order shell

Keeping in mind the above cross sections, we consider first a leading order shell, where $E_X = \Delta m_X^2/T_n$, and assume all the mass is gained in the PT with $\Delta m_X^2 \simeq \alpha_D v_\phi^2$. Then in the case of t-channel gauge boson exchange, such as Møller scattering, we have<fnref n="8" />

$$ \mathcal{P}_{\text{shell-wall}} \approx \text{Max} \left[ \frac{m_\phi T_n}{v_\phi^2}, \frac{\alpha_D}{\gamma_w} \right], \quad \text{and} \quad \mathcal{P}_{\text{shell-osc}} \approx \text{Max} \left[ \frac{\gamma_w m_\phi T_n}{v_\phi^2}, \gamma_w \alpha_D \right]. \quad (6.16) $$

This is applicable, e.g. in the case of a leading order shell of gauge bosons in a non-abelian theory. For interactions without t-channel gauge boson exchange, such as Compton scattering, we instead have

$$ \mathcal{P}_{\text{shell-wall}} \approx \frac{\alpha_D}{\gamma_w}, \quad \text{and} \quad \mathcal{P}_{\text{shell-osc}} \approx \gamma_w \alpha_D. \quad (6.17) $$

This is applicable, e.g. in the case of a leading order shell of gauge bosons in an abelian theory.

### 6.2.2 Gauge boson shell $\alpha_D \lesssim 1$

Here we have $E_X \approx \gamma_w \sqrt{\alpha_D} v_\phi$. Then in the case of an interaction involving t-channel gauge boson exchange, such as Møller scattering, we have

$$ \mathcal{P}_{\text{shell-wall}} \approx \text{Max} \left[ \frac{\alpha_D m_\phi}{\gamma_w^2 T_n}, \frac{\alpha_D^{3/2} v_\phi}{\gamma_w^2 T_n} \right], \quad \text{and} \quad \mathcal{P}_{\text{shell-osc}} \approx \text{Max} \left[ \frac{\alpha_D m_\phi}{\gamma_w T_n}, \frac{\alpha_D^{3/2} v_\phi}{T_n} \right]. \quad (6.18) $$

Again this is applicable for gauge bosons in a non-abelian theory. For Compton scattering, we instead have

$$ \mathcal{P}_{\text{shell-wall}} \approx \frac{\alpha_D^{3/2} v_\phi}{\gamma_w^2 T_n}, \quad \text{and} \quad \mathcal{P}_{\text{shell-osc}} \approx \frac{\alpha_D^{3/2} v_\phi}{T_n}. \quad (6.19) $$

Which is applicable for gauge bosons in an abelian theory. The results here does not depend on whether it is a reflected or transmitted gauge boson shell.

### 6.2.3 Discussion

For our representative examples, the regions where shells interact efficiently with either the gradient energy in bubble walls or the oscillating scalar field formed after the collision are shaded in figure 17. For the abelian case, we have assumed shells of gauge bosons, which do not exchange t-channel gauge bosons in their interactions with the condensate, so the scattering is Compton. For the non-abelian case, we have assumed t-channel gauge boson exchange, so the scattering is Møller.

In the case of leading order shells, we see our estimates show that in much of the parameter space, following the bubble collisions, the shell particles will interact significantly with the oscillating condensate. In the case of transmitted and reflected gauge bosons from

<footnote n="8">⁸Note the first term is non-zero in the limit $\alpha_D \to 0$, but this is not a problem, as we are calculating the scattering probability given the existence of a shell particle (which are not produced for $\alpha_D = 0$), and the cancellation occurs because the required momentum squared exchange $-\hat{t} \approx E_X^2 \propto \alpha_D^2$, cancels the $\alpha_D^2$ factor in the numerator in the $\alpha_D \to 0$ limit.</footnote>
bremsstrahlung, the free streaming regions after bubble collision are much larger. This is due to the higher energies of the shell particles. Note also, that in the case of radiated reflected particle shells, the shells coming from opposite bubbles always interact with each other first, before encountering the other bubble. Therefore, if one is interested in particle production from shell collisions as in bubbletrons [39], one can use the free-streaming predictions even inside the shaded areas for radiated reflected shells, while one cannot inside the shaded areas for LO and radiated transmitted shells, because these shells find the walls before finding each other. Other types of interactions and shells of table 1 can be dealt with in a similar fashion as the examples we have showcased above.<fnref n="9" />

Finally note that our results concerning wall-shells interactions have been derived following a rather simple approximation. For example, our shaded areas would be affected by taking into account that the center-of-mass frame of two bubble walls is not still in the CMB frame, but will realistically move because the colliding walls do not have exactly the same velocity. In addition, more complete treatment would require solving Boltzmann equations in the presence of a background field that varies both in time and in space, see e.g. [84] for a recent attempt in this direction, although in a completely different context and, as far as we understand, within approximations that do not straightforwardly apply to our case. Refining our predictions, while interesting, goes beyond the purposes of this paper.

## 7 Gravitational waves from free-streaming shells

Another application of the present study deals with the GW spectrum resulting from a strong first-order PT. In the free-streaming regime studied in this paper, the plasma dynamics can be reduced to the evolution of ultra-thin relativistic shells. As we argue in this section, this allows the application of established results [85, 86] concerning the GW spectrum from infinitely-thin bubble walls. This extends the applicability of the “bulk flow” model [85, 86], usually restricted to the run-away regime, to the friction-dominated regime provided that shells free stream.

**Energy budget.** When a bubble expand until a radius $R$, the latent heat initially stored in the false vacuum phase $E_{vac} = 4\pi R^3 \Delta V/3$ is converted into kinetic energy of the wall $E_{wall} = 4\pi R^2 \sigma \gamma_w$, where $\sigma$ is the wall surface tension. The fraction $\kappa_{\phi}$ of the latent heat energy stored in the scalar field gradient in bubble walls reads

$$ \kappa_{\phi} = \frac{E_{\text{wall}}}{E_{\text{vac}}} = \frac{3\gamma\sigma R_c}{\Delta V} = \frac{\gamma_{\text{coll}}}{\gamma_{\text{run}}}, \qquad (7.1) $$

where $\gamma_{\text{coll}}$ is the wall Lorentz factor at collision accounting for friction given by eq. (2.14), while $\gamma_{\text{run}}$ is the same quantity neglecting friction, given by eq. (2.11). The remaining fraction of the latent heat

$$ \kappa_{\text{plasma}} = 1 - \kappa_{\phi}, \qquad (7.2) $$

is converted into plasma excitations due to the presence of the friction pressure.

<footnote n="9">⁹For example, we can apply the above method to the Azatov-Vanvlasselaer shells of free streaming non-cold DM produced in a PT studied in ref. [26]. We find shell particle interactions with the oscillating condensate rule out the $T_n < T_{eq}$ case studied therein, while allowed regions of parameter space survive in the $T_n > T_{eq}$ case, where DM with a mass of $10^9$ GeV produced from a weak scale PT could possibly leave an observable signal both in the matter power spectrum and in GW.</footnote>
**Bubble wall contribution.** When $\kappa_{\phi} \simeq 1$ in eq. (7.1), the latent heat is dominantly kept in the scalar field gradients stored in bubble walls. GWs arising from scalar field gradients were initially calculated under the “envelope” approximation, where the assumption was that walls are infinitely thin and the parts that collided were disregarded [87–91]. Subsequently, the collided parts were included into the analysis through the “bulk flow” model both at the analytical [85, 92] and numerical levels [86, 93–95]. The bulk flow model relies on assuming that the collided parts can still be regarded as infinitely thin and expand freely [85]. Using conservation of the total bubble wall energy $E_{\text{wall}}$, this implies that the bubble wall stress-energy momentum tensor in the radial direction decreases as [85]

$$T_{rr} \sim \rho_{\text{wall}} \sim E_{\text{wall}}/(R^2 \times L_{\text{wall}}) \propto 1/R^2, \qquad (7.3)$$

until being dissipated, cf. eq. (6.8). We have assumed the wall thickness $L_{\text{wall}}$ in the bath frame to be constant after collision. The scaling in eq. (7.3) leads to an infrared enhancement of the GW spectrum, characterized by $\Omega_{\text{PT}} \propto f^1$, in contrast to the previously assumed $\Omega_{\text{PT}} \propto f^3$ [85, 86, 92–95]. The GW spectrum from the bulk flow model in the relativistic limit reads [86]

$$\Omega_{\text{PT}}h^2 \simeq \frac{10^{-6}}{(g_*/100)^{1/3}} \left(\frac{H_*}{\beta}\right)^2 \left(\frac{\alpha}{1+\alpha}\right)^2 S_{\text{PT}}(f)S_H(f), \qquad (7.4)$$

with the spectral shape $S_{\text{PT}}(f)$ peaked on $f_{\phi}$

$$S_{\text{PT}}(f) = \frac{3(f/f_{\text{PT}})^{0.9}}{2.1 + 0.9(f/f_{\text{PT}})^3}, \quad f_{\text{PT}} = \left(\frac{a_*}{a_0}\right) 0.8 \left(\frac{\beta}{2\pi}\right), \qquad (7.5)$$

and the redshift factor between percolation “*” and today “0”

$$a_*/a_0 = 1.65 \times 10^{-2} \text{ mHz} \left( \frac{T_{\text{eq}}}{100 \text{ GeV}} \right) \left( \frac{g_{\text{eff,reh}}}{100} \right)^{1/6} H_*^{-1}. \qquad (7.6)$$

We added the correction factor

$$S_H(f) = \frac{(f/f_*)^{2.1}}{1 + (f/f_*)^{2.1}}, \quad f_* = c_* \left(\frac{a_*}{a_0}\right) \left(\frac{H_*}{2\pi}\right), \qquad (7.7)$$

with $c_* = O(1)$ to impose an $f^3$ scaling for emitted frequencies smaller than the Hubble factor $H_*/(2\pi)$ as required by causality [96–99].

**Plasma contribution.** Instead, if $\kappa_{\text{plasma}} \simeq 1$ in eq. (7.2), then the friction pressure in eq. (2.12) prevents bubble walls from accelerating further, see the shaded regions in figure 3, and the latent heat is dominantly dumped into the plasma. Lattice calculations of the GW spectrum from fluid dynamics have been restricted to the regime of non-relativistic velocities $v_w \ll 1$, see e.g. [100, 101]. Apart from the very extreme case where bubbles expand in pure vacuum [95, 102], there is no lattice study of the GW spectrum from first-order PTs in the relativistic limit $v_w \simeq 1$ due to three technical limitations. At first, hydrodynamic equations become highly non-linear once relativistic corrections are added [103]. Second, Lorentz contraction makes the fluid profile thinner than the lattice spacing. Third, the
hydrodynamic regime is expected to break down as soon as the particle free-streaming length is larger than the fluid profile thickness, which applies to (but is not limited to) the case of free-streaming shells.

Instead, in the white region of figure 1 the plasma excitations take the form of relativistic shells which propagate freely, retaining their ultra-thin profile. This suggests that the energy-conserving scaling in eq. (7.3) is applicable not only to bubble walls but also to free-streaming shells, through the substitutions $E_{\text{wall}} \rightarrow E_{\text{shell}}$ and $L_{\text{wall}} \rightarrow L_{\text{shell}}$. From a gravitational perspective, an extremely spatially-walled momentum distribution within the plasma should be indistinguishable from an extremely spatially-walled momentum distribution within the scalar field, provided that the peak width is much smaller than the bubble radius. Therefore, we anticipate that the GW signal for relativistic bubble wall velocity in the combined friction-dominated and free-streaming regime should resembles the GW signal in the run-away regime, described by the bulk flow model in eq. (7.4). This is corroborated by the study [50] conducted in the moderately relativistic regime $\gamma_w \lesssim 10$ where the GW spectrum from free-streaming LO shells indeed resembles the one predicted in the bulk flow model. Hence in all the white regions of this paper, the GW signal is expected to be well approximated by eq. (7.4). Instead, in all the colored (or dashed) regions of this paper, the fate of those relativistic shells and the resulting GW spectrum, but also the friction pressure and bubble wall Lorentz factor, remain unknown. The calculation of GW spectra along with the potential applicability of the bulk flow model $T_{\text{tr}} \propto 1/R^2$, in regions where shells do not free stream, is left for future work.<fnref n="10" />

## 8 Conclusions

First-order cosmological phase transitions (PT) with relativistic bubble walls are predicted in several motivated extensions of the standard model and allow for unique opportunities to solve some of its most pressing mysteries, like the nature of dark matter or the origin of the baryon asymmetry. Importantly, they are the only PT (as opposed to those with slower walls) that result in a gravitational wave (GW) background which one could hope to observe with current and foreseen GW telescopes. A clear understanding of the dynamics associated with relativistic bubble walls has not yet been achieved by the community, possibly jeopardizing their predictions that attracted a lot of attention in recent literature, like those concerning dark matter, baryogenesis, and GW.

One key property of PTs with relativistic walls, that has so far received little attention, is that interactions between particles in the thermal bath and the background field that constitutes the wall necessarily generate secondary energetic particles [23, 26, 44–46, 48]. These form, and accumulate into, ultra-thin energetic shells, possibly very dense, on both sides of the walls. In this paper, we have performed the first systematic study of such shells. In particular:

1. In section 2.3 we have listed all shells ever mentioned in the literature and identified new ones that, to our knowledge, had never been considered before. We have then provided,

<footnote n="10">¹⁰In refs. [93, 94] and [104] it is found that the stress-energy momentum tensor after collision of, respectively, the scalar field gradient and fluid profile in the hydrodynamical limit (i.e. when shells do not free-stream), decreases as $T_{\text{tr}} \propto 1/R^3$. Such a result, if confirmed, suggests additional dissipation with respect to the naive expectation from energy conservation in eq. (7.3), and would imply the need to improve the bulk-flow model.</footnote>
for the first time, a systematic quantitative characterisation of their properties, see table 1 for a summary with pointers to more detailed derivations in the text. This constitutes a necessary ingredient for any subsequent study of shells.

2. We have then determined the regions, of parameter space of first-order PTs, where these shells free stream (meaning they propagate freely through the universe) and therefore where one can trust many of the commonly used predictions from first-order PTs with relativistic walls. A summary of our results is reported in figure 1: we discovered that regions of parameter space widely considered in the literature lies in the region where shells do not free stream and where many properties of the PT (wall velocities, particle production, GWs, etc) will need to be re-evaluated. This is for example the case for the PTs explaining the recently observed GW at pulsar timing arrays, as well as for basically all PTs associated with the spontaneous breaking of a gauge group, unless the gauge coupling is orders of magnitude smaller than one or the breaking scale is much above a PeV.

The processes that can prevent shells free-streaming, that we studied in this paper, are phase-space saturation effects (see section 3 and figure 6), interactions that affect the momentum of shell and/or bath particles, dominantly due to $2 \to 2$ processes (see section 4 and figures 10 and 11), interaction that change both the number and the momenta of shell and/or bath particles, dominantly due to $3 \to 2$ processes (see section 5 and figure 13) and shell interactions with the bubble walls (see section 6 and figure 17).

In exploring the evolution of shells, we obtained results that may be useful also in very different contexts. They include

◇ The computation of the rate of momentum loss of particles through interaction with a surrounding thermal bath in section 4, extending the methods of [26].

◇ An algorithm to simplify the phase-space integration of any 5-point amplitude, see appendix B.

◇ The computation of the $3 \to 2$ squared amplitudes and interaction rates for initial particles with very different momenta, which we report synthetically in section 5 and extensively in appendix C, where we also provide tabulated results for all $3 \to 2$ squared amplitudes of any theory involving initial vectors, fermions or scalars.

One implication of our results, concerning gravitational waves from phase transitions with relativistic walls, is an extension of the applicability of their bulk-flow modeling [85, 86] to some of the cases where GWs are dominantly sourced by the particle shells. GWs from shells are dominant over those from the background scalar field for friction-dominated walls, i.e. when walls do not run away but reach a terminal velocity. The computation of GWs from shells lies outside the regime of validity of hydrodynamics modelings employed so far. We have argued in section 7 that the bulk-flow modeling should provide a good description of GWs sourced by the shells, both in the friction-dominated and runaway regimes, provided that shells free stream.
The study carried out in this paper opens several new avenues of investigation. It identifies the first processes that need to be studied, in case one would want to determine how shells evolve in the regions where they do not free-stream. For example, for shells made of radiated particles (either transmitted or reflected) at a PT associated with the spontaneous breaking of an abelian gauge symmetry, the most important process in inducing a departure from free-streaming is given by $3 \to 2$ interactions, see figure 13. Determining the evolution of shells under these interactions would then allow to understand both the properties of those shells and the wall velocities, and in turn to obtain new predictions that depend on them like particle production via bubbletrons [39], as well as the GW spectrum predicted by the PT [105, 106].

## Acknowledgments

We thank Aleksander Azatov, Jacopo Ghislieri, Ryusuke Jinno, Felix Kahlhoefer, Henda Mansour, Guy Moore, Miguel Vanvlasselaer for very useful discussions. We also thank Miguel Vanvlasselaer for useful comments on the manuscript.

**Funding information.** The work of IB was supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No 101002846, ERC CoG “CosmoChart”. YG is grateful to the Azrieli Foundation for the award of an Azrieli Fellowship. The work of FS is partly supported by the Italian INFN program on Theoretical Astroparticle Physics (TAsP), by the French CNRS grant IEA “DaCo: Dark Connections” and by COST (European Cooperation in Science and Technology) via the COST Action COSMIC WISPers CA21106.

# A Shell production mechanisms

## A.1 Radiated vector bosons

### A.1.1 Number of radiated quanta

Upon acquiring a mass inside bubble walls, a charged particle $a$, either a scalar a fermion or a vector boson, can radiate bremsstrahlung radiation $a \to bc$ [44–46]. We denote by $m_{c,s}$ and $m_{c,h}$ the mass of the radiated vector boson in the symmetric and broken (“Higgs”) phase, by $x = E_c/E_a$ the parent-to-radiated energy ratio, and by $k_\perp$ the transverse momentum of the radiated particle. The parameter space allowed by kinematics can be approximated to [46]

$$ \sqrt{k_{\perp}^2 + m_c(z)^2}/E_a \le x \le 1, \quad 0 \le k_{\perp} \le E_a. \qquad (A.1) $$

The splitting probability reads [46]

$$ dP_{a \to bc} = \zeta_a \frac{dk_{\perp}^2}{k_{\perp}^2} \frac{dx}{x} \Pi(k_{\perp}), \qquad (A.2) $$

where $\Pi(k_\perp)$ contains the IR and UV suppression factors

$$ \Pi(k_{\perp}) \equiv \left( \frac{k_{\perp}^2}{k_{\perp}^2 + m_{c,s}^2} \right)^2 \left( \frac{m_{c,h}^2 - m_{c,s}^2}{k_{\perp}^2 + m_{c,h}^2} \right)^2, \qquad (A.3) $$
<figure><img src="image_20.png" /><figcaption>Figure 18. We display the number of radiated particle ‘c’ (per single incoming particle ‘a’) and their associated energy and emission angle (mean and standard deviation). All quantities are expressed in the bath frame except the energy $E_a$ of incoming particles which is expressed in the wall frame.</figcaption></figure>

and

$$ \zeta_a \equiv \frac{\alpha_D}{\pi} \sum_{b,c} C_{abc}, \qquad \alpha_D \equiv \frac{g_D^2}{4\pi}, \qquad (A.4) $$

where $C_{abc}$ is the charge factor [107]. In this work, we consider the mass in the broken and symmetric phase to be set by the Higgs mechanism and plasma effects respectively

$$ m_{c,h}^2 \approx 2\pi\alpha_D v_{\phi}^2 + \mu^2, \qquad m_{c,s}^2 \approx \mu^2, \qquad (A.5) $$

where $\mu$ is the IR cut-off sets by thermal corrections $\sim g_D T_n$ plus eventual finite-density corrections $\sim g_D \sqrt{n_c/p_c}$ in the non-abelian case, see eq. (2.38) for the precise expression. The number of reflected and transmitted vector bosons per incoming particle reads

$$ N_{a \to bc} = N_R + N_T, \qquad (A.6) $$

where $N_R$ and $N_T$ are the reflected and transmitted population. We compute

$$ N_R \approx \int_0^{E_a} dk_{\perp}^2 \int \frac{\sqrt{k_{\perp}^2 + m_{c,h}^2}}{E_a} dx \frac{dP_{a \to bc}}{dk_{\perp}^2 dx} \approx \zeta_a \log^2 \left( \frac{m_{c,h}}{m_{c,s}} \right), \qquad (A.7) $$

$$ N_T \approx \int_0^{E_a^2} dk_{\perp}^2 \int \frac{1}{\sqrt{k_{\perp}^2 + m_{c,h}^2}} dx \frac{dP_{a \to bc}}{dk_{\perp}^2 dx} \approx 2\zeta_a \log \left( \frac{m_{c,h}}{m_{c,s}} \right) \log \left( \frac{E_a}{m_{c,h}} \right), \qquad (A.8) $$
where we have expanded the results of the integrals for $m_{c,s} \ll m_{c,h} \ll E_a$. Note that eqs. (A.7) and (A.8) are consistent with eq. (59) in [46].

### A.1.2 Emitted energy and splitting angle

The mean energies of reflected and transmitted radiated particles in the bath frame read

$$ \langle E_R \rangle \approx \langle \gamma_w \beta x E_a + \gamma_w \sqrt{(xE_a)^2 - k_{\perp}^2 - m_{c,s}^2} \rangle \approx \frac{3.8\gamma_w m_{c,h}}{\log\left(\frac{m_{c,h}}{m_{c,s}}\right)}, \quad (A.9) $$

$$ \langle E_T \rangle \approx \langle \gamma_w \beta x E_a - \gamma_w \sqrt{(xE_a)^2 - k_{\perp}^2 - m_{c,h}^2} \rangle \approx \frac{0.6\gamma_w m_{c,h}}{\log\left(\frac{E_a}{m_{c,h}}\right)}, \quad (A.10) $$

where the limit $\gamma_w \to +\infty$ was taken. We also compute the standard deviation $\sigma_X \equiv \sqrt{\langle (X - \langle X \rangle)^2 \rangle}$

$$ \sigma_{E_R} \approx \frac{1.6\gamma_w m_{c,h}}{\log^{1/2}\left(\frac{m_{c,h}}{m_{c,s}}\right)}, \quad \sigma_{E_T} \approx \frac{0.4\gamma_w m_{c,h}}{\log^{1/2}\left(\frac{E_a}{m_{c,h}}\right)}. \quad (A.11) $$

The angular distributions read

$$ \langle k_{\perp, \mathcal{R}} \rangle^2 \approx \left( \frac{2m_{c,h}}{\log^2\left(\frac{m_{c,h}}{m_{c,s}}\right)} \right)^2 + T_n^2, \quad \langle k_{\perp, \mathcal{T}} \rangle^2 \approx \left( \frac{0.9m_{c,h}}{\log\left(\frac{m_{c,h}}{m_{c,s}}\right)} \right)^2 + T_n^2. \quad (A.12) $$

and their standard deviation are

$$ \sigma_{k_{\perp, \mathcal{R}}}^2 \approx \left( \frac{0.7m_{c,h}}{\log\left(\frac{m_{c,h}}{m_{c,s}}\right)} \right)^2 + T_n^2, \quad \sigma_{k_{\perp, \mathcal{T}}}^2 \approx \left( \frac{0.7m_{c,h}}{\log^{1/2}\left(\frac{m_{c,h}}{m_{c,s}}\right)} \right)^2 + T_n^2. \quad (A.13) $$

In eqs. (A.12) and (A.13), we have included a second piece of order the nucleation temperature $T_n$ to account for the transverse momentum $p_{a,\perp} \approx T_n$ in the thermal distribution of incoming quanta. All the analytical expressions above are valid in the large log limits and the numerical coefficient are fitted on the numerical calculation. Figure 18 shows a good agreement between the (semi-)analytical and numerical results. We conclude that vector bosons, which are radiated whenever the PT occurs in a gauge symmetry group, are produced with a rate controlled by the gauge coupling constant squared $N \propto g_D^2$, with a large momentum in the bath frame $p_X \approx \gamma_w m_{c,h}$ and highly collimated $k_{\perp} \ll p_X$. We refer to [46] for a calculation of the resulting friction pressure on bubble walls. In the main text, we account for effects from reflected and transmitted particles separately, with associated multiplicity factor $N$ and momentum $p_X$ in the wall frame

$$ N \approx (N_R, N_T), \quad p_X \approx (\langle E_R \rangle, \langle E_T \rangle). \quad (A.14) $$

In all the figures presented in this paper, we incorporate the standard deviations of the kinematic quantities $X = \{E_R, E_T, k_{\perp, \mathcal{R}}, k_{\perp, \mathcal{T}}\}$ by adding them directly to the means

$$ X \to \sqrt{\langle X \rangle^2 + \sigma_X^2}. \quad (A.15) $$

We now discuss a complication which arises when the vector bosons receive a large mass from finite-density effects.
## A.2 Large IR cut-off limit

The expression we have derived in the previous section, section A.1, are only valid in the limit $m_{c,h} \ll m_{c,s}$ (and also $E_a \gg m_{c,h}$). However as we discussed in section 2.4, for non-abelian gauge group finite-density corrections act as an IR cut-off ($X \equiv c$)

$$ \mu^2 \approx 2g_D^2 \frac{n_X}{p_X}, \qquad (A.16) $$

where we neglected the thermal part, see eq. (2.38). In the ultra-relativistic limit, shells are dense and the IR cut-off is larger than the vector boson mass in the broken vacuum $m_{c,h}^\infty \simeq g_D v_\phi / \sqrt{2}$, far away from the wall. We suppose that the vector boson is massless in the symmetric vacuum, so that far away from the wall and up to thermal effects $m_{c,s}^\infty = 0$. Close to the wall within the shells, the vector boson mass receives finite density corrections

$$ m_{c,s}^2 \approx \mu^2, \quad \text{and} \quad m_{c,h}^2 \approx \mu^2 + g_D^2 v_\phi^2/2. \qquad (A.17) $$

The numbers of reflected and transmitted radiated vector bosons — valid both in the limit $m_{c,h} \gg m_{c,s}$ in which eqs. (A.7) and (A.8) were derived and when $\Delta m_c \equiv m_{c,h} - m_{c,s} \lesssim m_{c,h}$ — are well approximated by the following fitting formula

$$ \mathcal{N}_R \approx \int_0^{E_a} dk_\perp^2 \int \frac{\sqrt{\frac{k_\perp^2+m_{c,h}^2}{E_a}}}{\sqrt{\frac{k_\perp^2+m_{c,s}^2}{E_a}}} dx \frac{dP_{a \to bc}}{dk_\perp^2 dx} \approx \zeta_a \left( \frac{\Delta m}{m_{c,s} + \Delta m / L_m^{2/3}} \right)^3, \qquad (A.18) $$

$$ \mathcal{N}_T \approx \int_0^{E_a} dk_\perp^2 \int \frac{1}{\sqrt{\frac{k_\perp^2+m_{c,h}^2}{E_a}}} dx \frac{dP_{a \to bc}}{dk_\perp^2 dx} \approx 2\zeta_a \left( \frac{\Delta m}{m_{c,s} + \Delta m / \sqrt{L_m L_E}} \right)^2, \qquad (A.19) $$

with $L_m = \log(m_{c,h}/m_{c,s})$ and $L_E = \log(E_a/m_{c,h})$. Those formula are valid for all values of $\Delta m_c \equiv m_{c,h} - m_{c,s}$. The radiation pressure becomes suppressed with the mass difference $\Delta m_c$ as

$$ \langle \Delta p_R \rangle \approx \zeta_a \int_0^{E_a} \frac{dk_\perp^2}{k_\perp^2} \int \frac{\sqrt{\frac{k_\perp^2+m_{c,h}^2}{E_a}}}{\sqrt{\frac{k_\perp^2+\mu^2}{E_a}}} \frac{dx}{x} \frac{k_\perp^4}{(k_\perp^2+\mu^2)^2} \left( \frac{m_{c,h}^2 - m_{c,s}^2}{k_\perp^2 + m_{c,h}^2} \right)^2 \times 2xE_a \qquad (A.20) $$

$$ \approx 4 \zeta_a \Delta m_c \left( \frac{\Delta m_c}{m_{c,s} + \Delta m_c / L_m^{1/3}} \right)^3, \qquad (A.21) $$

and

$$ \langle \Delta p_T \rangle \approx \zeta_a \int_0^{E_a} \frac{dk_\perp^2}{k_\perp^2} \int \frac{1}{\sqrt{\frac{k_\perp^2+m_{c,h}^2}{E_a}}} \frac{dx}{x} \frac{k_\perp^4}{(k_\perp^2+\mu^2)^2} \left( \frac{m_{c,h}^2 - m_{c,s}^2}{k_\perp^2 + m_{c,h}^2} \right)^2 \times \frac{k_\perp^2 + m_{c,h}^2}{2xE_a} \qquad (A.22) $$

$$ \approx \zeta_a \Delta m_c \left( \frac{\Delta m_c}{m_{c,s} + \Delta m_c / L_m^{2/3}} \right)^{3/2}, \qquad (A.23) $$

Hence both the number of particles produced and the friction pressure become suppressed in the limit $\Delta m_c \ll m_{c,s}$. We leave the implication of this results for the bubble wall velocity for future research.
## A.3 Radiated scalar bosons

### A.3.1 Number of radiated quanta

We now study the production of scalar bremsstrahlung radiation. We consider the vertex

$$ \phi \to \phi\phi: \qquad V(z) = \lambda\phi(z). \qquad (A.24) $$

In the high-energy limit $p_c^z \gg m_{c,h}$ and thin-wall limit, we can approximate the mode function of c by

$$ \chi_c(z) \approx \exp \left( i \int_0^z p_c^z(z') dz' \right) \approx e^{iE_c z} \exp \left( -\frac{i}{2E_c} \int_0^z (m_c^2(z') + k_\perp^2) dz' \right), \qquad (A.25) $$

and idem for a and b. We deduce

$$ \chi_a(z)\chi_b^*(z)\chi_c^*(z) = e^{i\Delta p_z z}, \qquad (A.26) $$

with

$$ \begin{aligned} \Delta p_z &= \frac{1}{2E_a} \left( -m^2 + \frac{m^2 + k_\perp^2}{1-x} + \frac{m^2 + k_\perp^2}{x} \right) = \frac{1}{2E_a} \left( \frac{k_\perp^2 + (1-x+x^2)m^2}{x(1-x)} \right) \\ &\approx \frac{1}{2E_a} \left( \frac{k_\perp^2 + m^2}{x(1-x)} \right). \end{aligned} \qquad (A.27) $$

The matrix element becomes

$$ \mathcal{M} = \int dz, V(z) \chi_a(z)\chi_b^*(z)\chi_c^*(z) \approx i \left( \frac{V_h}{\Delta p_{z,h}} - \frac{V_s}{\Delta p_{z,s}} \right) \approx i \frac{V_h}{\Delta p_{z,h}}, \qquad (A.28) $$

with $V_s = 0$ and $V_h = \lambda^2 v_\phi^2$. The matrix element squared becomes

$$ |\mathcal{M}|^2 \approx \frac{4E_a^2(1-x)^2 x^2 \lambda^2 v_\phi^2}{(k_\perp^2 + (1-x+x^2)m^2)^2}. \qquad (A.29) $$

The splitting probability reads

$$ dP_{\phi \to \phi\phi} = \frac{dk_{\perp}^2}{64\pi^2 E_a^2} \frac{dx}{x} |M|^2 = \frac{\lambda^2 v_{\phi}^2}{16\pi^2} \frac{dk_{\perp}^2}{(k_{\perp}^2 + (1-x+x^2)m^2)^2} x(1-x)^2 dx. \qquad (A.30) $$

Integrating over the phase space, we obtain the number of radiated scalar bosons (reflected and total)

$$ \mathcal{N}_R = \int_{\mu^2}^{E_a^2} dk_{\perp}^2 \int_{\sqrt{\mu^2 + k_{\perp}^2 / E_a}}^{\sqrt{m^2 + k_{\perp}^2 / E_a}} dx \frac{dP_{a \to bc}}{dk_{\perp}^2 dx} \approx \frac{\lambda^2 v_{\phi}^2}{32\pi^2 E_a^2}. \qquad (A.31) $$

$$ \mathcal{N}_{\text{tot}} = \int_{\mu^2}^{E_a^2} dk_{\perp}^2 \int_{\sqrt{\mu^2 + k_{\perp}^2 / E_a}}^{1} dx \frac{dP_{a \to bc}}{dk_{\perp}^2 dx} \approx \frac{\lambda^2 v_{\phi}^2}{192\pi^2 m_{c,h}^2}, \qquad (A.32) $$
### A.3.2 Emitted energy and splitting angle

The mean energies of reflected and transmitted radiated scalar in the bath frame can be expressed as

$$ \langle E_R \rangle \approx \langle \gamma_w \beta x E_a + \gamma_w \sqrt{(xE_a)^2 - k_{\perp}^2 - m_{c,s}^2} \rangle \approx 3\gamma_w m_{c,h}, \quad (A.33) $$

$$ \langle E_T \rangle \approx \langle \gamma_w \beta x E_a - \gamma_w \sqrt{(xE_a)^2 - k_{\perp}^2 - m_{c,h}^2} \rangle \approx 2\gamma_w \frac{m_{c,h}^2}{E_a} \log \left( \frac{E_a}{m_{c,h}} \right), \quad (A.34) $$

and the angular distributions are given by

$$ \langle k_{\perp, \mathcal{R}} \rangle \approx \langle k_{\perp, \mathcal{T}} \rangle \approx m_{c,h}. \quad (A.35) $$

We conclude that in the case of scalar Bremsstrahlung, which would be the leading process if the phase transition occurs in a global (not gauged) symmetry group, scalar bosons are radiated with a large momentum in the bath frame only if they are reflected. However, their abundance relative to the transmitted ones is suppressed $N_R/N_{tot} \approx 6(m_{c,h}/E_a)^2 \ll 1$.

### B 3 → 2 final phase space integration

In this appendix we describe an algorithm that allows one to perform the final phase space integral of 3 → 2 processes for any scattering amplitudes, where we define the momenta of incoming and outgoing particles as in the main text, i.e. $p_1 p_2 p_5 \to p_3 p_4$. We anticipate that our result, although the integration will be performed in a specific frame, is Lorentz invariant, and can therefore be expressed in terms of three independent scalar products $p_1 \cdot p_2, p_1 \cdot p_5$, and $p_2 \cdot p_5$.

The most complicated expressions of scalar products, in the denominators of squared amplitudes $|M|^2$, are those that come from interference terms. One 5-point amplitude, consisting of only 3-point or 4-point interactions, contains at most two denominators from propagators, which are different from each other. In an interference term this would be multiplied by two more denominators, all possibly different. Therefore we should be able to compute, and integrate over final phase space, terms with at least 4 different denominators, e.g.

$$ |M|^2 \supset \frac{1}{s_{13}^a s_{24}^b} \cdot \frac{1}{s_{23}^c s_{54}^d}, \quad (B.1) $$

where we have defined $s_{ij} = (p_i \pm p_j)^2$ with − for one incoming and one outgoing momenta, i.e. for $i, j = 3, 5$, and + for two incoming or outgoing momenta, i.e. for $i, j = 1, 5$ or $i, j = 3, 4$ (this corresponds to $s_{ij} = (p_i + p_j)^2$ if momenta were all defined as incoming). One possible approach to recast those integrals is the method of Feynman parameters. In loop-computations one usually completes the square to get rid of the linear terms, and then performs the loop integration. Here the approach would be to perform a SO(3) rotation to get rid of any azimuthal angle in the denominator, such that the integral over the polar angle can be performed. After doing so one ends up with an integral over 3 Feynman parameters. However, these integrals are not textbook integrals and evaluating them efficiently is not guaranteed.

In the following we will therefore describe an algorithm which is able to perform the final phase space integration for all possible 3 → 2 amplitudes.
## B.1 Partial fraction decomposition of scalar products

We start from the following expression for the squared $3 \to 2$ amplitude,

$$ |\mathcal{M}|^2 = \sum_n C_n \frac{\prod_{ij} s_{ij}^k}{s_{13}^a s_{14}^b s_{23}^c s_{24}^d s_{53}^e s_{54}^f} . \qquad (B.2) $$

In addition, $|\mathcal{M}|^2$ can also depend on combinations of $s_{12}$, $s_{15}$, $s_{25}$, $s_{34}$, but these do not depend on the final phase space integration. Therefore, we absorb them into the prefactors $C_n$.

The next step is to reduce the numbers of different factors in the denominator. We achieve this by the method of partial fraction decomposition and by exploiting energy momentum conservation. This can be achieved systematically, in an easy way, by multiplying with one non-trivial representation of the identity. We derive several non trivial representations of the identity by writing energy momentum conservation in the following way,

$$ p_k \cdot (p_1 + p_2 + p_5 - p_3 - p_4) = 0, \qquad (B.3) $$

which represents five independent equations for $k = 1, 2, 3, 4, 5$. We find, for example for $k = 1$,

$$ p_1 \cdot p_1 + p_1 \cdot p_2 + p_1 \cdot p_5 = p_1 \cdot p_3 + p_1 \cdot p_4, \qquad (B.4) $$

or equivalently

$$ s_{12} + s_{15} - 2m_1^2 - m_2^2 - m_5^2 - m_3^2 - m_4^2 = -s_{13} - s_{14}. \qquad (B.5) $$

Performing this manipulation for all $k$ we find several equations for the identity,

$$ 1 = \frac{s_{13} + s_{14}}{C_1}, \qquad (B.6) $$

$$ 1 = \frac{s_{23} + s_{24}}{C_2}, \qquad (B.7) $$

$$ 1 = \frac{s_{53} + s_{54}}{C_5}, \qquad (B.8) $$

$$ 1 = \frac{s_{13} + s_{23} + s_{53}}{C_3}, \qquad (B.9) $$

$$ 1 = \frac{s_{14} + s_{24} + s_{54}}{C_4}, \qquad (B.10) $$

where the constants $C_k$ do not depend on the final phase space integration.

More relevant identities, which can be derived in similar fashion, are

$$ 1 = \frac{s_{14} + s_{23} + s_{53}}{C_{14}}, \qquad (B.11) $$

$$ 1 = \frac{s_{13} + s_{24} + s_{54}}{C_{13}}, \qquad (B.12) $$

$$ 1 = \frac{s_{13} + s_{24} + s_{53}}{C_{24}}, \qquad (B.13) $$

$$ 1 = \frac{s_{14} + s_{23} + s_{54}}{C_{23}}, \qquad (B.14) $$

$$ 1 = \frac{s_{13} + s_{23} + s_{54}}{C_{54}}, \qquad (B.15) $$

$$ 1 = \frac{s_{14} + s_{24} + s_{53}}{C_{53}}. \qquad (B.16) $$
All C's are guaranteed to not introduce any additional singularities into the amplitude.

It is straightforward to see that by strategic multiplication with these identities one can step by step reduce the powers of particular denominators up to eliminate some (but not all) of them. After the full reduction has been completed, we end up with the following basic integrands in $|M|^2$,

$$ \mathcal{I}_1 \supset \frac{1}{s_{13}^a}, \qquad \mathcal{I}_2 \supset \frac{1}{s_{13}^a s_{24}^b}, \qquad \mathcal{I}_3 \supset \frac{s_{13}^m}{s_{53}^a}, \qquad (B.17) $$

where any 3 can also be a 4 and vice versa.

Note that any term in the numerator which depends on the final phase space integration was already present in the amplitude, before our reduction procedure. In other words, our reduction algorithm does not generate any additional term in the numerator.

## B.2 Parametrization of the integration region

It is most convenient to evaluate the two-particle final phase space integral in the center-of-mass frame, where after eliminating the delta-distribution ensuring energy momentum conservation we are left with an integral over two angles,

$$ \begin{aligned} R_2((p_3 + p_4)^2) &= \frac{d^3 p_3}{(2\pi)^3 2E_3} \frac{d^3 p_4}{(2\pi)^3 2E_4} \delta^{(4)}(p_1 + p_2 + p_5 - p_3 - p_4) \\ &= \frac{\sqrt{\lambda((p_3 + p_4)^2, m_3^2, m_4^2)}}{32\pi^2 (p_3 + p_4)^2} \int d\Omega, \end{aligned} \qquad (B.18) $$

where $\lambda(x, y, z) = x^2 + y^2 + z^2 - 2xy - 2yz - 2zx$ is the Källén function and $\int d\Omega$ is the standard integral over the $S^2$-sphere,

$$ \int d\Omega = \int_0^\pi d\theta_{13} \sin(\theta_{13}) \int_0^{2\pi} d\varphi_{13} = \int_{-1}^{+1} dz_{13} \int_0^{2\pi} d\varphi_{13} \qquad (B.19) $$

It is worthy to point out that the final phase space integral is invariant under the exchange of the two momenta of the final states, $p_3 \leftrightarrow p_4$. This is mostly useful for simplifying the amplitude in order to safe computational resources by not calculating the same integral twice, since $p_3 \leftrightarrow p_4$ implies $z_{13} \leftrightarrow -z_{13}$. The integral is also invariant under $O(3)$ rotations on the $S^2$-sphere, which we will exploit heavily in order to being able to compute some of the integrals.

We now state our choice for the parameterization of the five momenta,

$$ p_{1,\text{com}} = \begin{pmatrix} \sqrt{m_1^2 + p_I^2} \\ 0 \\ 0 \\ p_I \end{pmatrix}, \quad p_{2,\text{com}} = \begin{pmatrix} \sqrt{m_2^2 + p_J^2} \\ p_J \sin(\theta_{12}) \\ 0 \\ p_J \cos(\theta_{12}) \end{pmatrix}, \quad p_{5,\text{com}} = \begin{pmatrix} \sqrt{m_5^2 + p_K^2} \\ -p_J \sin(\theta_{12}) \\ 0 \\ -p_I - p_J \cos(\theta_{12}) \end{pmatrix}, \quad (B.20) $$

$$ p_{3,\text{com}} = \begin{pmatrix} \sqrt{m_3^2 + p_F^2} \\ p_F \sin(\theta_{13}) \cos(\varphi_{13}) \\ p_F \sin(\theta_{13}) \sin(\varphi_{13}) \\ p_F \end{pmatrix}, \quad p_{4,\text{com}} = \begin{pmatrix} \sqrt{m_4^2 + p_F^2} \\ -p_F \sin(\theta_{13}) \cos(\varphi_{13}) \\ -p_F \sin(\theta_{13}) \sin(\varphi_{13}) \\ -p_F \end{pmatrix}. \quad (B.21) $$
The on-shell condition $p_{5,com}^2 = m_5^2$ fixes $p_K = \sqrt{p_I^2 + 2p_I p_J \cos(\theta_{12}) + p_J^2}$. This allows us to replace the angle $\theta_{12}$ by the 'length' $p_K$, which can be interpreted as describing a triangle either by two lengths $p_I, p_J$ and an angle $\theta_{12}$, or by three lengths $p_I, p_J, p_K$. This constrains these parameters to the region $p_I \ge 0, p_J \ge 0, p_K \ge 0, p_K \le p_I + p_J, p_K \ge |p_I - p_J|$. Conservation of energy also fixes $p_F$ via the equation

$$ \sqrt{m_1^2 + p_I^2} + \sqrt{m_2^2 + p_J^2} + \sqrt{m_5^2 + p_K^2} = \sqrt{m_3^2 + p_F^2} + \sqrt{m_4^2 + p_F^2}. \qquad (B.22) $$

After the integration the result depends only three variables $p_I, p_J, p_K$. However, we can exploit the Lorentz invariance of the amplitude, and re-express the result in terms of scalar products. This is because $p_1 \cdot p_2 = p_{1,com} \cdot p_{2,com}$, etc. defines a solvable system of three coupled algebraic equations which lets us perform the substitution

$$ \{p_I, p_J, p_K\} \Rightarrow \{p_{1,com} \cdot p_{2,com}, p_{1,com} \cdot p_{5,com}, p_{2,com} \cdot p_{5,com}\} \Rightarrow \{p_1 \cdot p_2, p_1 \cdot p_5, p_2 \cdot p_5\}. \qquad (B.23) $$

This trick avoids having to work out the explicit Lorentz transformation from the center-of-mass frame to your favorite frame, since one just computes the three scalar products in the preferred reference frame and substitutes them into the result.

Note that because of $(p_3 + p_4)^2 = (p_1 + p_2 + p_5)^2$ the scalar product $p_3 \cdot p_4$ is not independent. This is consistent with the fact that integration over two integration variables should reduce the number of independent scalar products from five to three.

## B.3 Basis integrals

Let us now proceed with the computation of the basis integrals of equation (B.17). For a more transparent notation we are going to replace complicated expressions which do not depend on the final phase space angles with simple variables $c_l$. Also the substitution of integration variable $\cos(\theta_{13}) = z_{13}$ is going to be performed without extra comments.

The analytic integrals have been verified by comparing the results to those of numeric integration.

### B.3.1 $\mathcal{I}_1$

Consider the integral

$$ \int d\Omega \mathcal{I}_1 \supset \int d\Omega \frac{1}{s_{13}^a} \qquad (B.24) $$

$$ = \int_{-1}^{+1} dz_{13} \int_{0}^{2\pi} d\varphi_{13} \frac{1}{(c_0 + c_1 \cdot z_{13})^a}. \qquad (B.25) $$

The azimuthal integration evaluates trivially to $2\pi$, and the integration over $z_{13}$ is a textbook integral which evaluates either to a rational function or a logarithm, depending on the value of $a$.

Now consider an integral of similar structure, but different momentum,

$$ \int d\Omega \mathcal{I}_1 \supset \int d\Omega \frac{1}{s_{23}^a} \qquad (B.26) $$

$$ = \int_{0}^{\pi} d\theta_{13} \sin(\theta_{13}) \int_{0}^{2\pi} d\varphi_{13} \frac{1}{[c_0 + c_2 \sin(\theta_{13}) \cos(\varphi_{13}) + c_1 \cos(\theta_{13})]^a}. \qquad (B.27) $$
This integral can no longer be evaluated quickly. The solution is to exploit the invariance of the final phase space integral under $O(3)$-rotations, and rotate the vectors $p_{3,com}$, $p_{4,com}$ by an angle such that they are aligned with the direction of $p_{2,com}$. Note that because of the linearity of the integral one has to apply this rotation only to this term, not to all the others. However, scalar products in the same term are invariant under rotations, since they are part of the Lorentz group, and therefore we do not need to perform any transformations an scalar products independent of $p_{3,com}$ and $p_{4,com}$.

Therefore, the method of integrating structures of type $\mathcal{I}_1$ is performing an appropriate $O(3)$ rotation, and then evaluating the integral over a rational function.

### B.3.2 $\mathcal{I}_2$

Consider the integral

$$ \int d\Omega \mathcal{I}_2 \supset \int d\Omega \frac{1}{s_{23}^a s_{54}^b}. \qquad (B.28) $$

One possible way to proceed would be to perform a partial fraction decomposition with respect to the variable $z_{13}$, and then apply different $O(3)$ rotations to the different terms.

However, an easier way is Feynman parameterization, which is also easier to implement as a generally applicable algorithm. For two denominators we have

$$ \frac{1}{A^a B^b} = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \int_0^1 dx \frac{x^{a-1} \bar{x}^{b-1}}{(x A + \bar{x} B)^{a+b}}, \quad \bar{x} = 1-x. \qquad (B.29) $$

Then for the above example we have

$$ xA + \bar{x}B = c_0' - 2(xE_2E_3 + \bar{x}E_5E_4) + 2(x \vec{p}_2 \cdot \vec{p}_3 + \bar{x} \vec{p}_5 \cdot \vec{p}_4) \qquad (B.30) $$

$$ = c_0(x) + 2(x\vec{p}_2 - \bar{x}\vec{p}_5) \cdot \vec{p}_3 \qquad (B.31) $$

where $E_k = p_{k,com}^0$ are the energies in the center-of-mass frame, and for readability we omitted the index com indicating the evaluation in the center-of-mass frame. In the last step we used that $p_{3,com} = -p_{4,com}$.

In order to evaluate the integral over the azimuthal and polar angle we make again use of the invariance under $O(3)$-rotations and rotate in such a way that the vector $\vec{p}_3$ is parallel to the vector $x\vec{p}_2 + \bar{x}\vec{p}_5$. We can then write the scalar product as

$$ 2(x\vec{p}_2 - \bar{x}\vec{p}_5) \cdot \vec{p}_3 = 2|x\vec{p}_2 - \bar{x}\vec{p}_5||\vec{p}_3| \cos(\theta_{13}), \qquad (B.32) $$

such that the integral over $\mathcal{I}_2$ reduces to

$$ \int d\Omega \mathcal{I}_2 \supset \int d\Omega \frac{1}{s_{23}^a s_{54}^b} \qquad (B.33) $$

$$ = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \int_0^1 dx \int_{-1}^{+1} dz_{13} \int_0^{2\pi} d\varphi_{13} \frac{x^{a-1} \bar{x}^{b-1}}{[c_0(x) + c_1(x)z_{13}]^{a+b}}. \qquad (B.34) $$

The integration over the azimuthal angle $\varphi_{13}$ evaluates trivially to $2\pi$. Since this expression only appears for at least denominators, i.e. $a \ge 1, b \ge 1$, we have $a + b \ge 2$. Therefore,