Standard Model Baryon Number Violation at Zero Temperature from Higgs Bubble Collisions

TL;DR

This paper calculates baryon number violation from Higgs bubble collisions at zero temperature, showing it can match thermal sphaleron rates, with implications for baryogenesis.

Contribution

First computation of baryon number violation at zero temperature from Higgs bubble collisions using large-scale lattice simulations.

Findings

Baryon violation rate comparable to thermal sphalerons at electroweak temperatures.

Dependence of transition rates on scalar potential shape and bubble wall Lorentz factor.

Estimated baryon asymmetry assuming CP violation during phase transition.

Abstract

We compute for the first time baryon number violation at zero temperature from Higgs bubble collisions and find that it can be of the same order as that from thermal sphalerons in the symmetric phase at electroweak temperatures. We study the dependence of the rate of Chern--Simons number transitions on the shape of the scalar potential and on the Lorentz factor of the bubble walls at collision via large-scale (3+1)D lattice simulations of the Higgs doublet and SU(2) gauge fields. We estimate the resulting baryon asymmetry assuming some CP-violating source activated by the Higgs-field variation during the phase transition.

Equations

$${\Gamma }_{\mathrm{CS}}=\frac{1}{L^3}\frac{d\Delta N_{\mathrm{CS}}^2(t)}{dt}\equiv \frac{d\Delta n_{\mathrm{CS}}^2(t)}{dt},$$

$$\Delta N_{\mathrm{CS}}^2(t)\equiv \text{expectationvalue}{N}_{\mathrm{CS}}^2(t)-\text{expectationvalue}{N_{\mathrm{CS}}(t)}^2$$

$$\Delta n_{\mathrm{CS}}^2=c\cdot {\gamma }_{\star }^{a}f(t/R_{\star }),$$

$$\mathcal{L}=-|D_{\mu }\varphi {|}^2-\frac{1}{2}\mathrm{Tr}[W_{\mu \nu }{W}^{\mu \nu }]-V(\varphi ),V(\varphi )=\frac{m_h^2}{2v^2}{\left(|\varphi {|}^2-\frac{v^2}{2}\right)}^2+\frac{1}{{\Lambda }^2}{\left(|\varphi {|}^2-\frac{v^2}{2}\right)}^3+\dots .$$

$$\frac{v^4}{m_h^2}<{\Lambda }^2<\frac{3}{2}\frac{v^4}{m_h^2},$$

$$\Delta V=\frac{m_h^2{v}^2}{8}\left(1-1/\ell \right),ϵ=\frac{1}{4}\left(1+3/\ell \right){\left(2-3/{\ell }^2\right)}^2,\ell \equiv \frac{{\Lambda }^2{m}_h^2}{v^4}.$$

$$h^{\prime \prime }(r)+\frac{3}{r}{h}^{\prime }(r)={\partial }_{h}V(h/\sqrt{2}),$$

$${\varphi }_b(r)=\frac{1}{\sqrt{2}}{e}^{-i{\sigma }_i{\theta }_i}(0,h_c(r)),$$

$$D_{\mu }{D}^{\mu }\varphi -{\partial }_{{\varphi }^{\ast }}V(\varphi )=0,$$

$$D_{\nu }{W}^{\mu \nu }=J_a^{\mu }{T}_a,J_a^{\mu }=2g\mathrm{Im}[{\varphi }^{†}{T}_a{D}^{\mu }\varphi ].$$

$${\rho }_{\mathrm{tot}}\equiv |\dot{\varphi }{|}^2+|D_i\varphi {|}^2+V(\varphi )+\frac{1}{2}\sum _{a,i}{W_{0i}^a}^2+\frac{1}{2}\sum _{a,i<j}{W_{ij}^a}^2\equiv {\rho }_K+{\rho }_G+{\rho }_V+{\rho }_E+{\rho }_B$$

$$U_{\mu \nu }\equiv U_{\mu }{U}_{\nu ,+\mu }{U}_{\mu ,+\nu }^{†}{U}_{\nu }^{†}=e^{-igd{x}^{\mu }d{x}^{\nu }{W}_{\mu \nu }(1+O(d{x}^{\mu }))}$$

$$W_{\mu \nu }=\frac{i}{2gd{x}^{\mu }d{x}^{\nu }}(U_{\mu \nu }-U_{\mu \nu }^{†})+O(d{x}^2).$$

$$D_{\mu }^{\pm }\varphi =\pm \frac{1}{d{x}^{\mu }}\left(U_{\pm \mu }{\varphi }_{\pm \mu }-\varphi \right)+O(dx),$$

$$\varphi (t)$$

$${\pi }_{\varphi }(t+dt/2)$$

$$\ddot{\varphi }(t)=\frac{1}{d{x}^2}\sum _i(U_i{\varphi }_{+i}-2\varphi +U_{-i}{\varphi }_{-i})-{\partial }_{{\varphi }^{\ast }}V(\varphi ),$$

$$U_i(t)=$$

$${\pi }_{W,i}^a(t+dt/2)={\pi }_{W,i}^a(t-dt/2)+dt{\dot{\pi }}_{W,i}^a(t)+O(d{t}^2)$$

$${\dot{\pi }}_{W,i}^a={\dot{W}}_{0i}^a=(D_j{W}_{ji}{)}^a+2g\mathrm{Im}[{\varphi }^{†}{T}_a{D}_i\varphi ].$$

$${\dot{\pi }}_{W,i}^a=\frac{2i}{gd{x}^3}\mathrm{Tr}\left[\left(U_{ji}-U_{j,-j}^{†}{U}_{ji,-j}{U}_{j,-j}\right)T_a\right]+\frac{2g}{dx}\mathrm{Im}[{\varphi }^{†}{T}_a{U}_i{\varphi }_{+i}],$$

$$(D_i{W}_{0i}{)}^a=2g\mathrm{Im}[{\varphi }^{†}{T}_a\dot{\varphi }],$$

$$\sum _i\frac{{\pi }_{W,i}(t,x)-U_{i,-i}^{†}(t){\pi }_{W,i}(t,x-\hat{i})U_{i,-i}(t)}{dx}=2g\mathrm{Im}[{\varphi }^{†}(t,x)T_a{\pi }_{\varphi }(t,x)]T_a,$$

$${\Delta }_g\equiv \frac{⟨|\mathrm{RHS}-\mathrm{LHS}|⟩}{⟨|\mathrm{RHS}+\mathrm{LHS}|⟩}.$$

$$⟨{\rho }_E⟩=\sum _{i,a}\frac{1}{2}⟨({\pi }_{W,i}^a{)}^2⟩,⟨{\rho }_B⟩=\sum _{i,a}\frac{1}{2}⟨(B_i^a{)}^2⟩=\sum _{i,j,a}\frac{1}{4}⟨(W_{jk}^a{)}^2⟩,$$

$$N_W=\frac{1}{24{\pi }^2}\int d^{3}x{ϵ}_{ijk}{\Phi }_{ijk}=\frac{1}{8{\pi }^2}\int d^{3}x\left({\Phi }_{123}-{\Phi }_{132}\right),{\Phi }_{ijk}\equiv \mathrm{Tr}[{\Phi }^{†}({\partial }_i\Phi ){\Phi }^{†}({\partial }_j\Phi ){\Phi }^{†}({\partial }_k\Phi )],$$

$$j_{\mathrm{CS}}^{\mu }\equiv \frac{g^2}{16{\pi }^2}{ϵ}^{\mu \nu \rho \sigma }\mathrm{Tr}\left[W_{\nu }{W}_{\rho \sigma }+\frac{2}{3}ig{W}_{\nu }{W}_{\rho }{W}_{\sigma }\right].$$

$$n_{\mathrm{CS}}\equiv j_{\mathrm{CS}}^0=\frac{g^2}{16{\pi }^2}{ϵ}^{ijk}\mathrm{Tr}\left[W_i{W}_{jk}+\frac{2}{3}ig{W}_i{W}_j{W}_k\right].$$

$$N_{\mathrm{CS}}(t)=\int d^{3}x{n}_{\mathrm{CS}}(t,x).$$

$$n_{\mathrm{CS}}=\frac{g^2}{64{\pi }^{2}d{x}^3}{ϵ}^{ijk}\text{Tr}\left[(U_i-U_i^{†})\left((U_{jk}-U_{jk}^{†})-\frac{1}{3}(U_j-U_j^{†})(U_k-U_k^{†})\right)\right]+O(d{x}^3).$$

$${\Gamma }_{\mathrm{CS}}\equiv \underset{L,t\to \infty }{\lim }\frac{⟨N_{\mathrm{CS}}^2(t)⟩-⟨N_{\mathrm{CS}}(t){⟩}^2}{L^{3}t},$$

Full text

††preprint: DESY-25-120

Introduction– Baryon number ($B$) is a global symmetry of the Standard Model (SM) whose violation has never been observed experimentally. It is however efficiently violated by the chiral anomaly at high temperature $T$ through $SU(2)$ vacuum transitions producing integer Chern–Simons (CS) number variations $\Delta N_{\rm CS}\propto\Delta B$, the sphaleron processes, which play a key role in essentially all models of baryogenesis, whether at the electroweak scale Kuzmin et al. (1985) or well beyond, as in leptogenesis Fukugita and Yanagida (1986); Davidson et al. (2008). These transitions can also occur at zero temperature when the Higgs field is driven far from equilibrium, e.g. if the SM electroweak phase transition (EWPT) is quenched, i.e. happens much faster than when thermally induced. This has been exploited in models of (local) cold baryogenesis where inflation takes place at the EW scale and the inflaton coupling to the Higgs quenches the Higgs mass parameter, leading to a tachyonic instability that triggers $\Delta N_{\rm CS}\neq 0$ at every point in space Garcia-Bellido et al. (1999); Krauss and Trodden (1999); Cornwall and Kusenko (2000); Copeland et al. (2001); Garcia-Bellido et al. (2004); Enqvist et al. (2010). The rate of $B$ violation from such a quenched EW crossover (i.e. not first-order EWPT) has been studied on the lattice in Smit and Tranberg (2002); Tranberg and Smit (2003); van der Meulen et al. (2006); Tranberg and Smit (2006); Tranberg et al. (2007); Tranberg (2011); Tranberg and Wu (2012); Mou et al. (2015, 2017a, 2017b, 2018b, 2018a).

In this letter, we compute for the first time the rate of $B$ violation at $T=0$ arising from Higgs bubble collisions. We show that this new source of $B$ violation can be of the same order as the one from thermal sphalerons. This opens up the possibility for a new mechanism of EW baryogenesis in models where the reheat temperature associated with a first-order EWPT never exceeds the 130 GeV sphaleron freeze-out temperature D’Onofrio et al. (2014). In this context, EW baryogenesis could be potentially realised even for very strongly-first order supercooled EWPT i.e. for bubble wall velocity exceeding the speed of sound, a regime where the standard non-local EW baryogenesis charge transport mechanism is highly suppressed Dorsch et al. (2021). Such frameworks are particularly motivated in Composite Higgs models where a light dilaton drives the EWPT Konstandin and Servant (2011a, b); Servant (2014); von Harling and Servant (2017); Bruggisser et al. (2018, 2023b, 2023a). See also Refs. Baldes et al. (2021); Azatov et al. (2021); Dichtl et al. (2024); Cataldi and Shakya (2024) for baryogenesis with ultra-relativistic walls but additional sources of $B$ violation. More generally, this new source of $B$ violation may impact common models of EW baryogenesis as a new source of $B$ washout.

The production of CS number in this work can be attributed to the dynamics of the electroweak textures. These are configurations where the scalar field is in the vacuum manifold $M=SU(2)\times U(1)/U(1)$ everywhere in space, but defines a nontrivial map from $M$ onto the three-sphere $S^{3}=\mathbb{R}^{3}\cup\{\infty\}$, and are classified by third homotopy group $\pi_{3}(M)=\mathbb{Z}$. Thus, the energy of a texture is entirely due to the field gradient. In this paper, we consider $SU(2)$ as a simplified version of the SM gauge group, i.e. we set the hypercharge gauge coupling $g^{\prime}=0$, as this captures the main physics of interest. Textures can arise from the different $SU(2)$-phase orientations of the Higgs doublet $\phi$ across space after spontaneous symmetry breaking. A texture is labeled by the (integer) Higgs winding number $N_{W}=\frac{1}{24\pi^{2}}\int d^{3}x\,\epsilon_{ijk}\text{Tr}(\partial_{i}\Phi)\Phi^{\dagger}(\partial_{j}\Phi)\Phi^{\dagger}(\partial_{k}\Phi)\Phi^{\dagger}$, where $\Phi=(i\sigma_{2}\phi^{*},\phi)/|\phi|^{2}\in SU(2)$, and the gauge field Chern–Simons number $N_{\rm CS}(t)-N_{\rm CS}(0)=\frac{g^{2}}{16\pi^{2}}\int_{0}^{t}dt^{\prime}\int d^{3}x\,\text{Tr}[W^{\mu\nu}\tilde{W}_{\mu\nu}]$, where the dual of the $SU(2)$ field strength $W^{\mu\nu}$ is $\tilde{W}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}W^{\rho\sigma}$. A prototypical texture with $N_{W}=1$ has $W_{\mu\nu}=0$ and $\phi(x)=e^{-i\eta(r)x_{i}\sigma_{i}/r}(0,v)$, where $\eta(r)=\pi\tanh(r/r_{t})$, $\sigma_{i}$ are Pauli matrices, $v=246\,$GeV is the Higgs expectation value and $r_{t}$ is the radius. The minimal-energy (vacuum) configuration has $N_{\rm CS}=N_{W}$, while a texture has $N_{\rm CS}-N_{W}\neq 0$. Textures with characteristic length scale smaller than the inverse $W$ mass $m_{W}^{-1}$ decay by changing $N_{W}$, while those larger than $m_{W}^{-1}$ get ‘dressed’ by the gauge fields and decay by altering $N_{\rm CS}$ Turok and Zadrozny (1991). Thus, crucially, the formation and subsequent decay of $SU(2)$-textures can lead to a change in the CS number $\Delta N_{\rm CS}$, which eventually leads to baryon number production via the quantum anomaly, $\Delta B=N_{F}\Delta N_{\rm CS}$, where $N_{F}=3$ is the number of flavours.

Higgs bubble collisions generate Higgs inhomogeneities that produce $N_{W}$ and $N_{\rm CS}$, possibly directly or through formation and decay of textures. Our simulations indeed observe such production, see Fig. 1. To quantify the change in the CS number, we define the CS number diffusion rate

[TABLE]

with Chern-Simons variance

[TABLE]

where $\langle.\rangle$ is the statistical average over realizations and $L^{3}$ is the volume. Notice that $\langle N_{\rm CS}\rangle=0$, as we do not include sources of CP violation. It will be useful to compare Eq. (1) with the SM thermal sphaleron rate in the symmetric phase Moore (1998, 2000); D’Onofrio et al. (2014), $\Gamma_{\rm th}\sim 25\alpha_{w}^{5}T^{4}$, and with the rate from a quenched EW crossover Garcia-Bellido et al. (1999); Tranberg and Smit (2003); van der Meulen et al. (2006), $\Gamma_{\rm tach}\equiv\alpha_{w}^{4}T_{\rm eff}^{4}$, where $T_{\rm eff}$ parametrizes the efficiency of baryon number violation induced by the tachyonic instability and depends on the amount of quenching of the EW crossover.

Real-time simulations– We rely on SM-physics only, except that we assume a first-order EWPT triggered by a new light scalar, with the Higgs field tracking its bubble profile. For simplicity, however, we model the transition using a modified single-field Higgs potential in which $\phi=0$ is a local minimum, see Suppl. Mat. D. We solve the Higgs and gauge boson equations of motion numerically on a periodic cubic grid with $N_{x}^{3}\in[500^{3}-1500^{3}]$ points, setting $m^{2}_{W}/m^{2}_{h}=0.4$ as in the SM (neglecting hypercharge). We expect our main results on the CS rate to be barely affected when the full $SU(2)\times U(1)$ is considered. The main model-dependent parameter is the degeneracy parameter $\epsilon$, which characterizes the scalar potential shape, and plays an important role in the field dynamics during and after bubble collisions Konstandin and Servant (2011b); Jinno et al. (2019); Gould et al. (2021). It is defined by $\epsilon=\frac{V_{\rm barr}}{V_{\rm barr}+|\Delta V|}\in(0,1)$, where $V_{\rm barr}$ is the barrier height and $\Delta V$ is the energy difference between the $\phi=0$ and $\phi=v/\sqrt{2}$ vacua Jinno et al. (2019). We have $\epsilon\to 1$ for degenerate vacua, $|\Delta V|\ll V_{\rm barr}$, and $\epsilon\to 0$ for non-degenerate vacua, $|\Delta V|\gg V_{\rm barr}$. For large $\epsilon$, the walls tend to reflect upon collision, temporarily recreating a region of false vacuum that eventually shrinks after a several iterations. Instead, for $\epsilon\ll 1$ the barrier is not effective in trapping the field in the false vacuum, and the true minimum is eventually reached after just a few oscillations.

We nucleate $N_{b}=10$ bubbles simultaneously, but verify that our results are independent of $N_{b}$, see Suppl. Mat. D. Their profile $h_{c}(r)$ is the $O(4)$ solution of the equations of motion characterized by a critical radius $R_{c}$ and wall width $l_{0}$, with $R_{c}>l_{0}\sim m_{h}^{-1}$. Contrary to tachyonic instability studies, no (quantum) fluctuations are introduced. Bubbles are nucleated with a random position and orientation $\theta_{i}$ in the $SU(2)$ space: $\phi=\frac{1}{\sqrt{2}}e^{-i\theta_{i}\sigma_{i}}(0,h_{c}(r))\,.$ The different orientations allow the scalar field to acquire a nonzero winding number, generating a current that ultimately sources the gauge fields (initially set to zero): $D_{\mu}W^{\mu\nu}=2g\,\text{Im}\,(\phi^{\dagger}D^{\mu}\phi)$. After collision, bubbles fill the whole volume with the true vacuum. Defining the bubble radius at collision $R_{\star}=\gamma_{\star}R_{c}$, we estimate the Lorentz factor of the wall at collision, $\gamma_{\star}$, via $L^{3}\sim N_{b}\left(2/{\sqrt{3}}\right)^{3}(4\pi/3)\left(R_{c}\gamma_{\star}\right)^{3}$. Fig. 2-left shows the time evolution of the volume-averaged Higgs norm during the transition, for varying $\gamma_{\star}$. Our analysis applies to EW bubbles that collide in the runaway regime, see Sec. B of the End Matter.

The evolution of the volume-averaged energy density components normalized by the initial total energy $\rho_{\rm tot}^{in}$ is shown in Fig. 2-right ($\rho_{\rm tot}$ is conserved). Most of the potential energy is converted into Higgs kinetic and gradient energy, with a non-negligible fraction transferred to electric and magnetic fields. This picture does not change notably for different values of $\epsilon$ from what shown in Fig. 2.

Non-zero $N_{W}$ and $N_{\rm CS}$ are produced as the bubbles collide (see Fig. 1). Once the phase transition completes, $N_{\rm CS}$ relaxes toward $N_{W}$ in the vacuum configuration; see Fig. 5 in the End Matter for an example. Note that the collision of only two bubbles does not produce $N_{W}$ due to the cylindrical symmetry. The CS variance $\Delta n^{2}_{\rm CS,tach}$ in the original cold baryogenesis from a tachyonic instability with instantaneous quench is shown for comparison in Fig. 2 left.

**Dependence of Chern-Simons rate on $\gamma_{\star}$ and $\epsilon$ – ** We study the dependence of the CS rate, $\Gamma_{\rm CS}$, on the EWPT parameters: i) $\gamma_{\star}$, the average Lorentz-boost factor of the bubble walls at collision, which in a cosmological context is proportional to $\beta/H$, where $\beta^{-1}$ is the phase transition duration and $H$ is the Hubble parameter; ii) the degeneracy parameter $\epsilon$, encoding the potential shape and the dynamics of bubble wall collisions. We perform simulations varying $\gamma_{\star}$ (via changing $L$ and thus $R_{\star}$), and for several values of $\epsilon$, with $N_{b}=10$ and $l_{\star}/dx\gtrsim 2$ points per bubble wall (see Suppl. Mat. D for details). The statistical average in Eq. (2) is done over 30 independent simulations. Our results for $\Gamma_{\rm CS}$ have to be extrapolated to the physical limit of large $\gamma_{\star}$, or small $\beta/H$, for which there is a large hierarchy between the bubble size and the microscopic scale set by the Higgs mass. We are able to test values of $\gamma_{\star}\lesssim 8$, and since $HR_{\star}=\gamma_{\star}HR_{c}=(8\pi)^{1/3}v_{w}(\beta/H)^{-1}$ with $H\sim v^{2}/M_{\rm Pl}$ and $v_{w}\sim 1$ the wall velocity, this corresponds to $\beta/H\gtrsim 10^{15}$, which are much larger than the typical ones for a strong first-order PT, namely $\beta/H\sim 10^{1}-10^{3}$.

In Fig. 3-left we show the time evolution of the CS variance, $\Delta n_{\rm CS}^{2}$, and the CS rate $\Gamma_{\rm CS}$, defined by Eqs. (1), (2), for $\epsilon=0.32$. At a fixed $\epsilon$, the CS variance may be parametrized by

[TABLE]

where $c$ is a constant, the exponent $a$ encodes the possibility of a dependence on $\gamma_{\star}$, either increasing or decreasing, $f$ is the profile observed in Fig. 3. In a timescale set by the completion of the EWPT, ${\cal O}(R_{\star}$), the CS variance has grown to a sizable value $\Delta n^{2}_{\rm CS}/m_{h}^{3}\sim 10^{-6}$. This value corresponds to $\Gamma_{\rm CS}$ being roughly of the same order as the thermal sphaleron rate in the symmetric phase, $\Gamma_{\rm th}$, at $T\sim 50\,\text{GeV}$, active for the whole duration of the transition, see inset of Fig. 3-left. Despite the small range of $\gamma_{\star}$ considered for this value of $\epsilon$, as we discuss next, we see some indication for a dependence compatible with $a=1$ in Eq. (3). Since $\Gamma_{\rm CS}\propto\gamma_{\star}^{a-1}$, this would suggest that the (peak of the) CS rate approaches a constant at $\gamma_{\star}\gg 1$ for this large $\epsilon$, see the inset. The characteristic timescale associated with $B$ violation, when $\Gamma_{\rm CS}\neq 0$ and reaches macroscopic values, is $\sim R_{\star}$ (i.e. fractions of Hubble time), in contrast with the tachyonic instability in original cold baryogenesis where the relevant timescale is the time of equilibration of the Higgs modes, which is microscopic.

In Fig. 3-right we show the asymptotic value of $\Delta n_{\rm CS}^{2}$ after the EWPT completes as a function of $\gamma_{\star}$ for different $\epsilon$, normalized to the variance from a tachyonic phase transition, $\Delta n^{2}_{\rm CS,tach}$, shown in Fig. 2 for comparison. For small $\epsilon\lesssim 0.1$, $\Delta n_{\rm CS}^{2}$ decreases with $\gamma_{\star}$ (or equivalently $R_{\star}$), namely $a<0$ in Eq. (3). For $\epsilon=10^{-4}$ the numerical results are compatible with $a=-1$. This dependence is consistent with CS production being purely a surface effect, which vanishes as $R_{\star}^{-1}$ in the physical limit of macroscopic bubbles (with $N_{\rm CS}$ produced only around the collision points). On the other hand, larger values of $\epsilon\gtrsim 0.2$ suggest a CS variance that becomes constant, $a=0$, or grows, $a>0$, in the physical limit of $m_{h}R_{\star}\to\infty$ and $\gamma_{\star}\rightarrow\infty$, as noticed before for $\epsilon=0.32$. In the tested dynamical range, the CS variance from bubble collisions is smaller than the one from tachyonic instability with instantaneous quench by one or two orders of magnitude. However, we stress again that these results need to be extrapolated to the physical point with $\gamma_{\star}\gg 1$.

The regime $\epsilon\gtrsim 0.2$, where interestingly $N_{\rm CS}$ does not appear to decrease with $\gamma_{\star}$, is challenging to study numerically at $\gamma_{\star}\gg 1$ because $R_{c}$ is parametrically larger than the wall width. Since we only test a relatively small range of $\gamma_{\star}$ and results have still quite large statistical uncertainties, the data does not allow us to reconstruct $\Delta n_{\rm CS}^{2}$ reliably at macroscopic $\gamma_{\star}$. In any case, for $\epsilon\gtrsim 0.3$ the data favors a growing variance $\Delta n_{\rm CS}^{2}\propto\gamma_{\star}^{a}$, with $a\gtrsim 1$. The best-fit values are $a\simeq 1.3(2)$ and $a\simeq 1.8(4)$ for $\epsilon=\{0.32,0.45\}$ respectively. Taken together, these results suggest that $\Gamma_{\rm CS}\propto\gamma_{\star}^{a-1}$ is likely not volume-suppressed at large $\gamma_{\star}$. Stronger evidence for this trend, possibly suggesting the value $a\simeq 1$, appears in the (1+1)D $U(1)$ gauge-field case (see Fig. 6 and the discussion in Sec. C of the End Matter), which probes larger $\epsilon$ and $\gamma_{\star}$. The numerical treatment of such (1+1)D system is discussed in more detail in Suppl. Mat. E. Despite these limitations, to the best of our knowledge, this work represents the current state of the art in real-time (3+1)D simulations of first-order PT in $SU(2)$ gauge theory (simulations of 2-bubble collisions with(out) $U(1)$-gauge field have been carried out up to $\gamma_{\star}\sim 5$ in Ref. Lewicki and Vaskonen (2021) ($\gamma_{\star}=20$ in Ref. Gould et al. (2021)), while multiple bubble collisions without gauge fields were simulated in (3+1)D with $\gamma_{\star}$ as large as $\gamma_{\star}\simeq 6$ in Cutting et al. (2018) and $\gamma_{\star}\simeq 50$ in (1+1)D in Jinno et al. (2019)).

Our main conclusion is that $\epsilon$ is a crucial parameter determining whether $B$ can be produced over a large volume. If the vacua are highly non-degenerate (small $\epsilon$), bubble collision is limited to the surfaces. Instead, for nearly degenerate vacua (large $\epsilon$), walls pass through each other, extending the collision region and the production of $N_{\rm CS}$ far out from the initial point. CS transitions then extend to a sizable fraction of the whole volume due to the subsequent transfer of energy from the wall kinetic energy to other gauge and Higgs field configurations (see third panel of Fig 1).

Chern–Simons variance from the spectrum– The CS number density spectrum $P_{\rm CS}(k)$ provides a more detailed information on the length scales at which the CS number is produced. This is defined by $\langle n_{\rm CS}(k)n_{\rm CS}^{*}(k^{\prime})\rangle=(2\pi)^{3}\delta^{3}(k-k^{\prime})\frac{2\pi^{2}}{k^{3}}P_{\rm CS}(k)$, where $n_{\rm CS}(k)=\int d^{3}xe^{-ikx}n_{\rm CS}(x)$ and $N_{\rm CS}=\int d^{3}x\,n_{\rm CS}$. Since the field fluctuations are uncorrelated beyond the typical bubble separation, we expect $P_{\rm CS}(k)$ to have the white-noise behavior $P_{\rm CS}(k)=Ck^{3}/(2\pi^{2})$ for $k\ll 2\pi/R_{\star}$. From Eqs. (1) and (2) one readily verifies that the coefficient $C$ is in fact the CS variance, $C=\Delta n_{\rm CS}^{2}$, see Suppl. Mat. D.

Fig. 4 shows power spectra for $\epsilon=\{10^{-2},0.45\}$ and $\gamma_{\star}=\{3,6\}$. Most of the CS number is generated at scales of order $m_{h}^{-1}$, since both spectra peak at $k\simeq m_{h}$. However, for $\epsilon=0.45$, the spectra also develop a small IR shoulder at $k\simeq 2\pi/R_{\star}$. This becomes increasingly pronounced as $\gamma_{\star}$ grows, at the expense of the $m_{h}$ peak, which is progressively reduced at larger $\gamma_{\star}$. At lower momenta, the spectrum roughly follows the expected white-noise behavior. These properties of the CS spectrum are much clearer in our (1+1)D simulations, as shown in Fig. 6 and further detailed in Suppl.Mat. E. Taken together, these results suggest that in the $\gamma_{\star}\gg 1$ limit, the CS spectrum peaks at $k\simeq 2\pi/R_{\star}$ for large $\epsilon$; consequently, the coefficient $C=\Delta n_{\rm CS}^{2}$ grows with $\gamma_{\star}$, and the production occurs over macroscopic scales $\sim R_{\star}$. The IR shoulder is instead absent for $\epsilon=10^{-2}$, explaining the different behavior for the CS variance with $\gamma_{\star}$ in Fig. 3-right.

**Estimate for the baryon asymmetry – ** This letter computes a universal effect of the EW theory, which can have implications for a large range of Standard-Model extensions. In contrast, a baryon asymmetry calculation requires the introduction of CP-violating interactions and a model-dependent analysis. We can nevertheless argue how this new effect could lead to the correct baryon asymmetry. In a thermal environment, the baryon number density $n_{b}$ obeys a Boltzmann equation that combines the rate of $B$ violation $\Gamma$, the effect of CP violation in producing a net asymmetry $\mu$, and washout terms $\mathcal{A}$: $\frac{\partial}{\partial t}n_{b}(x,t)=\Gamma\left(\xi-\mathcal{A}\,n_{b}\right),$ where $\xi\equiv\mu/T$ is the effective chemical potential for the baryons. $\Gamma$ and $\mu$ are in general space-time dependent, while we will neglect $\mathcal{A}$ for simplicity. In standard electroweak baryogenesis (EWBG) for instance, the chemical potential is non-zero only in the vicinity of the expanding bubble, while sphalerons are active only in the symmetric phase ahead of the wall, $\Gamma=\Gamma_{\rm th}$. On the other hand, the scenario of cold baryogenesis introduced in Ref. Garcia-Bellido et al. (1999) does not rely on a first-order PT for driving the system out of equilibrium, but rather on the tachyonic instability at the top of the SM Higgs potential. This dynamics can be considered homogeneous in space, at least over scales that are much larger than the typical size of the SM textures Garcia-Bellido et al. (2004). The baryon asymmetry can be estimated by referring to a volume-averaged sphaleron rate and chemical potential. The former is extracted from numerical simulations via the diffusion of the CS number, Eq. (1). CP violation is usually introduced via the following dimension-six operator Garcia-Bellido et al. (1999, 2004); Tranberg and Smit (2006) (see however  Tranberg et al. (2010) and  Servant (2014) for other choices), $\mathcal{O}=\frac{g^{2}}{32\pi^{2}\Lambda^{2}}\phi^{\dagger}\phi\,W_{\mu\nu}^{a}\tilde{W}^{\mu\nu\,a}$, where $\Lambda$ corresponds to the new physics scale at which $\mathcal{O}$ is generated. Current bounds on $\Lambda$ are dominated by those from the electron electric dipole moment, see e.g. Cirigliano et al. (2019), and lead to $\Lambda\gtrsim 6$ TeV. $\mathcal{O}$ generates an effective chemical potential $\mu=N_{F}^{-1}\frac{1}{\Lambda^{2}}\frac{d}{dt}\langle\phi^{\dagger}\phi\rangle$, as can be seen by integrating by parts $\int d^{4}x\,\mathcal{O}$ and taking volume averages, $\langle\dots\rangle$, see e.g. Dine et al. (1992); Garcia-Bellido et al. (1999).

The scenario considered in this paper is somewhat in between EWBG and the original cold baryogenesis setup, as it relies on a first-order EWPT but also on the subsequent dynamics leading to a re–equilibration of the Higgs and gauge fields, as advocated in Konstandin and Servant (2011b); Servant (2014). We may estimate the baryon asymmetry following the approach adopted in cold baryogenesis. We stress however that this estimate is based on volume-average quantities, and it may differ significantly in case the CS transitions are localized around the (collided) walls. By considering the same source of CP violation from the operator $\mathcal{O}$, we can derive the effective chemical potential by referring to Fig. 2–left and introducing the following parametrization for the Higgs average variance: $\langle|\phi|^{2}\rangle=\frac{v^{2}}{2}f(b\,t/R_{\star})$, where $b\sim\mathcal{O}(1)$ accounts for a possibly different time scale compared to the CS variance in Eq. (3) (in fact, the Higgs norm converges faster than the CS variance). The total baryon asymmetry is then given by $n_{b}=\int dt\,\Gamma_{\rm CS}(t){\mu(t)}/{T_{\rm eff}(t)}$, where $T_{\rm eff}=\Gamma_{\rm CS}^{1/4}/\alpha_{w}$. We then obtain for the baryon yield to entropy density ratio ${n_{b}}/{s}\simeq 10^{-10}\gamma_{\star}^{3(a-1)/4}\left(\frac{30\,{\rm GeV}}{T_{\rm rh}}\right)^{3}\left(\frac{10\,{\rm TeV}}{\Lambda}\right)^{2}\,$, where $T_{\rm rh}$ is the reheat temperature after the EWPT. In deriving this expression we have taken $R_{c}=10\,m_{h}^{-1}$ as a typical value, and used that $\frac{b}{2}\int_{0}^{\infty}dxf^{\prime}(x)f^{\prime}(bx)\approx 0.1$. This estimate implies that $a\geq 1$ is required for successful baryogenesis (i.e. $n_{b}/s\gtrsim 10^{-10}$); otherwise, $n_{b}/s$ is suppressed at large $\gamma_{\star}$. However, $a\geq 0$ may suffice if the CS transitions occur locally, following the motion of the walls after collision, as suggested by the (1+1)D case in Fig. 6 (bottom panel). In this case the previous estimate based on volume–average quantities would not apply, and the convolution of a local CS rate with the regions where CP violation occurs may still lead to unsuppressed baryon asymmetry, similarly to EWBG. The analysis of the (3+1)D system including a dynamical CP-violating operator (e.g., $\mathcal{O}$) is however needed to ultimately establish the viable parameter space (e.g., $\epsilon$ and $\Lambda$) for successful cold baryogenesis. For such analysis, the key quantity will be the (non–zero) average $\langle N_{\rm CS}\rangle$, which is directly related to baryon number, rather than the variance $\langle N_{\rm CS}^{2}\rangle$ as considered here.

Conclusion– SM Higgs bubble collisions provide a new sizeable source of Chern–Simons number. We calculated from first-principles (3+1)D real-time simulations the corresponding CS rate in the runaway regime and showed its sensitivity on the shape of the scalar potential driving the first-order EWPT. For $\epsilon\gtrsim 0.3$, there is no volume suppression. The next step will be to implement CP violation on the lattice to evaluate precisely the net baryon number. This may provide an alternative realization of electroweak baryogenesis in the large bubble wall velocity limit, which does not rely on the charge transport mechanism. More generally, it will be important to investigate the impact of this potentially new washout effect on the standard EW baryogenesis mechanism. An interesting extension is to include the $U(1)$ gauge field and calculate magnetic-field production in the spirit of e.g. Vachaspati and Brandenburg (2025); Di et al. (2021, 2024, 2025). These simulations can also predict the resulting gravitational-wave background.

Acknowledgements– We thank Anders Tranberg and Hyungjin Kim for discussions. This research was supported in part through the Maxwell computational resources operated at Deutsches Elektronen-Synchrotron DESY, Hamburg, Germany. This work, and MC in particular, is supported by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy—EXC 2121 “Quantum Universe”—390833306. NB acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under the DFG Emmy Noether Grant No. PI 1933/1-1. The work of MG is supported by the Alexander von Humboldt foundation. GS thanks the theory group at Berkeley for hospitality in the very last stages of completion of this work as well as the European Union’s Horizon Europe programme under Marie Skłodowska-Curie Actions – Staff Exchanges (SE) grant agreement No 101086085 -ASYMMETRY.

*End Matter *

Appendix A Winding and CS number

In Fig. 5 we show the evolution of the winding number and the CS number during the collision of three (left panel) and two (right panel) runaway bubbles. The two quantities converge at the end of the EWPT, as expected for vacuum configurations where indeed $N_{\rm CS}=N_{W}$. This represents a non-trivial cross check of our numerical simulations, and shows how winding configurations play a crucial role in triggering CS transitions. Bubble collisions initially generate a non-vanishing $N_{W}$, after which $N_{\rm CS}$ asymptotes toward $N_{W}$. However, $N_{\rm CS}$ does not coincide exactly with $N_{W}$, since part of the Chern–Simons number can be carried away by propagating $W$-boson waves. Conversely, we have checked that the collision of two bubbles leads to identically vanishing $N_{W}$ due to the axial symmetry, and we have consistently observed that the total CS number is $N_{\rm CS}\ll 1$ oscillating around zero. This behavior can actually be seen in the first snapshot of Fig. 1, where a two–bubble collision produces approximately equal and opposite CS number density.

Appendix B Runaway bubbles

As mentioned in the main text, our simulations apply to electroweak bubbles that collide in the runaway regime. For this case to be physically realized, the vacuum pressure given by the zero-temperature potential needs to overcome the friction from particles gaining mass across the wall in the ballistic regime Bodeker and Moore (2009). Moreover, the boost factor at collision, $\gamma_{\star}$, needs to be smaller than the one corresponding to the terminal velocity imposed by balancing the next-to-leading order pressure from particle splitting, $\mathcal{P}_{\rm NLO}\sim\alpha_{w}\gamma_{w}m_{W}T_{n}^{3}$ Bodeker and Moore (2017) (see also Gouttenoire et al. (2022); Azatov and Vanvlasselaer (2021); Azatov et al. (2024); Baldes et al. (2024) for recent work), where $T_{n}$ is the nucleation temperature, with $\Delta V\sim v^{4}$, as typical for a supercooled phase transition. Runaway collisions take place then if $\gamma_{w}\sim 4\pi^{2}(v/gT_{n})^{3}$ is larger than $\gamma_{\star}=R_{\star}/R_{c}\sim vR_{\star}$, with $R_{\star}\simeq(8\pi)^{1/3}\beta^{-1}v_{w}$. Taking $H\sim v^{2}/M_{\rm Pl}$ for the Hubble expansion, one obtains $\beta/H\gtrsim 0.1\,(T_{n}/{\rm MeV})^{3}$. This condition may be satisfied for low nucleation temperatures or large values of $\beta/H$, i.e. for fast phase transitions. For instance, $T_{n}\sim 10$ MeV requires $\beta/H\sim 10^{2}$, while $T_{n}\in[0.1-1]$ GeV leads to $\beta/H\in[10^{5}-10^{8}]$. The latter (large) values of $\beta/H$ could in fact be realized in scenarios with a seeded electroweak phase transition as explored in Refs. Blasi and Mariotti (2022); Agrawal et al. (2024); Blasi et al. (2023), where $\beta/H$ is effectively replaced by the number of defects in the Hubble volume at the time of bubble nucleation, $\xi$, which can be as large as $\xi\sim 10^{8}$ for a domain wall network in the friction–dominated regime at the electroweak scale.

Appendix C Results from (1+1)D simulations of Higgs bubble collisions

To explore large collision boost factors, $\gamma_{\star}\gtrsim\mathcal{O}(10)$, desirable for extrapolation to the physical point, as well as scalar potentials with $\epsilon\gtrsim 0.1$ that produce large critical bubbles, we performed simulations of runaway bubbles in (1+1)D. These significantly reduce computational cost compared to the (3+1)D simulations presented in the main text, allowing us to extend the simulation box more efficiently. For our (1+1)D simulations we have considered the abelian Higgs model with a spontaneously broken $U(1)$ gauge symmetry already adopted in Ref. Garcia-Bellido et al. (1999), as this captures the main topological properties of the SM in (3+1)D. The Higgs potential is modified in order to have a first order phase transition, and we have implemented a scalar potential analogue to the one used in the (3+1)D simulations, as presented in Suppl. Mat. D.

Our main results are displayed in Fig. 6. As shown in the top-left panel, the CS power spectrum exhibits markedly different IR behavior depending on $\epsilon$. For $\epsilon=0.32$ (solid lines), a peak appears at $k\sim 2\pi/R_{\star}$, which shifts toward the IR as collisions go from $\gamma_{\star}=32$ to $\gamma_{\star}=64$, consistently tracking $2\pi/R_{\star}$ (recall that $R_{\star}/R_{c}=\gamma_{\star}$). Since the CS variance is proportional to the amplitude of the IR tail at $k\ll 2\pi/R_{\star}$, which scales as $\propto k$ by causality in (1+1)D, the CS variance is observed to increase with $\gamma_{\star}$ (top-right panel of Fig. 6). In contrast, for $\epsilon=0.01$, the IR power spectrum does not increase with $\gamma_{\star}$, leading to a constant or decreasing CS variance. As in (3+1)D, a smaller and smaller peak appears at $k\sim m_{h}$, but it does not contribute to the CS variance.

Despite the clear differences between the (1+1)D and (3+1)D systems, the Higgs potential shape, encoded in $\epsilon$, strongly influences CS production in both cases. In (1+1)D, the variance is significantly larger than that from the tachyonic instability and grows approximately linearly with $\gamma_{\star}$, as indicated by the $a=1$ line in the top-right panel of Fig. 6 (using the same parameterization of Eq. (3) for the CS variance). Whether such linear growth occurs in (3+1)D remains an open question.

Finally, we comment on the spatial dependence of CS transitions in this lower-dimensional setup, as seen in the bottom panel of Fig. 6. CS number appears to be predominantly produced along the trajectories of collided bubble walls and is not reprocessed afterward. This corresponds to straight lines along the time direction with the same value of the CS number density, $n_{\rm CS}$. Since the collided walls eventually sweep an $\mathcal{O}(1)$ fraction of the volume, CS production is not purely a surface effect. The (3+1)D analogue of Fig. 6-bottom panel, namely Fig. 1, shows a less clear indication for the spatial dependence of the CS production (i.e., whether it follows precisely the collided walls), but we expect that similar dynamics may take place in (3+1)D. Additionally, for $\epsilon=0.32$ (left panel), regions with sizable CS number are much wider than for $\epsilon=0.01$ (right panel). For $\epsilon=0.32$ one can also notice the re-formation of pockets of the false vacuum after the initial collision, as seen for the left-most collision in bottom panel of Fig. 6 (left).

*Supplemental Material *

Appendix D Setup and details on the numerical simulation

We only consider the $SU(2)$ gauge group and set $g^{\prime}=0$. The Higgs field $\phi$ is an $SU(2)$ doublet. During the phase transition the Higgs $\phi$ and the $SU(2)$ gauge bosons $W_{\mu}^{a}$ follow their classical equations of motion (EoM), derived from the bosonic part of the SM Lagrangian:

[TABLE]

Here, $D_{\mu}\phi=\partial_{\mu}\phi-igW_{\mu}\phi$ is the $SU(2)$-covariant derivative with gauge coupling $g$, $W_{\mu}=W_{\mu}^{a}T_{a}$ with $T_{a}=\sigma_{a}/2$, and $W_{\mu\nu}=\partial_{\mu}W_{\nu}-\partial_{\nu}W_{\mu}-ig[W_{\mu},W_{\nu}]$ is the $W$ boson field strength.

As a single-field toy model that induces a first-order electroweak phase transition (EWPT), we consider a sextic potential – i.e., we retain only the terms shown in Eq. (4). The potential is minimized at the vacuum expectation value $|\phi|^{2}=v^{2}/2$, and leads to the canonically normalized Higgs $h=\sqrt{2}|\phi|$ with mass $m_{h}$. The mass of the $W$ boson is $m_{W}=gv/2$.

The regime of interest is characterized by a low cutoff $\Lambda$, i.e.,

[TABLE]

where the potential develops an additional local minimum at the origin $h=0$ even at zero temperature, $T=0$. For $\Lambda$ larger than the upper limit in Eq. (5), the origin is a maximum, while for $\Lambda$ smaller than the lower limit, $\sim 484$ GeV, the origin becomes the global minimum. The potential’s shape depends only on the dimensionless cutoff $\Lambda^{2}m_{h}^{2}/v^{4}$. The potential difference $\Delta V$ between vacua and the degeneracy parameter $\epsilon$ that determines the barrier height with respect to $\Delta V$ are

[TABLE]

$\epsilon$ takes value within $1>\epsilon>0$ for $\Lambda$ in Eq. (5), with $\epsilon=1$ and [math] corresponding to degenerate vacua and no barrier, respectively. In Fig. 7 we show $V$ for different $\Lambda$ and related $\epsilon$, as well as $\Lambda^{2}m_{h}^{2}/v^{4}\to\infty$, i.e. the SM. We will use $\epsilon$ or $\Lambda$ interchangeably according to Eq. (6).

The EWPT with the potential in Eq. (4) has been studied in several works, see e.g. Refs. Grojean et al. (2005); Delaunay et al. (2008); de Vries et al. (2018). The finite-temperature corrections to $V$ are non-trivial and phase transition occurs via thermal tunneling and bubble nucleation. In principle, the transition completes only for sufficiently large cutoffs, $\Lambda\gtrsim 550$ GeV. The nucleation temperature in this case is relatively low, $T_{n}\gtrsim 50$ GeV, and the phase transition strength is $\beta/H\gtrsim 100$ Grojean et al. (2005), with friction likely playing a significant role on the bubble expansion (see Sec. B in the End Matter). For smaller $\Lambda$, the false vacuum is too deep resulting in a too low nucleation rate, and the Universe inflates. Instead, for $\Lambda\gtrsim 1$ TeV the thermal transition is purely second order. We stress that the potential in Eq. (4) is a technical trick to simplify the simulation. A more realistic modeling of the first-order EWPT that we plan to investigate in the future, will be to include a second scalar field $\sigma$ coupling to the Higgs. A particularly natural way to obtain a strongly supercooled phase transition is to have a nearly conformal potential in the $\sigma$ direction while the potential can remain quite standard in the Higgs direction, see e.g. Bruggisser et al. (2023b). In this case, we can expect qualitatively different dependence of the CS production on the global shape of the scalar potential than the one obtained for the simplified single-scalar-field model.

In the following, we neglect that the transition occurs at a finite temperature and instead analyze bubble nucleation, collision, gauge boson production and associated effective sphaleron rate assuming the zero-temperature potential in Eq. (4). In this regard, let us however mention a possible realization where bubble nucleation is actually compatible with a large potential barrier at $T=0$, corresponding to $\epsilon\sim\mathcal{O}(1)$. While this barrier would typically prevent a fast enough nucleation rate, the presence of topological defects at the time of the phase transition can exponentially enhance the rate in their vicinity, as shown in Refs. Blasi and Mariotti (2022); Agrawal et al. (2024); Blasi et al. (2023) for the simplest extension of the SM featuring a singlet scalar field. While the potential barrier around the defects can be effectively smaller and lead to successful nucleation, the dynamics of bubbles at collision will be instead controlled by the actual shape of the potential with $\epsilon\sim\mathcal{O}(1)$, which allows for a larger baryon asymmetry according to the results of our work. This would reconcile efficient bubble nucleation with the optimal type of wall-wall collisions for the production of CS number.

D.1 Initial conditions and critical bubble

For transitions with strengths $\beta/H\gtrsim 1$, the Hubble expansion rate is comparable to or slower than the bubble expansion rate. As a result, the effects of cosmic expansion on the bubble evolution, collision and Chern–Simons number production are expected to be only of order one or smaller. We thus study the EWPT in flat spacetime. To a good approximation the nucleated bubbles are critical Higgs configurations $h_{c}(r)$. These are the $O(4)$ solutions of the Euclidean action’s equation of motion with potential $V$ in Eq. (4), i.e.

[TABLE]

with boundary condition $h^{\prime}(0)=0$ and $h(r\to\infty)=0$. One has $h_{c}(0)\simeq v$ so that they interpolate between the true and false vacua. We define their typical radius $R_{c}$ as $h_{c}(R_{c})=h_{c}(0)/2$. The left panel in Fig. 8 shows the bubble profile for different values of $\epsilon$. The bubble radius $R_{c}$ is of order $m_{h}^{-1}$, however as $\epsilon$ tends to $1$ and the vacua at $h=0$ and $h=v$ become degenerate, $R_{c}$ becomes parametrically larger than $m_{h}^{-1}$, see right panel of Fig. 8; indeed, $R_{c}=3S_{1}/\Delta V$ where $S_{1}$ is the bubble tension. This means that lower $\gamma$ factors can be probed for a fixed lattice spacing. The thickness of critical bubble wall, $l_{c}$, is also set parametrically by $m_{h}^{-1}$, and $l_{c}\sim 3m_{h}^{-1}$ for the values of $\epsilon<1$ shown in Fig. 8.

Bubbles are expected to nucleate at random spatial positions with random orientations in the Higgs vacuum manifold $SU(2)\simeq S^{3}$, i.e.

[TABLE]

where $\theta_{i}$ are three random angles, different for each bubble. To a good approximation, all the bubbles are nucleated at the same time $t=0$. The scalar field associated to the initial condition is simply the sum of single bubble configurations $\phi_{b}$ at rest. The gauge field and its derivatives vanish, $W_{\mu}=W_{\mu\nu}=0$, consistently with the Gauss constraint (see below).

D.2 Higgs and gauge fields evolution

After nucleation, the bubble radius increases as $R(t)=\sqrt{t^{2}+R_{c}^{2}}$. The wall velocity is $v_{w}=\dot{R}$ and its $\gamma$ factor is $\gamma=(1-v_{w}^{2})^{-1/2}$, so that $R=\gamma R_{c}$ and the wall quickly becomes relativistic. At the same time, the effective wall thickness $l_{w}$ decreases as $l_{w}=l_{c}/\gamma$, and gets smaller than $m_{h}^{-1}$. As the bubbles collide, the evolution turns nonlinear following the EoM from Eq. (4),

[TABLE]

We solve these numerically on a cubic lattice with length $L$ and $N_{x}^{3}=500^{3}-1500^{3}$ points, and periodic boundary conditions. Upon redefining time and distances as $(t,x)\to m_{h}^{-1}(t,x)$, and fields as $\phi\to v\phi$ and $W_{\mu}\to m_{h}W_{\mu}/g$, Eqs. (9) and (10) depend only on $\epsilon$ and the ratio $m_{W}^{2}/m_{h}^{2}$, fixed to 0.4 based on the SM parameters. Although $m_{h}/v$ does not enter the EoM, we will use $m_{h}^{2}/v^{2}=1/2$ when showing observables.

The initial conditions have $N_{b}$ bubbles, so that their radius at collision is approximately $R_{\star}=\frac{\sqrt{3}}{2}L\big(\frac{3}{4\pi N_{b}}\big)^{1/3}$ and correspondingly $\gamma_{\star}=R_{\star}/R_{c}$. The physical input parameters are $\epsilon$ and $\gamma_{\star}$ (or equivalently $R_{\star}$). During the evolution we calculate the volume-average Higgs variance $\langle h^{2}\rangle\equiv L^{-3}\int d^{3}x\,h^{2}(x)$, as well as the average components of the total energy density $\rho_{\rm tot}=T_{00}$, which is conserved:

[TABLE]

where $\rho_{E}={E_{i}^{a}}^{2}/2$ and $\rho_{B}={B_{i}^{a}}^{2}/2$ are the energy densities of the electric and magnetic fields, $E_{i}^{a}=W_{0i}^{a}$ and $B_{i}^{a}=\epsilon_{ijk}W_{jk}^{a}/2$. We also calculate the Higgs winding number in the volume $L^{3}$, defined in the main text.

D.3 Link-plaquette discretization, numerical algorithms and observables

Eqs. (9) and (10) are discretized with standard link-plaquette techniques, which preserve gauge invariance at finite lattice spacing $dx=L/N_{x}$, and evolved via a second-order-in-time staggered leapfrog algorithm, with time-step $dt$. For more details see Ref. Figueroa et al. (2021), which we follow.

Links and plaquettes. We define the unitary link field matrix $U_{\mu}(x)\equiv e^{-igdx^{\mu}W_{\mu}(x)}$ going from the space-time point $x+\hat{\mu}$ to $x$, with $dx^{0}=dt$ and $dx^{i}=dx$. Although its argument is $x$, $U_{\mu}$ is to be understood as residing at the midpoint between these two points, namely at $x+\hat{\mu}/2$. Instead, $U_{\mu}^{\dagger}$ goes in the opposite direction, i.e. from $x$ to $x+\hat{\mu}$. We denote link field $U_{\mu}$ at the point $x$ shifted by $\hat{\nu}$ as $U_{\mu}(x+\hat{\nu})\equiv U_{\mu,+\nu}$. The matrix

[TABLE]

represents a plaquette centered at a point, with the factors corresponding to its edges ordered counter-clockwise. The second equality can be readily verified by substituting the definition of $U_{\mu}$. Notice that $U^{\dagger}_{\mu\nu}=U_{\nu\mu}$. From above, the field strength reads

[TABLE]

For simplicity, we define $U_{-\mu}\equiv U_{\mu,-\mu}^{\dagger}$, corresponding to the link field entering in $x$ from $x-\hat{\mu}$. Similarly, $\phi_{\pm\mu}\equiv\phi(x\pm\hat{\mu})$ is the scalar field $\phi$ evaluated at $x\pm\hat{\mu}$.

We solve the EoM in the temporal gauge $W_{0}^{a}=0$, so that $U_{0}=\mathbf{1}$, and $U_{i}=e^{-igdxW_{i}}$ are the only relevant link fields. The field’s convariant derivative appearing in Eq. (9) can be discretized either in the forward or backward directions as

[TABLE]

and should be thought as evaluated at $x\pm\hat{\mu}/2$ (as $U_{\pm\mu}$ are). A generic unitary matrix can be expressed in the basis $\bar{\sigma}_{\nu}\equiv(\mathbf{1},i\vec{\sigma})$ as $U_{i}=\sum_{\nu=0}^{3}c_{i\nu}\bar{\sigma}_{\nu}$, where the coefficients $c_{i\nu}$ are subject to the normalization $\sum_{\nu=0}^{3}c_{i\nu}^{2}=1$.

Evolution equations. We evolve the 4 real scalar fields $\varphi_{i}$, defined via the complex doublet $\phi=(\varphi_{1}+i\varphi_{2},\;\varphi_{3}+i\varphi_{4}),$ and the 12 components $c_{i\nu}$ of the link variables $U_{i}$. Along with these, we also evolve the 4+9 components of the field’s conjugate momenta $\pi_{\phi}\equiv\dot{\phi}$ and $\pi_{W,i}^{a}\equiv W_{0i}^{a}=\dot{W}_{i}^{a}$. The last equality follows from $W_{0}^{a}=0$, so that $D_{0}W_{i}^{a}=\dot{W}_{i}^{a}$. We use a staggered leapfrog evolution scheme, where at each time step the field variables $\phi$ and $U_{i}$ are evaluated at time $t$, while the momenta $\pi_{\phi}$ and $\pi_{W,i}^{a}$ at time $t+dt/2$.

• •

The evolution of $\phi(t-dt)\to\phi(t)$ and $\pi_{\phi}(t-dt/2)\to\pi_{\phi}(t+dt/2)$ is obtained expanding $\dot{\phi}(t-dt/2)$ as $(\phi(t)-\phi(t-dt))/dt+O(dt)$, etc.:

[TABLE]

$\ddot{\phi}(t)$ in Eq. (16) is then replaced with the discretized version of Eq. (9):

[TABLE]

where right hand side (RHS) is calculated at time $t$. The second covariant derivative $D_{i}^{2}\phi$ has been evaluated as $D_{i}^{+}D_{i}^{-}\phi$.

• •

The evolution of $U_{i}$ (and $\pi_{W,i}^{a}$) is slightly more complicated. For the former, we evaluate $\dot{U}_{i}=-igdx\pi_{W,i}U_{i}$ with $\pi_{W,i}\equiv\pi_{W,i}^{a}T_{a}$ at the time $t-dt/2$. In the RHS, we then substitute $\dot{U}_{i}(t-dt/2)=(U_{i}(t)-U_{i}(t-dt))/dt+O(dt)$ and in the left hand side (LHS) we approximate $U_{i}(t-dt/2)$ with the time average $(U_{i}(t)+U_{i}(t-dt))/2$. This leads to the evolution rule:

[TABLE]

We simplify the product of the square brackets to $[1-igdtdx/2\,\pi_{W,i}(t-dt/2)]^{2}$, since the difference is ${O}(dt^{2})$.

• •

Similarly to Eq. (16), the evolution of $\pi_{W,i}=W_{0i}=\dot{W}_{i}$ is given by

[TABLE]

where $\dot{\pi}_{W,i}^{a}(t)$ follows from the gauge field’s EoM, Eq. (10), evaluated for $\mu=i$:

[TABLE]

Specifically, $\dot{\pi}_{W,i}^{a}$ in Eq. (19) is substituted with the discretized version of Eq. (20):

[TABLE]

where $U_{ji,-j}$ denotes $U_{ji}$ at $x-\hat{j}$. The first term in Eq. (21) is the backward difference approximation of the gauge covariant derivative $D_{j}W_{ji}$. To derive it, we used that $W_{\mu\nu}$ transforms in the adjoint representation, along with Eq. (13), noting that both terms in the latter contribute equally. The second term in Eq. (21) originates from $J_{a}^{\mu}$ and is derived using $D_{i}^{+}\phi$ from Eq. (14), noting that only the first term contributes since $\phi^{\dagger}T_{a}\phi$ is real.

A time step consisting of Eqs. (15)+(18), followed by Eqs. (16)+(19), defines a fully algebraic update rule that evolves initial data (given by $\phi$ and $U$ at time $t-dt$, and $\pi_{\phi}$ and $\pi_{W,i}$ at time $t-dt/2$) forward by a time step $dt$.

Note that this algorithm requires the conjugate momenta to be evaluated half a time step ($dt/2$) ahead of the fields. However, our initial conditions are specified at $t=0$: the field $\phi(t=0)$ is given by a superposition of bounce solutions $\phi_{b}(x)$ as defined in Eq. (8), and the initial momenta and gauge fields vanish, i.e. $\pi_{\phi}(t=0)=0$, $W_{i}(t=0)=0$, and $\pi_{W,i}(t=0)=0$. To obtain the values $\pi_{\phi}(t=dt/2)$ and $\pi_{W,i}(t=dt/2)$ required to initialize the evolution, we apply Eqs. (16) and (19) with the RHSs evaluated at $t=0$ and using $dt\to dt/2$. This corresponds to a single forward step of $O(dt)$.

We work with the dimensionless variables and fields obtained via the rescalings defined below Eq. (10). In this way, only the ratio $m_{W}^{2}/m_{h}^{2}$ and $\epsilon$ enter the evolution equations Eqs. (15)+(18) and (16)+(19). We express the evolution equations in terms of $\varphi_{i}$ and $c_{i\nu}$.

Gauss constraint. The $\mu=0$ component of the gauge field EoM, which we have not yet used,

[TABLE]

is a constraint equation – commonly referred to as the Gauss constraint – since it is first order in time derivatives of $W_{i}$. This constraint arises from gauge invariance. If it is satisfied by the initial conditions (as it trivially is), it is preserved by the evolution provided the EoM hold.

As the simulation progresses, we check the validity of the discretized Gauss constraint:

[TABLE]

where we used the fact that $W_{\mu\nu}$ transforms in the adjoint representation and approximated the covariant derivative $D_{i}W_{0i}$ using a backward finite difference in space. The conjugate momenta $\pi_{\phi}$ and $\pi_{W,i}$ are evaluated at $t+dt/2$ during the evolution (and not $t$), so to evaluate Eq. (23) one needs to make a half-step backwards for $\pi$ via Eqs. (16) and (19). We estimate the degree of preservation of the Gauss constraint by calculating the volume-averaged difference and sum of the left- and right-hand side of Eq. (23):

[TABLE]

Given the gauge invariance of the discretized system, $\Delta_{g}$ is expected to remain at the level of machine precision – typically around $10^{-7}$ in single precision and $10^{-14}$ in double precision.

Observables

Field variance and energies. The computation of the volume-averaged Higgs field variance, $\langle|\phi|^{2}\rangle$, as well as the volume-averaged Higgs kinetic and potential energies, $\langle\rho_{K}\rangle=\langle|\pi_{\phi}|^{2}\rangle$ and $\langle\rho\rangle=\langle V(\phi)\rangle$, is straightforward. The Higgs covariant gradient energy is obtained as $\langle\rho_{G}\rangle=dx^{-2}\sum_{i}\langle\left|U_{i}\phi_{+i}-\phi\right|^{2}\rangle$. The electric and magnetic field energies are

[TABLE]

where $W_{ij}^{a}=\text{Tr}[iT_{a}(U_{ij}-U_{ji})]/(g\,dx^{2})$, see Eq. (13).

Higgs winding number. The Higgs winding number is

[TABLE]

where $\Phi=(i\sigma^{2}\phi^{*},\phi)/|\phi|^{2}$. To compute $\Phi_{123}$ and $\Phi_{132}$, we first express all the spatial derivatives of $\Phi$ in terms of $\partial_{i}\phi$, which are then discretized as $\partial_{i}\phi=(\phi_{+i}-\phi)/dx$, as shown in Eq. (14).

Note that cylindrically symmetric configurations (e.g., along the $z$-axis) satisfy $\Phi_{ijk}=\Phi_{jik}$ and lead to a vanishing $N_{W}$. This occurs, for example, in the case of two bubbles or three bubbles with centers aligned along a straight line.

Chern–Simons number. The topological term ${\rm Tr}[W_{\mu\nu}\tilde{W}^{\mu\nu}]$ that enters $N_{\rm CS}$ can be written as the total derivative $\frac{g^{2}}{16\pi^{2}}{\rm Tr}[W_{\mu\nu}\tilde{W}^{\mu\nu}]=\partial_{\mu}j^{\mu}_{\rm CS}$, with the Chern–Simons current

[TABLE]

The time component of this is the Chern–Simons number density

[TABLE]

Using that $\int d^{3}x\,\partial_{i}j_{\rm CS}^{i}=0$ under periodic boundary conditions (or assuming $j_{\rm CS}^{i}$ vanishes at spatial infinity), the Chern–Simons number indeed reduces to

[TABLE]

We calculate $N_{\rm CS}$ by integrating $n_{\rm CS}$ in space as above. Using $W_{i}=\frac{i}{2g\,dx}(U_{i}-U_{i}^{\dagger})+O(dx)$ and $W_{ij}$ in Eq. (13), the latter is calculated from the link and plaquette fields as

[TABLE]

Chern-Simons rate. The Chern-Simons rate measures the number of transitions per unit time and volume between vacua with different Chern–Simons numbers. Formally, it is defined by

[TABLE]

where the denominator $L^{3}t$ accounts for the fact that the variance of the Chern–Simons number grows linearly with both spatial volume and time in a process that proceeds throughout all space and persists indefinitely. The brackets represent an ensemble average. Because bubble collisions occur over a timescale of order $R_{\star}$, it is more appropriate to consider the time-dependent sphaleron rate during this transient period:

[TABLE]

In simulations, the $L\to\infty$ limit in Eq. (31) is effectively realized by ensuring $R_{\star}\ll L$. We calculate the ensemble averages in Eq. (31) by averaging $N_{\rm CS}^{2}$ and $N_{\rm CS}$ over several different simulations with statistically similar initial conditions. The Chern–Simons variance $\Delta n_{\rm CS}^{2}$ is computed first, followed by evaluation of its time derivative (this smoothed over times $>m_{h}^{-1})$.

In presenting the results in Fig. 3 in the main text, the statistical errors over $\Delta n_{\rm CS}^{2}$ and $\Gamma_{\rm CS}$ are estimated as errors on the variance as follows. First, the Chern–Simons variance $\Delta n_{{\rm CS},i}^{2}$ per each run $i=1,\dots,N$ (usually $N=30$) is computed exploiting the white-noise behavior of the Chern–Simons power spectrum for small momenta $k\ll 2\pi/R_{\star}$ (see discussion in the main text and in the following paragraph), i.e. fluctuations of $N_{\rm CS}$ are uncorrelated for lengths larger than $R_{\star}$. In this way, the CS variance can also be computed from only one simulation by dividing the whole simulation volume $L^{3}$ in $n^{3}$ smaller boxes of volume $(L/n)^{3}$, such that $L/n\gg R_{\star}$ (we use $n=4$). The results presented in the main text are obtained using this sub-box division technique. Specifically, the central value of $\Delta n_{\rm CS}^{2}$ is given by the average $N^{-1}\sum_{i}\Delta n_{{\rm CS},i}^{2}$, while its error is estimated as the standard deviation of the values $\Delta n_{{\rm CS},i}^{2}$

[TABLE]

Power spectra. The power spectrum $\mathcal{P}_{\varphi}$ of a real field $\varphi$, e.g. the Chern–Simons number density $n_{\rm CS}$, is defined in analogy to $P_{\rm CS}$ and satisfies $\langle\varphi^{2}\rangle=\int\log k\mathcal{P}_{\varphi}(k)$, as can be immediately seen from its definition. We compute $\mathcal{P}_{\varphi}$ by comparing this last equation to

[TABLE]

which leads to

[TABLE]

where $d\Omega_{k}$ is the solid angle. Here, as before and in the main text, we use the notation $|\mathbf{k}|=k$.

For a white-noise power spectrum $\mathcal{P}_{\rm CS}=Ck^{3}/(2\pi^{2})$, one can show that $C=\langle N_{\rm CS}^{2}\rangle/V$, where $V=4\pi R^{3}/3$ is a spherical volume. Indeed

[TABLE]

(Strictly speaking, the integral over $k$ should be a discrete sum, which is well-approximated by an integral in the large-volume limit.) Using the defintion of $\mathcal{P}_{\rm CS}$ one gets

[TABLE]

Considering a white-noise spectrum peaked at $k_{p}$, $\mathcal{P}_{\rm CS}=Ck^{3}/(2\pi^{2})$ for $k<k_{p}$ and $\mathcal{P}_{\rm CS}=0$ otherwise, and using the change of variable $p=kR$, we compute

[TABLE]

In the LHS we assume that the volume is large enough to resolve the white-noise peak, i.e., $k_{p}R\gg 1$, in which case the integral evaluates to $\pi/6$. From Eq. (37) it follows $C=\langle N_{\rm CS}^{2}\rangle/V$. The same result can be obtained by noticing that $N_{\rm CS}$ actually corresponds to the $|k|\rightarrow 0$ modes of $\tilde{n}_{\rm CS}(k)$.

D.4 Systematic uncertainties

Potential sources of systematic uncertainties in our results are: (1) finite spatial resolution, i.e lattice spacing, set by $m_{h}dx$, (2) finite time resolution, set by $m_{h}dt$, and (3) finite volume effects, set by $R_{\star}/L$ or equivalently by the number of bubbles $N_{b}$ in the box. For the main runs we use values of these parameters that are free from these systematics.

(1) Spatial resolution – resolving the wall width. The grid spacing must be sufficiently small to resolve the bubble wall. Since the initial wall width is $l_{c}\sim 3m_{h}^{-1}$, this requires $m_{h}\,dx\lesssim 1$. As the bubble expands, the wall contracts. To resolve the wall width at the time of collision, $l_{\star}=l_{c}/\gamma_{\star}$, the grid must contain few points $l_{\star}/dx$ per contracted wall. This imposes a more stringent condition:

[TABLE]

We do not calculate the boost factor $\gamma_{\star}$ explicitly. Instead, we estimate it as

[TABLE]

where we used that the mean collision radius is approximately $R_{\star}=\frac{\sqrt{3}}{2}L\big(\frac{3}{4\pi N_{b}}\big)^{1/3}$.

We investigate how our results – primarily $N_{\rm CS}$ – depend on the number of points per wall width $l_{\star}/dx$ by running simulations with $N_{b}=3$ bubbles with fixed position and Higgs orientations $\theta_{i}$, corresponding to a fixed $\gamma_{\star}\simeq 3$. We vary the lattice spacing $dx$ (by changing $N_{x}$) while keeping all other parameters fixed. Fig. 9 shows that $l_{\star}/dx\gtrsim 2$ is sufficient for accurately resolving the wall. We adopt these values for our main runs, except in Fig. 4 (showing $P_{\rm CS}$ at $\epsilon=0.45$), where achieving a relatively large $\gamma_{\star}$ required $l_{\star}/dx\simeq 1$.

As an additional check of bubble wall resolution, we examine the conservation of the total average energy density $\langle\rho_{\rm tot}\rangle$ (Eq. (11)). In our run-away bubble scenario, the majority of the energy is carried by the kinetic energy of the bubble walls at collision. Total energy conservation is thus a robust indicator for the validity of the overall evolution. In Fig. 10, we show the variation of the total energy of 3-bubble simulations, for benchmarks $\epsilon=\{0.01,0.45\}$. $\langle\rho_{\rm tot}\rangle$ is conserved better than per-mille level, with the largest loss occurring at the time of the first bubble collision. If the number of points per wall, $l_{\star}/dx$, is less than 2, $\langle\rho_{\rm tot}\rangle$ is however conserved only at the per-mille level.

(2) Time resolution. In each simulation, we fixed the time-step as $dt/dx=0.3$. The errors in both $\langle\rho_{\rm tot}\rangle$ and $N_{\rm CS}$ from not using a smaller $dt$ are negligible compared to those arising from finite lattice spacing.

(3) Finite volume. The Chern–Simons number $N_{\rm CS}$ and the Chern-Simons rate (Eq. (31)) may be impacted by finite volume effects when the box size $L$ is not significantly larger than the largest characteristic scale of field fluctuations, set by the bubble radius at collision $R_{\star}$. In particular, periodic boundary conditions and a limited number of bubbles can introduce systematic biases in the generation of $N_{\rm CS}$. Since two bubbles alone do not produce $N_{\rm CS}$, repeated collisions of the same pair – enabled by periodic boundaries – can suppress $N_{\rm CS}$ production. To avoid these effects, we ensure that the ratio $R_{\star}/L\simeq 0.5\,N_{b}^{-1/3}$ is sufficiently small – that is, we use a large enough number of bubbles, $N_{b}$. Fig. 11 shows that $N_{\rm CS}$ is statistically indistinguishable for $N_{b}=10$ and $N_{b}=30$, suggesting that $R_{\star}/L\lesssim 1/4$ is already close to the infinite volume limit. In our main simulations, we use $N_{b}=10$.

Finally, Fig. 12 shows the Gauss constraint violation $\Delta_{g}$ (Eq. (24)), which is a strong check of the accuracy of the gauge field evolution. We find that $\Delta_{g}$ is conserved at the $10^{-4}$ and $10^{-13}$ levels when using single- and double-precision variables, respectively. For our main simulations, we use single precision, as the associated truncation errors are negligible compared to the systematic uncertainties arising from finite lattice spacing effects.

Appendix E $U(1)$-gauge field in (1+1)$D$

The 3+1-dimensional $SU(2)$ simulations are constrained to relatively small grid sizes (up to $N_{x}^{3}\sim 1500^{3}$), limiting the accessible Lorentz factors at collision to $\gamma_{\star}\lesssim 10$. To gain insight into the behavior of the Chern–Simons rate at higher $\gamma$, we instead study a simpler 1+1-dimensional $U(1)$ model, which still supports texture configurations. This lower-dimensional setup allows for significantly lattice sizes up to $N_{x}\sim 10^{7}$ points, enabling the exploration of regimes with $\gamma_{\star}\sim 200$. This setup has been considered in the original work on cold baryogenesis  Garcia-Bellido et al. (1999).

We consider a $U(1)$ classical gauge theory in 1D (spatial dimension) with gauge field $A_{\mu}$, $\mu=0,1$, and a complex scalar $\psi$ in the fundamental representation. The Lagrangian is

[TABLE]

with field strength $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and $V$ is the same potential as in Eq. (4). In 1D, the mass dimensions are $[\psi]=[A_{\mu}]=[v]=0$, $[g]=[m_{h}]=1$ and $[\Lambda]=-1$, while $\ell$ in Eq. (6) is still dimensionless. The potential difference between vacua $\Delta V$ and the barrier parameter $\epsilon$ are still given by Eq. (6). Around the true vacuum, the physical Higgs is $h\equiv\sqrt{2}|\psi|$ and has a mass $m_{h}$, while the gauge boson mass is $m_{A}=gv$; note the factor of 2 difference compared to the non-Abelian case. The vacuum manifold is $M=U(1)$. Thus, the theory admits texture configurations labeled by the maps $M=U(1)\to S^{1}=\mathbb{R}\cup\{\infty\}$, where $\psi$ is always in the vacuum. As one-dimensional bubbles collide, we may therefore expect analogous dynamics in the generation of winding number and Chern–Simons number. The equations of motion (EoM) are

[TABLE]

Upon redefining $(t,x)\to m_{h}^{-1}(t,x)$, and $\psi\to v\psi$ and $A_{\mu}\to m_{h}A_{\mu}/g$, they still depend only on $\epsilon$ and the ratio $m_{A}^{2}/m_{h}^{2}=(gv/m_{h})^{2}$, fixed to 0.4. We solve them numerically in a periodic 1D grid with length $L$ and $N_{x}\lesssim 10^{6}-10^{7}$ points.

Bubbles in 1D and initial conditions. Higgs bubbles are segments where the field $h$ remains near the true vacuum $h=v$ over a distance $2R$, with $R$ the bubble radius, and is otherwise at $h=0$. The critical bubbles nucleated are the $O(2)$-symmetric bounce solutions of the Euclidean EoM, $h_{c}(x)$,

[TABLE]

with boundary conditions $h^{\prime}(0)=0$ and $h(x\to\infty)=0$. We define their critical radius $R_{c}$ and wall thickness $l_{c}$ in analogy to the three-dimensional case. Both are set parametrically by $m_{h}^{-1}$ and $l_{c}\sim 3m_{h}^{-1}$ for the values of $\epsilon$ we studied. During the expansion, the bubble radius follows $R(t)^{2}-t^{2}=R_{c}^{2}$. Similarly to 3D, the wall velocity is $v_{w}=\dot{R}=\sqrt{1-R_{c}^{2}/R^{2}}$ and the boost factor is $\gamma=R/R_{c}$.

Note that the bounce solution in Eq. (42) only sets $|\psi|$ and not the phase of $\psi$. As initial conditions, we use a superposition of $N_{b}=10$ critical bubbles $\psi_{c}(x)=e^{i\theta}h_{c}(x)/\sqrt{2}$, each placed at a random position in the domain $L$ and assigned a random phase $\theta$, independent for each bubble. In addition, we set $A_{\mu}=F_{\mu\nu}=0$.

Winding and Chern–Simons numbers. The winding and Chern–Simons numbers are

[TABLE]

The Chern–Simons number density is simply $n_{\rm CS}=gA_{1}/(2\pi)$. As mentioned, the theory in Eq. (40) supports texture configurations with a nontrivial $N_{W}$, e.g. $\psi(x)=e^{i\chi(x)}$, with $\chi(x)=\pi\tanh(x/r_{t})$ has $N_{W}=1$. The Chern–Simons variance and rate are

[TABLE]

Similar to 3D, we also compute the Chern–Simons power spectrum $\mathcal{P}_{\rm CS}$ defined by

[TABLE]

As in 3D, for $k<2\pi/R_{\star}$, the spectrum takes a white-noise form, $\mathcal{P}_{\rm CS}(k)=Ck/\pi$, where, as shown by a calculation similar to that in Suppl. Mat. D.3, the coefficient is given by $C=\Delta n_{\rm CS}^{2}$.

E.1 Discretization and numerical evolution

We solve Eq. (41) with the same link-plaquette techniques and staggered Leapfrog algorithm as the 3D case in Suppl. Mat. D.3. Here we highlight only the differences. Although the spatial index is $i=1$, for generality we express the formulas for an arbitrary number of spatial dimensions. By running large grids in (1+1)D, we can use a particularly fine lattice spacing $dx$, ensuring that the number of points per wall width at collision satisfies $l_{\star}/dx\gtrsim 8$. As discussed below, this keeps the main quantities close to the continuum limit.

Links and plaquettes. The link field $V_{\mu}\equiv e^{-igdx^{\mu}A_{\mu}(x)}$ that connects $x+\hat{\mu}$ to $x$ is now a phase, and $V_{\mu}^{*}$ connects the two points in the opposite direction. The plaquette $V_{\mu\nu}$ is

[TABLE]

Notice that $V^{*}_{\mu\nu}=V_{\nu\mu}$. We define $V_{-\mu}\equiv V_{\mu,-\mu}^{*}$, and the forward and backward covariant derivatives are as before $D_{\mu}^{\pm}\phi=\pm\left(V_{\pm\mu}\phi_{\pm\mu}-\phi\right)/dx^{\mu}+O(dx)\,.$

Evolution equations. We use the gauge $A_{0}=0$, leading to $V_{0}=1$ and $V_{i}=e^{-igdxA_{i}}$. We evolve the real and imaginary parts of $\psi$ and of $V_{i}$, as well as of their conjugate momenta $\pi_{\psi}\equiv\dot{\psi}$ and $\pi_{A,i}\equiv F_{0i}=\dot{A}_{i}$. The evolution of $\psi$ is as in Eq. (15):

[TABLE]

where, from the scalar field EoM in Eq. (41),

[TABLE]

As in Eq. (18), from $\dot{V}_{i}=-igdx\pi_{A,i}V_{i}$, it follows the evolution of $V_{i}$:

[TABLE]

Similar to Eq. (19), for the evolution of $\pi_{A,i}$ we have

[TABLE]

where $\dot{F}_{0i}=\partial_{j}F_{ji}+2g{\rm Tr}[\psi^{*}\partial_{j}\psi]$ follows from the EoM in Eq. (41). Its discretized version is

[TABLE]

Gauss constraint. The Gauss constraint $\partial_{i}F_{0i}=\partial_{i}\pi_{A,i}=2g{\rm Im}[\psi^{*}\dot{\psi}]$ is discretized as:

[TABLE]

where as in 3D we make a step backwards for $\pi_{\psi}$ and $\pi_{A,i}$. $\Delta_{g}$ in Eq. (23) is again used to evaluate the preservation of the Gauss constraint. We use double-precision variables and observe conservation of $\Delta_{g}$ at the $10^{-10}$ level.

Observables. The average energy is $\langle\rho\rangle=\int\frac{d^{3}x}{V}\rho$ and

[TABLE]

with $E_{i}=F_{0i}=\pi_{A,i}$, $B_{i}=\frac{1}{2}\epsilon^{ijk}F_{jk}$. Note that $B_{i}=0$ in 1D and $|D_{i}\psi|^{2}=dx^{-2}|V_{i}\psi_{+i}-\psi|^{2}$. We compute the 1D winding and Chern–Simons numbers of Eq. (43) as

[TABLE]

Here, $x_{n}=n\,dx$, and the sums run over all grid points. Note that $\mathrm{Arg}[V_{1}]=\tan^{-1}(\mathrm{Im}[V_{1}]/\mathrm{Re}[V_{1}])\simeq\mathrm{Im}[V_{1}]$ to a good approximation. With our choice of $dx$, the total energy is conserved to the $10^{-6}$ level.