Bubble-wall velocity in local thermal equilibrium: hydrodynamical simulations vs analytical treatment

We perform real-time hydrodynamical simulations of the growth of bubbles formed during cosmological first-order phase transitions under the assumption of local thermal equilibrium. We confirm that pure hydrodynamic backreaction can lead to steady-state expansion and that bubble-wall velocity in such case agrees very well with the analytical estimates. However, this is not the generic outcome. Instead, it is much more common to observe runaways, as the early-stage dynamics right after the nucleation allow the bubble walls to achieve supersonic velocities before the heated fluid shell in front of the bubble is formed. This effect is not captured by other methods of calculation of the bubble-wall velocity which assume stationary solutions to exist at all times and would have a crucial impact on the possible generation of both baryon asymmetry and gravitational wave signals.


Introduction
Phase transitions are present in a variety of particle-physics models.If they are first order they could create an environment for the generation of baryon asymmetry [1][2][3][4] and production of a stochastic gravitational wave background, which could be potentially observed with the next generation of detectors [5][6][7][8][9][10].Cosmological phase transitions are typically induced by quantum tunnelling or thermal fluctuation of a scalar field between non-degenerate minima of its potential and proceed via nucleation of bubbles of the energetically favourable phase.Once the bubbles are nucleated, they start to grow, driven by the potential difference between the different phases inside and outside the bubble.On the other hand, particles crossing the bubble wall exert velocity-dependent friction on it.These two forces may balance each other, leading to steady-state expansion with constant velocity.Otherwise, the bubble wall will continue to accelerate and reach a velocity very close to the speed of light 1 Correct distinction between these two situations and accurate estimation of the bubble-wall velocity is crucial for both predictions of baryon asymmetry production in theories of Electroweak Baryogenesis and modelling of gravitational wave spectra produced in the transitions.
Backreaction coming from the particles interacting with the bubble wall was estimated by solving the Boltzmann equation for all the species in the plasma together with the equation of motion for the scalar field [15][16][17][18][19][20][21][22][23][24][25][26].However, this method is quite challenging and the non-equilibrium part of the friction in the equation of motion is instead often parameterized with an effective phenomenological coefficient treated as a free parameter in hydrodynamical simulations [27,28].Recently, it has been shown that purely equilibrium hydrodynamic backreaction can inhibit the accelerating expansion [29][30][31] and a simple estimate that can be interpreted as the upper limit on bubble-wall velocity for a given model has been derived.However, this approach, similarly to methods based on the Boltzmann equation, assumes that the plasma profiles are at any time fully developed into steady-state solutions with a given wall velocity [32].The mechanism preventing the further acceleration of the wall is connected with the heating of the plasma outside of the bubble and relies on this assumption heavily.In this work, we verify this assumption using hydrodynamic simulations without non-equilibrium friction.
As a benchmark model, we use a singlet scalar extension of the Standard Model, which is a simple and well-known example of a theory in which electroweak symmetry breaking may proceed as a first-order phase transition [33][34][35][36][37][38][39][40][41][42][43][44][45][46].We use the temperature-dependent potential to find the solution corresponding to the nucleating bubble and simulate its evolution to observe as it develops into a self-similar solution.We compare the results with the bubble-wall velocities and plasma profiles obtained with analytical methods.We have found that if a steady state expanding slower than the Jouguet velocity is reached, the two methods agree with very good accuracy.Such scenarios are, however, rare and fine-tuned.In the absence of non-equilibrium friction typical bubble accelerates beyond the Jouguet velocity before heating the plasma as the self-similar profile would suggest.As a result, most bubbles develop into very fast detonations which in the absence of nonequilibrium friction continue to accelerate corresponding to the so-called runaway solution.Thus, we find that velocity estimates assuming steady-state evolution may not be valid in the majority of cases and predictions of baryon asymmetry and gravitational wave emission based on them would be incorrect.
The paper is organised as follows.In section 2 basic parameters describing cosmological phase transitions are defined, while in section 3 steady-state dynamics of the bubble expansion presenting different expansion modes are discussed.We review recent estimates of velocity [31], to which we refer further.In Section 4 we introduce the singlet scalar extension of the Standard Model and compute a high-temperature approximation of the effective potential of this theory.Then, we perform a scan of the parameter space determining transition parameters for different realisations of that model and estimate bubble-wall velocity from the matching equations.Section 5 is devoted to the derivation of equations of motion adjusted to the scalar singlet model used in real-time simulations.Finally, in section 6 we show the results of the real-time simulations and compare them with the estimations mentioned above.We summarize our results in section 7. Details of our code and numerical methods used in the simulations are reviewed in the Appendix A.

Transition parameters
Cosmological first-order phase transitions proceed via nucleation of bubbles of the broken phase from the background state of the symmetric, metastable phase.The probability of tunnelling per unit time and volume at temperature T is given with the bubble nucleation rate [47][48][49][50] where for finite temperatures the Euclidean action S = S 3 T and A(T ) = T 4 S 3 2πT 3 2 .Nucleation temperature is defined as the value such that the probability of a true vacuum bubble forming within a horizon radius grows close to unity [51], i.e.
where T c denotes the critical temperature in which both minima are equally deep.Assuming that the transition is fast, one can neglect the change of the expansion rate of the Universe (assume that Hubble parameter H(t) ≈ const).Then, eq.(2.2) reduces to which for temperatures close to the electroweak scale gives S 3 /T n ≈ 140 [6].To quantify the level of supercooling of the phase transition, we introduce the temperature ratio T n /T c .An important parameter describing the strength of the transition is the amount of latent heat released during the transition, typically normalised to the energy density of the cosmological background, resulting in [9,52] where ∆ denotes the difference between symmetric and broken phases, and V is temperaturedependent potential.For model-independent studies, one needs to define a generalized transition strength based on the trace of the energy-momentum tensor θ [53] where w s denotes enthalphy in the symmetric phase while e and p correspond to energy and pressure.Taking into account that the speed of sound might be different from the standard c s = 1/ √ 3 and assuming only weak dependence on T , one can introduce pseudotrace θ [54].We then have a more accurate prescription for the transition strength where c b is the speed of sound in the broken phase.

Late-time expansion
Late-time evolution of the bubble walls was studied in [32] using a hydrodynamical approximation in which it is assumed that the cosmic plasma can be modelled as the relativistic perfect fluid, therefore its energy-momentum tensor is given by where u µ is the four-velocity of the plasma, while p and w are respectively the pressure and the enthalpy.Conservation of energy-momentum projected onto direction u µ along the flow and the perpendicular one ūµ = γ(v, v/v) gives the following equations with ūµ u µ = 0 and ū2 = −1.As the steady-state profile has no characteristic length scale, the solution should depend only on the self-similar variable ξ = r/t, where r denotes the distance from the centre of the bubble and t is the time since nucleation.Variable ξ can be interpreted as the velocity of a given point in the profile, while the plasma at the point described by ξ moves with velocity v(ξ).Under this assumptions, equations (3.2) and (3.3) can be rewritten as and using the definition of the speed of sound in the plasma c s ≡ dp dT / de dT , they can be combined into the single hydrodynamic equation describing plasma velocity profile v(ξ) in the frame of the bubble centre with µ = ξ−v 1−ξv denoting the Lorentz-transformed fluid velocity.To proceed further, one has to define the equation of state for the plasma.The most popular choice is the so-called bag model [32,55], where the speed of sound in both phases is constant and equal c 2 s = 1/3, however, more complex options were recently considered [30,31].
Solutions of equation (3.6) with the bag-model equation of state depend only on the transition strength α and bubble-wall velocity in the stationary state ξ w .Possible types of solutions are subsonic deflagrations or supersonic detonations and hybrids schematically depicted in Fig. 1.The velocity at which the shell around the bubble disappears and the solution shifts from hybrid to detonation is given by the Chapman-Jouguet value [31]: For a detailed discussion of bag model solutions see [32].
In order to estimate bubble-wall velocity, it is necessary to estimate the total force exerted on a bubble wall from the plasma surrounding the bubble and compare it with a driving force coming from the pressure difference between the false and true vacuum.The equation of motion of a scalar field reads where the index i runs over all the species in the plasma and δf i is the non-equilibrium part of the particle distribution function as the equilibrium part is already included in the effective potential [30].Assuming a planar wall expanding with a constant velocity in z direction, the equation above can be written as and integrated leading to The left-hand side can be interpreted as the pressure difference between the true vacuum inside the bubble and the false vacuum outside.The right-hand side is the backreaction force and can be separated into the equilibrium part coming from heated plasma at the bubble front (first term) and typically subdominant [25] non-equilibrium friction (second term).Steady-state bubble-wall velocity can be determined by solving the equation of motion for the scalar field (3.8) together with the Boltzmann equations for all the particle species in the plasma.This approach is computationally demanding and typically the so-called fluid ansatz is used and solutions are found by looking for the configuration for which equation of motion has two vanishing moments (see [19,22,24,25,56] for more details).Full computations of the bubble-wall velocity including out-of-equilibrium friction are challenging, therefore to simplify the problem, it is often assumed that δf i = 0, which is called the local thermal equilibrium (LTE) scenario.
It was recently shown that assuming the LTE approximation and using the conservation of entropy it is possible to determine the bubble wall velocity in a steady state in a much simpler way.The first attempt based on the bag equation of state [30], has been already generalized to a more broad class of the equations of state and beyond the planar wall limit [31].It has been also shown, that the bubble-wall velocity can be estimated based on a few parameters evaluated at the nucleation temperature: enthalpy ratio between phases ψ N = ω b ωs , transition strength αθ and speed of sound in both phases c s and c b [31].Nevertheless, both non-equilibrium Boltzmann methods as well as simplified equilibrium approximations are based on the assumption, that stationary hydrodynamical behaviour sets in just after nucleation and do not take into account the early evolution of the system, as all the calculations are performed for stationary profiles.

Benchmark model: SM + real scalar singlet
Probably the simplest and most studied model in which a strong first-order electroweak phase transition can be realized is the real scalar singlet extension of the Standard Model [33][34][35][36][37][38][39][40][41][42][43][44][45][46], where in addition to the Standard Model Higgs doublet, the scalar sector includes the Z 2 −symmetric scalar field s.The tree-level scalar potential of this theory in a unitary gauge is given by The mass term for the Higgs µ h and the quartic coupling λ h are fixed in such a way, that in the electroweak vacuum is (h, s) = (υ, 0) with υ = 246.2GeV and the physical mass of the Higgs is m h = 125.09GeV, which leads to leaving scalar singlet mass m s , its quartic coupling λ s and the portal coupling between the Higgs and the scalar singlet λ hs as three free parameters of the model.To study the phase transition, thermal corrections to the potential need to be included: where the sum runs over all particle species of the model, and for all species n i is the number of degrees of freedom, m i (h, s) is its field-dependent mass, and J b/f is the thermal function given by where the upper (lower) sign is for bosons (fermions).Expanding the thermal function in the relativistic regime (T ≫ m) gives Therefore, the high temperature expansion (4.4) yields the well-known result with c i = 1/2 for fermions, c i = 1 for bosons and g * denoting the effective number of relativistic degrees of freedom, namely Finally, the whole effective potential of the model can be written in compact form as the tree-level potential (4.1) with temperature-dependent mass terms with where g and g ′ are electroweak couplings and y t is the Yukawa coupling for the top quark.This leads to the high-temperature approximation of the effective potential in a compact form (4.9) In this simplified treatment we neglect the presence of the Coleman-Weinberg corrections to the effective potential.To prepare a sample for our study we performed a scan over the scalar singlet mass m s and portal coupling λ hs fixing the quartic coupling to λ s = 1 with the use of the CosmoTransitions package [57].We determine the critical and nucleation temperatures as well as the transition strength for each point.
Next, using approximation derived in [31], we determine the bubble-wall velocities under the assumption of preserving local thermal equilibrium, finding a variety of deflagration and hybrid solutions, as well as runaway scenarios for the strongest transitions from the investigated sample.The results are presented in Fig. 2, where the shades of blue indicate the wall velocity while the red points mark the runaway solutions.

Real-time dynamics of the system
In hydrodynamical treatment, we describe the plasma as a perfect fluid with temperature T , characterized by the internal energy density e, pressure p and enthalpy density w.As the effective potential V eff (h, s, T ) can be interpreted as the free energy density F of the system, we can define2 p(h, s, T ) = −V eff (h, s, T ), (5.1) w(h, s, T ) = −T dV eff (h, s, T ) dT . (5. 3) The total energy-momentum tensor of the system is a sum of energy-momentum tensors for the fields3 : and the one of the fluid given by (3.1).The energy-momentum tensor of the system is conserved (∇ µ T µν = 0), however, both contributions are not conserved separately: Note that as we are interested in the local thermal equilibrium scenario, here the energy transfer between the scalar field and the plasma is possible only through the temperaturedependent effective potential, and there is no additional effective friction typically added to capture non-equilibrium effects [58].The left equality of (5.5) is satisfied if fields h and s follow standard wave equations which in spherical coordinates take the form Due to the spherical symmetry of our problem, we assume that the four-velocity of the perfect fluid is of the form u = (γ, γv, 0, 0) T with γ := (1 − v 2 ) −1/2 .We will determine the equations governing the evolution of two parameters v and p considering temporal (ν = 0) and radial component (ν = 1) of Eq. (5.5).Introducing new variables Z := wγ 2 v and τ := wγ 2 − p we get ) The system of equations for the fields (5.6) and for the plasma (5.7), together with the equations of state (5.1) and (5.2) define the dynamics of the growing bubble and was solved numerically on the lattice.For details of numerical treatment, see Appendix A.

Results of the simulations
We use our code [58] to perform real-time simulations of the bubble-wall expansion for all the points shown in Fig. 2. Each simulation was initialized with fields profiles corresponding to the nucleated critical bubble, while plasma velocity and temperature were fixed as spatially uniform, v(r) = 0 and T (r) = T n respectively.This gives a good approximation of homogeneous plasma in which nucleation takes place and allows us to carefully follow the early stages of the evolution of the bubbles and formation of the fluid profiles.This is not possible with different methods of determining bubble-wall velocities, which look for the stationary states only.Fig. 4 highlights the difference between our results in the left panel and the matching method of [31] in the right.For majority of points where the matching method predicts deflagrations and hybrids, we observe a very different behaviour.The bubble walls accelerate beyond the Jouguet velocity before the heated fluid shell in front of the bubble is formed.The corresponding evolution for a benchmark is shown in the left panel of the Fig. 3.This phenomenon is easy to understand as the bubble growth in local thermal equilibrium is very rapid.During the early stages of the evolution, the amplitude of plasma profiles is much smaller, and profiles are less sharp than the ones predicted by the steady-state solutions.This leads to a significantly lower backreaction force (see eq. (3.10)), which allows the bubble wall to accelerate further.If the profile accelerates above the Jouguet velocity before the steady state profile is reached, the stationary state will not be achieved and the bubble will continue to grow in the runaway regime.This occurs for the vast majority of the points from the scan where the matching method can be used to compute the wall velocity.In most of the parameter space where matching equations predict a deflagration or hybrid solution the simulation results in a runaway.In those cases, the simulated bubble accelerates beyond the Jouguet velocity before the heated fluid shell around the bubble is formed while analytical methods assume a steady state at all times where the shell is heated enough to cease the acceleration. .The low number of points results from the fact that the two methods agree on the final steady state solution in a very small part of the parameter space where T n /T c and αθ are tuned.The result of the simulation is much more often a runaway (see Fig. 4).
For moderately strong transitions which at the same time are not supercooled, the temperature of the plasma around the bubble front reaches more quickly the critical temperature T c .This stops further acceleration of the bubble wall, as larger wall velocity (up to Chapman-Jouguet velocity) would demand higher temperature in the peak exceeding T c which would locally tend to reverse the ongoing transition.This mechanism is called the hydrodynamical obstruction [55] and has already been demonstrated in hydrodynamical simulations [58] leading to velocity gaps within possible expansion modes.From this moment the bubble starts to asymptote to a constant velocity and the plasma profiles evolve towards stationary states predicted by the matching conditions.The evolution of the plasma velocity profile towards a stationary state compared with the analytically predicted profile for a benchmark transition is shown in the right panel of Fig. 3.While this is the case only in a very small part of the parameter space, it is worth noting that bubble-wall velocities as well as plasma profiles obtained in the simulations are in very good agreement with those found with the matching method of [31] wherever they both find non-runaway behaviour.We show a direct comparison of the predictions for this set of points in Fig. 5 finding the difference to be of the order of a few per cent at most.
Ref. [59] used a similar setup and solved the equations of motion for the very early stages of evolution (t max ≈ 0.6 GeV −1 , compared with t max ≃ 100 GeV −1 in this work.).It was shown the profiles began to form, however, the dynamical range was not large enough to verify if they reached steady-state solutions.It is also important to point out that the authors focused on a single benchmark with T n /T c ≈ 0.994 and observed the impact of hydrodynamical obstruction, which is consistent with our results.Ref. [60] also observed relaxation of the fluid to form a shell, although, in that work, the code used a fixed value for the wall velocity and the formation of a steady-state profile was only the result of this prior assumption.

Conclusions
We have studied the growth of cosmological bubbles nucleated during first-order phase transitions using real-time hydrodynamic simulations.Focusing on the bubble-wall velocity in the local thermal equilibrium, we have confirmed that pure equilibrium backreaction can lead to a steady-state solution.If this is the case and the stationary state is realised in our simulation, the final velocity agrees very well with recent predictions based on the matching equations [31].Such scenarios are however very rare and demand fine-tuning of nucleation temperature as the stationary expansion is the outcome only in the absence of any supercooling (T n /T c ≈ 1).
We have found that generically without a non-equilibrium friction bubbles expand as runaways.In the very early stages of their evolution, the nucleated bubbles accelerate quickly as the heated fluid shell around them is just beginning to form.Typically the walls accelerate past the Jouguet velocity and start to develop into detonations before the fluid around them forms into a familiar steady state profile.Once this occurs the bubble will behave as a runaway despite the fact that a solution exists for which the steady-state profile would match the vacuum pressure and cease the acceleration at a lower velocity.
We have focused on a very simple extension of the Standard Model with a neutral singlet to illustrate the prevalence of the early runaway scenario.We have found that nearly the entire parameter space of the model predicting a first-order transition would result in a runaway if one neglects the non-equilibrium friction.The only exceptions are tuned scenarios where the transition is relatively strong yet nucleation occurs extremely close to the critical temperature.
We have shown that the early stages of the evolution of bubbles in cosmological phase transitions can have a crucial impact on the resulting phenomenology.In particular, assuming a local thermal equilibrium at those early times generically leads to a runaway solution even in cases where analytical methods using fully developed steady-state profiles suggest a small terminal velocity.This would have serious consequences for the predictions of models concerning the possible production of baryon asymmetry in electroweak baryogenesis, as well as on the gravitational wave signals generated in first-order phase transitions.

A Numerical treatment
Equations (5.6) and (5.7) describing the dynamics of the model under consideration are highly nonlinear and the exact solution of the system is not known.In order to get a detailed understanding of the evolution of bubbles described by the set of equations we solve the system numerically extending the code used in our previous studies [58] to the case of two scalar fields.The extension is straight forward and to keep this manuscript self-contained we will now only briefly describe used methods.
Equations (5.6) and (5.7) have to be supplemented by the proper boundary conditions.We chose Neumann boundary conditions for both scalar fields at the centre of the bubble which correspond to r = 0.In general, the Dirichlet boundary conditions corresponding to the false vacuum expectation values should be imposed at an infinite distance r = ∞ from the nucleation site.In order to treat the problem numerically, we had to limit our computational domain to be a finite interval [0, R] with R large enough to be outside the light cone of the nucleated bubble, so imposing boundary conditions at r = R is physically equivalent to r = ∞ for a finite period of time simulated in our code.
The initial conditions for our simulations are critical bubble profiles calculated using CosmoTransitions code [57].The critical bubble which is at the boundary between collapsing and expanding bubbles is an unstable static solution of equations of motion.Due to numerical imperfections, the approximate profile may collapse instead of growing and a minimal stretching of the initial profile is necessary to guarantee the expansion of the bubble.The conserved variables Z and τ describing the state of the plasma are initialized to the values corresponding to the plasma with temperature T equal to nucleation temperature, staying at rest v = 0 with respect to the nucleation site.
We used the Galerkin discontinuous method to discretize equations (5.6) and (5.7) in space.We used the so-called mixed formulation of the equations of motion for the fields (5.6) (which are nonlinear wave equations), so the gradients of the fields and their values are discretized as independent variables.Both field values and quantities Z and τ describing plasma are interpolated using piecewise constant functions.Piecewise linear functions are used for gradients of fields.The integrals in weak forms of discretized equations are approximated by numerical quadratures which are exact for elements of the base of space of interpolation functions.The numerical fluxes were chosen such that the obtained method is the generalization of the central finite difference scheme (for planar bubble walls described in Cartesian coordinates the method would recover the central finite difference scheme of second order).analogously as described in [58].Equations of state (5.1) and (5.3) can be combined to form τ + p(h, s, T ) − 1 2 w(h, s, T ) + w(h, s, T ) 2 + 4Z 2 = 0. (A.1) For given values of τ , Z (and h, s) equation (A.1) can be numerically solved (we used Raphson-Newton method) with respect to T .The velocity of the plasma can be recovered from the definition of Z = wγ 2 v and (5.3) for a given temperature T .
We demonstrate the convergence of our numerical methods on two benchmarks from our scans in Fig 6 .The columns correspond to finer lattice spacings which we see have no impact on the results.The upper row shows a point with T n /T c ≈ 1 in which the wall reaches a finite velocity and a steady state profile forms.The lower row represents a point where the wall passes the Jouguet velocity in its early evolution and runs away.

Figure 1 .
Figure 1.Schematic representation of three different types of expanding bubbles.Colour saturation denotes the value of the plasma velocity, while black circles represent the position of the bubble wall.

Figure 2 .
Figure2.Bubble-wall velocities for the scalar singlet model in local thermal equilibrium determined using the matching method derived in[31] which assumes a steady state solution at all times.Deflagrations and hybrids are represented by blue dots with the shade depicting the wall velocity, while runaways are marked with red.No stationary detonations were found which gives rise to the gray gap in velocities.

Figure 3 .
Figure 3. Two possible scenarios for the growing bubble in local thermal equilibrium approximation: rapid expansion beyond Chapman-Jouguet velocity leading to a runaway scenario (left panel) and evolution toward a stationary state predicted by matching conditions (right panel).Solid blue curves represent the results of the real-time simulations, where darker shades correspond to later times.Red, dashed curves denote predictions of the equilibrium methods.

Figure 4 .
Figure 4. Results for the bubble wall velocity in LTE from the hydrodynamical simulations (left panel) and from analytical calculations [31] (right panel).The stationary states are shown in shades of blue corresponding to the wall velocity in both panels while runaway scenarios are shown in red.In most of the parameter space where matching equations predict a deflagration or hybrid solution the simulation results in a runaway.In those cases, the simulated bubble accelerates beyond the Jouguet velocity before the heated fluid shell around the bubble is formed while analytical methods assume a steady state at all times where the shell is heated enough to cease the acceleration.

Figure 5 .
Figure 5.Comparison of the bubble wall velocity from matching method estimates v match w [31] vs those measured in hydrodynamic simulation v sim w .The lower panel shows the difference between these two values ∆v w = v sim w − v match w