Simulations of a bubble wall interacting with an electroweak plasma

We perform large-scale real-time simulations of a bubble wall sweeping through an out-of-equilibrium plasma. The scenario we have in mind is the electroweak phase transition, which may be first order in extensions of the Standard Model, and produce such bubbles. The process may be responsible for baryogenesis and can generate a background of primordial cosmological gravitational waves. We study thermodynamic features of the plasma near the advancing wall, the generation of Chern-Simons number/Higgs winding number and consider the potential for CP-violation at the wall generating a baryon asymmetry. A number of technical details necessary for a proper numerical implementation are developed.


Introduction
A first order phase transition in the early Universe proceeds through the nucleation of bubbles of the low temperature phase inside the background high temperature phase. These bubbles grow and coalesce while interacting with the ambient plasma, thereby generating gravitational waves (see for instance [1]), with a spectrum that may be detectable by the upcoming LISA mission [2]. If the first order phase transition is associated with electroweak spontaneous symmetry breaking, baryogenesis may take place in the presence of CP-violating interactions near the bubble wall [3,4]. A first order electroweak phase transition does not arise in the Minimal Standard Model, but is a fairly generic result of extending the scalar sector with additional fields [5][6][7][8][9][10][11][12][13][14][15][16].
The precise dynamics of the bubble wall is complicated by the interactions with the plasma, and the transport of any currents and energy-momentum away from or through the wall is a highly non-equilibrium phenomenon, ultimately requiring a non-perturbative real-time treatment. One hope to simplify things is that the propagating wall has enough time to reach a steady-state regime, wherein the Higgs profile of the wall and the velocity of the wall are roughly constant. And furthermore, that this steady-state regime persists long enough to dominate any creation of a baryon asymmetry. This is simpler to treat analytically than, say, a fast evolving transient stage prior to such a steady-state or while -1 -(or even after) the walls are colliding, the latter being the primary source of gravitational waves [1].
Ideally, one would wish for a real-time numerical simulation of the bubble nucleation process, then the growth and evolution of bubbles followed by the steady-state regime, the bubble collisions, coalescing of domains of low temperature phase and thermalization of the released latent heat into the plasma. Also, ideally, such simulations would include all the degrees of freedom of the (extension of the) Standard Model, gauge fields, scalar fields and fermions.
The tools for this exist (see for instance [17][18][19][20]), but are numerically much too expensive to simultaneously include all degrees of freedom, with the fermions inherently quantum mechanical, the very large volumes necessary for multiple bubbles to fit, and for the microscopically quite long times involved in nucleation, growth, collision and thermalization. One solution is to split up the problem into smaller, more tractable, components some of which may be treated in or near thermal equilibrium. A huge analytic and numerical effort over decades has been put into different aspects of this problem ( for reviews, see for instance [21,22]).
In this work, we make the following simplifications: We include only the bosonic part of the Standard Model, a single Higgs field and the SU(2)×U(1) gauge field, and we treat their dynamics classically. We assume that the bubble has been nucleated at some earlier stage, and consider a mostly planar Higgs wall. This wall is advancing not as a result of the thermodynamic pressure in the physical transition, but through the driving force of an external current, which we insert by hand. We are then free to choose this current to enforce a wall speed and profile (more or less) of our choosing.
We then simulate this wall sweeping through the plasma of quasi-thermal Higgs and gauge field fluctuations, and may compute various components of the energy-momentum tensor, as well as other field observables, including topological ones (winding number, Chern-Simons number). This will allow us to monitor the development of an out-ofequilibrium envelope near the wall, which we expect is the region of most relevance to baryogenesis. We can follow the creation of Higgs winding number and estimate the possibility of an asymmetry being created. We will also consider the introduction of CP-violation in a particular form, but although we did perform first-principles simulations including such CP-violation, they were inconclusive and are postponed for upcoming work.
The paper is structured as follows: In section 2 we will describe the field theory model, and introduce notation and basic observables. In section 3 we develop our numerical setup, including the initial conditions of the plasma and the driving of the wall. We will do some basic tests of the setup, wall speed and shape. In section 5 we perform simulations focusing on the thermodynamical observables, and the transport of energy-momentum near the wall. In section 6 we instead emphasize topological transitions (Higgs winding), and connect to baryogenesis. In section 7 we comment on improvements, caveats and possible future avenues for such simulations, and briefly discuss the introduction of CP-violation through a bosonic dimension 6 operator. Finally, we conclude.
The simplest extension of the Standard Model viable for electroweak baryogenesis is the 2-Higgs Doublet model (2HDM) [5,6,9,13,14,16], where the Minimal Standard Model is augmented by an additional Higgs doublet field 1 . This new field couples to SU(2) and U(1) gauge fields exactly as the first one, and different models allow for different Yukawa coupling patterns between Higgs and fermion fields (see [23,24] for model descriptions and experimental constraints).
Depending on the shape of the Higgs potential, the quasi-order parameter of the transition (the operator that is close to zero outside the bubble, and large inside) may be the original Higgs field, the new Higgs field or a combination of the two. The bubble wall is then a field profile interpolating between the inside value and the outside value, minimizing the (free) energy.
In this first investigative work, we will stick to an order parameter consisting of a single Higgs field, and postpone the case when the fields mix as we move through the wall. We will ignore SU(3) gluons and to a first approximation neglect the fermion degrees of freedom. Fermions colliding with the wall give an important contribution to the force determining the wall dynamics. But because we are moving the wall by hand, the backreaction on the wall is not our primary concern. Also, the numerical methods available to simulate real-time fermions are prohibitively expensive [25].
This means that we will approximate the electroweak plasma by the SU(2) and U(1) gauge fields, coupled to a single Higgs doublet. The doublet is fully dynamical, and interacts with the gauge fields but, as we will describe, an external current forces the Higgs to have a wall-like profile and drives it through the plasma with an overall shape and velocity of our choosing.
Our SU(2)-U(1)-Higgs model is then described by the continuum action where φ is the Higgs scalar doublet, W a µ are the SU(2) gauge fields and B µ the U(1) hypercharge gauge field. We have defined and the covariant derivatives where the Higgs field has hypercharge Y = −1/2. We have introduced the parameters

4)
1 A first order phase transition can also be realised by just adding a singlet field, but then insufficient CP-violation is a problem (see for instance [7, 8, 10-12, 14, 15] The term ∆L 2 represents a potential CP-violating term, which we would like to introduce to bias the dynamics. We will return briefly to this in the conclusion.

The wall and initial conditions
In a first order electroweak phase transition, the high temperature global minimum of the effective potential at |φ| = 0 becomes degenerate with a second minimum at 0 < |φ| < v/ √ 2 as the temperature is lowered to some critical T c . Lowering the temperature further, the second minimum takes over as the global minimum and the minimum at |φ| = 0 becomes first a local minimum, and eventually ceases to be a minimum altogether.
If the system starts out in the first minimum |φ| 0, there is a period of metastability while this minimum still exists, but is not the global minimum. This becomes more and more pronounced as the (free) energy difference between the two minima increases with decreasing temperature. At some nucleation temperature, T N , random thermal fluctuations have a sufficient probability of generating a large enough bubble of the low-temperature phase. If the bubbles are large enough, they will dynamically continue to grow rather than collapse back, and the phase transition is triggered.
The walls of the bubble interpolate between the metastable minimum (outside) and the global minimum (inside), with a profile so as to minimize the (free) energy of the wall. In vacuum this profile may be found by shooting methods, or one may find an approximate expression for its shape by invoking a thin wall approximation or a thick wall approximation. Ultimately, one may solve for the profile numerically, also in a thermal plasma (see for instance [17]). It will depend on the supercooling T c −T , the Higgs potential and the thermal fluctuations through the free energy associated with a unit area of wall, the wall tension. If the bubble is not very large relative to the wall thickness, curvature effects may also play a role.
To tune the parameters and the temperature to the exact nucleation point, and to realise this transition in individual classical simulations is difficult. We will instead force a moving wall into our system by adding a current R(x, t), so that the potential for the Higgs field is time and space dependent, Whenever and wherever R = 0, we recover the Standard Model potential, which at zero temperature has a global minimum at |φ| = v/ √ 2. At leading order, thermal fluctuations at some temperature T will contribute a term analogous to R T 2 and raise the global -4 -minimum, and at high enough temperature T > T c , the global minimum is at |φ| = 0. This thermal contribution is the same inside and outside the bubble, and so in a true phase transition, an "outside" and an"inside" is established not because R varies in space, but because the temperature (and hence potential) is so finely tuned that it has two minima and both phases are allowed simultaneously. As we will describe below, we will operate at a lower temperature than the critical one, and instead drive the wall through a choice of R In this way, the potential has a minimum at |φ| = 0 in the region of space we want to be outside the bubble and minimum at a non-zero value of the field in the region we want to be inside. By changing R, we can make the bubble grow, by first having R non-zero in all of space, and then flipping it to zero in an ever expanding region. We emphasize that we do not directly stipulate the wall, but rely on the field dynamics to generate the wall itself in the background of this current. In particular, the Higgs field is free to have physical fluctuations around the overall wall shape. Moreover, R is a gauge singlet, and so introducing does not affect Gauss's law, which is an important consideration in the numerical simulations. It turns out that for the initial conditons we have in mind (to be described shortly) it is important that the zero-temperature mass outside (or prior to the appearance of) the bubble is small. By choosing the current R to be 1 2 m 2 H outside the bubble this mass is zero, and we avoid that the random fluctuations, defined by the Bose-Einstein distribution below, are cut off by a large mass.

Initialization of Higgs and gauge fields
Our simulations emulate a region of space with a thermal plasma in the high-temperature phase, where a bubble wall sweeps through as it expands. In our numerical setup, we have a finite 3-D box, of dimension L x L y L z , where z denotes the direction of motion of the wall, and L z L x , L y . The box has periodic boundaries, a point to which we shall return below.
For the situation prior to the arrival of the bubble wall, we wish to introduce fluctuations in the Higgs and gauge fields that are quasi-thermal , and may provide an ensemble of random classical field realisation. This is in principle straightforward for a free theory, but because the fields interact, the true thermal state is non-linear and must obey Gauss's law. This may in turn be resolved by brute force, for instance by Monte Carlo methods [26], but this is numerically heavy. Since we are not attempting to reproduce an exactly thermal state with a specific temperature, we allow ourselves some short cuts. These may later be ameliorated.
For each of the four real component of the Higgs scalar doublet we set φ a = 0, a = 1, 2, 3, 4. The momentum variables π a = ∂ t φ a are taken to be non-zero, and to satisfy the free-field two-point correlations π a (p)π b (p ) = 2ω p n p (2π) 3 δ ab δ 3 (p − p ). (3.4) The total energy of a free field system is For an equipartitioned system the two terms in the bracket give the same contribution, and if we insist on not initialising φ (the second term), we can mimic the correct initial condition by putting twice as much energy in the momentum component π (the first term). This explains the overall factor of 2 in (3.4). The initial state is then determined by the particle number n p and the frequency ω p , which we take to follow a Bose-Einstein distribution, We have introduced the one-loop thermal mass, m 2 eff (T ), to take into account (at leading order) that the system is interacting [27,28]. We find the thermal mass self-consistently by iterating Eqs. (3.6), (3.7). Since we have chosen the zero-temperature mass to be zero outside the bubble, only the thermal mass enters in the dispersion relation. This is convenient, since it is fairly small and allows for sizable occupation numbers in the IR, while still cutting off the IR divergence.
By choosing this initialisation, where φ a = 0, we ensure that the initial local Higgs charge vanishes, rendering the Gauss law much simpler to solve. Obviously, already one time-step further down the line, both φ a and π a will be non-zero. However, the ensemble does not start out thermal, since only the momentum variables are initialized. Rather quickly, the energy will equipartition to some near-equilibrium state, but the parameter T is no longer the actual physical temperature. We should therefore treat it as a parameter fixing the energy density, which in turn is related to the actual temperature of the system. For more on real-time initialization of Higgs fields on the lattice, see [29].
We initialize the gauge fields in a similar way, but now we set the momentum variables, the electric fields, to zero E a i = 0. The field variables A i are initialized as a free-field-like thermal ensemble, again with an overall factor of 2, so that the total initial energy is approximately correct. This applies to the three SU (2) gauge fields and single U (1) gauge field. In practice, we write where ξ a (p, λ) is complex random number with variance (ξ a (p, λ)) * ξ a (p, λ) = 2 and the two transverse polarization vectors are with r a random vector. The initial particle number again follows the Bose-Enstein distribution, In the symmetric phase, we take all the gauge fields to be massless, and so the zero mode is singular. We set it to zero initially. With this initialization, Gauss's law is explicitly obeyed everywhere on the lattice initially, and the classical equations of motion ensure that this is then also true at all later times. The procedure for preparing our initial condition is described in the following, and is depicted in Fig. 1: • We start by initialising the fields throughout the lattice at a temperature T i (see section 3.1) in the symmetric phase with R(x, t) = 1 2 m 2 H . We will generate a whole ensemble of such classical realisations, and observables quoted below are averages over such an ensemble of real-time simulations. We leave the system to equilibrate for a time m H t thermal = 50, as seen in Fig. 1, top left.
• We then establish the wall by changing the current R(x, t), so that it depends on time and space, but in the z-direction only We will think of this as smoothed-out step functions, and will establish the profile of the wall through the parameter d.
We are creating two walls, one centered at z r and one at z l . This is so that we may have periodic boundary conditions in the z-direction throughout. The x-and y-directions along the wall are also periodic. We may later decide whether both walls move outwards or just one of them. z r − z l determine the size of the initial bubble, prior to expansion. Hence before we make the walls move, part of the volume is already in the broken phase.
The time dependent parameter c(t) changes from 0 to 1 and determines how fast the bubble wall is established. If we do this very fast, the wall is shocked into existence and this may bring in unwanted, and unphysical, transient effects. In practice, we find that choosing is convenient 2 , with m H τ wall = m H τ thermal = 50. This is seen in Fig. 1, top right.
• After the current is established (c(t) = 1), we leave the wall to settle into its shape and for any transients to damp away. We do this for a time m H τ stable = 50. We are then ready for the simulation proper to commence. Fig. 1, bottom left, and bottom right show these phases.

Running the wall
After the walls have been established and have settled down, we may control their motion by further changing the current R. We use the form with starting time, τ move = τ thermal + τ wall + τ stable , taken to be = 150. v is the speed of the wall, and we have included the Lorentz factor γ = 1/ √ 1 − v 2 for the moving (scalar) wall. Note that we may choose to move both walls, but to take maximum advantage of the available lattice volume, we move only one, keeping the other at rest. In Fig. 2 we show the initial wall and walls after a time m H t = 500 with different speeds. Dots indicate the lattice points, and give an impression of the discretization. Clearly, by v = 0.99 our resolution of the wall is no longer satisfactory.  To summarize: On a lattice of size N x = N y = 64 N z = 1000, physical size L i = N i a with lattice spacing a, am H = 0.5, we initialize Higgs and gauge fields using a parameter T , and let it settle. We then grow a pair of walls spaced by z r − z l with a width parameter dm H = 15, and let it settle. We then drive one of those walls through the plasma with a speed parameter v, until we run out of box.
With this many parameters, we have had to make some choices, and although we have tried quite a few combinations to find some sensible values, a complete sweep has not been done. For the largest possible volumes N x , N y , N z , we focus on the dependence on v for a few values of T .

Observables
Given the initial conditions presented above, the simulation is ready to go. We are interested in monitoring a broad set of observables, as we carry out the simulation. We monitor -9 -total Gauss law and total energy over the lattice to check our numerics. But the physically most interesting quantities are space-, and in particular z-dependent. We will therefore introduce a number of quantities averaged over x-y spatial slices, such as the Higgs field This may further be integrated over z, to get average Higgs field over the whole volume φ 2 (t), if required. Similarly, for the observables introduced in the following, we will consider local, slice and global versions as appropriate.

Thermodynamical observables
We may compute diagonal components of the energy-momentum tensor. For the scalar field the energy density is The three diagonal space components are and where the equality of the latter two is ensured by the symmetry of the system. Similarly, for the gauge fields we have while the three diagonal components are and For all of these quantities, we will present results averaged over a slice perpendicular to the z-direction (similar to (4.1)).
It is not unreasonable, as a first approximation, to assume that the plasma acts as a perfect fluid. In that case, the energy-momentum tensor can be written as (4.8) -10 -Under the assumption that the fluid moves only in the z-direction, we can write Inserting this into (4.8), we find so that We see that if the space-like diagonal components are the same, i.e. b = u √ 1−u 2 = 0, then the energy density is simply the 00-entry. We will see below, however, that because of the push of the wall, the zz-component is different from the other two, and there is a net space-dependent fluid velocity, which we extract as indicated in Eq. (4.11).  Electroweak baryogenesis is closely related to the chiral anomaly, relating baryon number to the Chern-Simons number of the underlying gauge field.

Topological observables
For the degrees of freedom included in our simulations, only Chern-Simons number is available to us, one for the SU(2) field and one for the U(1) hypercharge field, 14) and A particularly useful proxy for the baryon asymmetry turns out to be the Higgs winding number where the matrix form of the Higgs field φ is defined by Close to equilibrium and at low temperatures, N w N cs N cs,2 . The first relation follows from minimization of the covariant derivative which is achieved at low energies, when the gauge fields are close to a "pure gauge" vacuum. It is not an identity, and Chern-Simons number and winding number may be very different in a violent non-equilibrium environment. The second relation follows from the field space of the U(1) field, which has a single vacuum at N cs,1 = 0. In contrast, the SU(2) field has infinitely many degenerate vacua, enumerated by integer values of N cs,2 . We will consider Higgs winding number "density" in z, corresponding to the integrand of the expression for N w , but which is then integrated over x and y. Outside the bubble, where φ 0, the winding number density observable is very noisy, and integrating it up as described is demonstrated in Fig. 3. We will therefore choose to show only the winding number inside the bubble for either the whole bubble volume (which grows over time as the bubble expands), or a smaller, constant volume immediately behind the wall.  Looking first at the bottom left plot (T /m H = 0.5, v = 0.25) the positions of the initial walls are around z = 30, and as the bubbles come up in the time interval 50-100 a small spike in energy appears. Then after settling down, the right-hand wall starts moving to the right from time 150. We see that the straight line with speed 0.25 corresponds to the upper edge of the red wedge. With excess energy density coloured red, this means that the energy is deposited in front of the wall, spreading out far ahead of it to create a large region of activity. From this plot we can see the effect of the periodic boundary conditions, with the energy flow from the moving bubble wrapping around the box. This will need to be borne in mind when viewing the profiles of thermodynamic quantities. Furthermore, we note that two energy "beams" from the initial growth of the walls travel in both directions around the lattice. These are an artefact of the lattice implementation , and may be reduced by an even slower wall initialization.

Thermodynamics of the plasma
The simulation ends at a time 1000, which is 850 after the wall started moving, Indeed, in this plot the wall has moved to 850v from the initial wall position to m H z 250.
Increasing speed (moving up in the left-most column of plots), the energy wedge becomes more pronounced, but narrower, as the wall is able to more and more keep up with the released energy. Whereas for v = 0.5 all the energy is still outside the bubble, as we increase the speed to v = 0.75 some of the energy is overtaken by the wall, so that energy is also deposited inside the bubble. For the fastest speed of v = 0.99, the wall and energy profile are very localised and, as indicated earlier, we begin to be wary of discretization and resolution effects. We also remark that the current R, responsible for the bubble, is stopped when it reaches the edge of the simulation. For wall speeds of v = 0.75 and v = 0.99 this happens before the end of the simulation, as can bee seen in the top three rows of Fig. 4.
Moving instead to the right in the array of plots, temperature increases, and we see a similar picture in each case, but with much larger numbers and contrasts (note the difference in colour scale). In particular, for faster speeds, energy density is definitely overtaken by the wall and is deposited both inside and outside the bubble. Another important feature in all these plots is the scaling, or self-similar, behaviour. We can clearly see that the energy density forms a wedge-like shape as the bubble expands, which is expected as the time since nucleation is the only relevant macroscopic variable in the simulation. This means that at late times, thermodynamic profiles are functions of ξ = z/(t − τ move ), as is expected from relativistic hydrodynamics [30].
Finally, in a relativistic fluid made up of massless particles we have an equation of state P = 1 3 ρ, leading to a sound speed (squared) of c 2 s = ∂P ∂ρ = 1 3 . In the high temperature simulations (T = m H ) of Fig. 4, for wall speeds of v = 0.25 and 0.5 , we have included a curve corresponding to this sound speed . As we can see, the outer part of the fluid envelope is propagating at a speed consistent with this sound speed.  In addition to the energy density, we also show the other diagonal components of the energy momentum tensor, Fig. 5. All quantities have been averaged over the x-y plane, and are shown as a function of z at particular times, m H t = 600 for v = 0.25, 0.375, 0.5, 0.625 0.75 and m H t = 500 for v = 0.99. We also show the inferred fluid velocity, assuming a perfect fluid form of the energy-momentum tensor, as described in section 4.1. We must bear in mind that the effects of the periodic boundary conditions seem particularly pertinent for low temperatures, as can be seen from Fig. 4, and so some care is needed when interpreting Fig. 5. The blue dashed line is the external current R (divided by 4R max as measured on the right-hand vertical axes), and this describes the position of the wall as it moves outwards.
Beginning again in the bottom left plot (low temperature, 0.5 m H , low speed, 0.25), we see that all of the observables have a clear jump in magnitude, coinciding with the position of the wall. This reflects the observation in Fig. 4 that at low speeds, the energy is deployed outside the bubble, ahead of the wall. The z-component is somewhat different from the x-and y-components, but the anisotropy is rather weak at these low speeds. This corresponds to a quite low fluid velocity around 0.1.
Picking up speed (v = 0.5, upwards in the column of plots), we see that the picture persists, but with a larger overall amplitude, a much larger anisotropy and hence a much larger fluid velocity. All observables decay smoothly with the distance from the wall, as they interact with the ambient plasma. At a glance, the dumped energy-momentum reaches a distance of about zm H ∼ 150(m H t/600) ∼ 1 4 m H t from the wall, where the factor of m H t/600 comes from the self-similar behaviour, allowing us to extrapolate to later times. The fluid velocity is still much smaller than the wall speed.
At a speed of v = 0.75, the wall is able to overtake some of the thermodynamic envelope, before it has time to thermalize with the plasma. As a result, there is an energy peak, for which the maximum coincides with the position of the wall. The tails of the peak stretch both inside and outside the bubble. The peak is asymmetric, reaching a distance of about zm H ∼ 50(m H t/600) ∼ 1 12 m H t outside the wall and as far as zm H ∼ 200(m H t/600) ∼ 1 3 m H t inside the wall. The anisotropy is even more pronounced than before, because of the high wall speed.
For asymptotically high speeds (v = 0.99), the picture is even more extreme. The wall is now firmly ahead of most of the released energy, and no longer coincides with the maximum of the energy peak. This peak is now quite broad and stretches a distance of at least zm H ∼ 100(m H t/500) ∼ 1 5 m H t inside the bubble. We also see an anisotropy between z and x, y components is a factor of 4, giving a large fluid velocity. Outside the bubble, in contrast, the plasma is fairly quiet.
Proceeding to larger temperatures (second and third column of plots), both the overall relative magnitudes of the observables increase substantially. The qualitative picture is the same, but quantitative differences exist. As an example, for the largest temperature T /m H = 1 and a speed of v = 0.75, the energy peak reaches only about half as far into plasma outside the bubble. This likely means that a baryogenesis mechanism relying on out-of-equilibrium conditions and CP-bias emerging from the bubble wall will have less volume to act.
Also, it seems that although the (maximum) fluid velocity is always smaller than the wall speed, the two are approximately proportional to each other, with the fluid velocity roughly inversely proportional to temperature (Table 1 ).

Dynamics of topological observables
Single configurations Electroweak baryogenesis relies on a CP-asymmetric fermion current biasing near-equilibrium sphaleron processes in the bulk. As the bubble wall collides with the plasma, CP-violating interactions produce a net left/righthanded fermion current through/reflected from the wall into the bubble/back into the symmetric phase plasma in front off the wall. If the sphalerons are sufficiently suppressed in the broken phase inside the bubble, baryon number is only generated outside the bubble in response to the CP-violating currents.
Since we do not include the fermion degrees of freedom, the role of a baryon number is played by Chern-Simons number and Higgs winding number. In order to potentially be able to react to a CP-bias, we need to have sufficiently high temperature for sphaleron processes (and more general processes changing Chern-Simons number and winding number) to readily take place, at least outside the bubble.
In Fig. 6 we show the winding number inside the bubble as a function of time for two different sets of realisations, with two different combinations of temperature and wall speed. The three blue dashed lines denote the set-up stages of the wall (τ thermal + τ wall + τ stable ). The last dashed line corresponds to the wall reaching the end of the lattice and stopping.
The winding number is first recorded as the bubble wall comes up at around m H t = 50, both because winding number is generated, but also since that is when an "inside" of the bubble begins to exist. Once the bubble moves, more and more of the volume is included inside the bubble, and one would expect that the total winding number increases, assuming that there is a random distribution of winding number nuclei distributed throughout the volume. For the cold temperature (left-hand plot) this does not happen, and the integer winding number for each realisation remains at its initial value. This shows that the system is too cold to spontaneously generate winding number nuclei.
For the larger temperature (right-hand plot) we see that there are instances of "flips", where the winding number increases as a symmetric phase becomes subsumed by the bubble. The majority of realisations remain in their initial value, but a few add up to other integer values. At even higher temperature, the activity is even higher, but so is the noise on the observable.

Diffusion of winding number
Another way to quantify the activity of winding number changing processes outside the bubble is to compute the average of winding number square, where the average is over the ensemble of realisations. On one hand, if N w would evolve via diffusion in a fixed volume, one would expect this quantity to increase linearly with time as a random walk. The linear slope is the diffusion rate, well-known from the Sphaleron rate for Chern-Simons number [31]. On the other hand, integrating a fixed set of configurations over an ever increasing bubble volume, if the winding nuclei are randomly placed in space, we would again expect a linear increase as a function of volume, because bubble volume grows linearly in time with the speed of the wall. When combining the two effects, a linearly increasing volume, spreading into a region where diffusion is taking place, we expect a time dependence ofN 2 w which is faster than linear. -19 -We can try and construct a simple model of this. Let us assume that outside the bubble, there is a constant diffusion rate Γ. But as soon as a region is swallowed by the bubble, this diffusion shuts off. Also, the winding number density is assumed to be uncorrelated in space. Then we can then writē where A is the cross-section area of the box (x-y plane), and Z is the distance covered by the wall up to time t, Z = vt. Then we have the simple expression for the winding number squared inside the bubble as a function of time: proportional to the wall speed and to the diffusion rate, which is temperature-dependent, and notably quadratic in time.
In Fig. 7 we show this observable for the same set of speeds and temperatures as before. We see that for low temperatures, there is no activity at all.N 2 w is constant (or marginally increasing) as a function of time. Only for larger temperatures T /m H = 1 do we see activity in the winding number. It is indeed faster than linear in time, and it increases with wall speed. We remain wary of the results at the largest speed 0.99, but for the next-to-fastest speed v = 0.75 we also see some activity.
As a complementary measurement to establish the presence of topological transitions, we may again consider winding number squared, but rather than integrating over the entire volume inside the bubble, we only integrate over a cubic volume (z-range ∆Z as large as the x, y-ranges), immediately on the inside of the wall. The upshot is, that as the bubble passes by, the winding number freezes in. And so we record a snapshot of the winding number diffusion process in a fixed-sized volume moving through space. Modelling this again, we haveN but now the volume integral is only over A∆Z, and we find which for large enough times is proportional to Γ, independent of speed v and linear in time.
We show this in Fig. 8, where we see that at least for small enough wall speeds, the time dependence is linear, independent of v and increasing with increasing temperature. For speeds larger than 0.5, our model breaks down, and the registered activity is larger than predicted. This may be because the passage of the wall itself and the heating up from the deployed latent heat influences the diffusion rate.
-20 - Electroweak baryogenesis relies on out-of-equilibrium, CP-violating interactions between an advancing Higgs bubble wall, fermions coupled to it, and the SU(2) and U(1) gauge fields together making up the finite-temperature plasma near the walls.
Much work has gone into understanding the dynamics of this wall-plasma interaction and the transport of energy (latent heat) and fermion currents from the wall and into the plasma, where sphaleron processes generate the baryon asymmetry. This is extremely challenging, and a complete framework connecting all the different sub-processes has not been achieved.
It is appealing to attempt to simply simulate the whole thing numerically on the lattice. The dynamics of the bosonic fields is essentially classical, and the fermions are quantum mechanical but with linear equations of motion. One would hope and expect that a steady state would establish itself, where the bubble wall is driven by the thermodynamics to sweep through the plasma at a constant speed, gradually generating the baryon asymmetry. This would require that we tune the temperature to precisely the nucleation temperature; that the lattice volume is large enough for the wall to run until it reaches the steady state; that the lattice is also large enough to hold both the bubble-dynamics region and the sphaleron transition region; that the entire set of fermion modes is included, as well as the UV-completion (Hard Thermal Loops) necessary to get the bosonic UV dynamics correct.
All these techniques exist on the market, although the combined numerical effort is daunting. In this work, we have focussed on how to set up the system and carry out the simulations, with the minimal set of dynamical degrees of freedom, and trying to bypass the difficult fine-tuning of the temperature. The bubble wall is driven by an external current, and we carefully investigated how to initialize the fields, the wall, the lattice implementation and what observables to consider. We studied the extent of the out-of-equilibrium region near the bubble wall and its dependence on the wall speed. We also did a number of tests of the bubble width and shape and optimizing lattice sizes and spacings, of which the main outcome is that the physical transverse lattice size should not be smaller than used here, and the discretization of the wall not coarser. There are some transient effects of turning on the current and establishing the wall, and one may consider whether careful use of damping and energy drains could improve on this.
We have studied the (diagonal components of the) energy momentum tensor and how the energy released by the wall is transferred to the plasma both in front of the wall and inside the bubble. The fluid profile was confirmed to be self-similar and we observed, in particular, that how far outside and inside the bubble energy is projected depends strongly on the speed of the wall and the temperature of the plasma. In a simulation with fermions, it would be very interesting to similarly compute the distance that left-and righthanded currents are able to propagate before being thermalized by the plasma.
An obvious next step is making the connection to the fluid models of [1], to see whether the first-principles field theory simulations map in a sensible way to a hydrodynamical approach. Most likely, this requires a substantial scaling up of the volumes, times and statistics achieved here. This is work in progress.
We also studied the activity of winding number creation in the region near the wall, which is consistent with diffusion. The diffusion rate seems to depend somewhat on the wall speed, most likely through the heating of the plasma immediately in front of the wall. This sets the stage for including CP-violation. We did in fact implement a particular realisation of CP-violating dynamics in some test runs It is the simplest CP-violating operator involving the fields included in our model, and is a stand-in for CP-violation generated by the fermion degrees of freedom, either through the CKM-matrix [32], or for a two-Higgs-doublet model when the relative complex phase of the Higgs doublets varies through the wall (see for instance [33]). This effective term is local in space, and only sizeable near the wall (where ∂ t φ is large). In contrast, the full mechanism of electroweak baryogenesis generates a handed current near the wall, which then propagates through and back away from the wall into the plasma, where sphaleron processes take place. It is hence quite non-local and so even if the effective term may be computed in the SM and applied to the dynamics, the resulting mechanism is very different as no propagating fermion currents enter. Including CP violation in this way is quite computing heavy (a factor of 4), and we were unable to see any net asymmetry created, given the volumes and statistics available to us. We intend to return to this problem in due course. The holy grail is to include fermions dynamically. At first as a probe, in the sense of computing the fermion currents in the background of the other fields, and studying the transport of the asymmetric currents into the bulk. The implementation is underway, but scaling to a lattice volume comparable to the one used here will require massive computing resources. Eventually, fermion current backreacting on the gauge fields will have to be included as well, in order to bias the sphaleron processes.