Vorticity, kinetic energy, and suppressed gravitational wave production in strong first order phase transitions

We have performed the first 3-dimensional simulations of strong first-order thermal phase transitions in the early Universe. For deflagrations, we find that the rotational component of the fluid velocity increases as the transition strength is increased. For detonations, however, the rotational velocity component remains constant and small. We also find that the efficiency with which kinetic energy is transferred to the fluid falls below theoretical expectations as we increase the transition strength. The probable origin of the kinetic energy deficit is the formation of reheated droplets of the metastable phase during the collision, slowing the bubble walls. The rate of increase in the gravitational wave energy density for deflagrations in strong transitions is suppressed compared to that predicted in earlier work. This is largely accounted for by the reduction in kinetic energy. Current modelling therefore substantially overestimates the gravitational wave signal for strong transitions with deflagrations, in the most extreme case by a factor of $10^{3}$. Detonations, however, are less affected.

The Laser Interferometer Space Antenna (LISA) is scheduled for launch in 2034, and will open up the mHz band of the emerging field of gravitational wave astronomy [1]. One of the most exciting goals of LISA is to probe the early universe through the search for gravitational wave signals from a first-order phase transition.
An important parameter of a first-order phase transition is the trace anomaly difference, which quantifies the energy available for conversion to shear stress, and hence the power of the gravitational wave signal. If the trace anomaly difference is comparable to the radiation energy density of the universe, we call the transition strong. We denote the ratio of the trace anomaly to the thermal energy α, in which case a strong transition has α ∼ 1. We call α 1 very strong; our results do not access this parameter range. Substantial progress has been made in understanding the gravitational wave production from first-order transitions with weak (α ∼ 10 −2 ) to intermediate (α ∼ 10 −1 ) strength using numerical simulations [27][28][29][30], as well as modelling [31][32][33]. While the fluid motion is welldescribed as a linear superposition of sound waves after a weak transition [27], rotational modes and turbulence are expected in stronger transitions [34,35], which could substantially affect the gravitational wave signal [36][37][38][39][40].
At the same time, investigation of the underlying particle physics models has indicated that intermediate to strong transitions are common in conservative extensions of the Standard Model [41,42], and very strong transitions are possible in models of composite Higgs and nearly conformal potentials [12][13][14][15][16][17][18]. It is also clear that LISA will be most likely to observe transitions where nonlinear effects like shocks and turbulence become important [30]. Recent work tackling the non-linear regime includes gravitational wave production from magnetohydrodynamic turbulence [43] and studies of shock collisions using a mixture of 1-dimensional simulations and modelling [44].
In this paper, we present results from the first numerical simulations of strong first-order phase transitions. We measure the fraction of the fluid kinetic energy in rotational modes, as traced by the mean-square velocity. As we increase the strength of the transition, this proportion grows substantially for deflagrations, with up to 65% of the mean square velocity found in rotational motion. The rotational proportion is far less for detonations, and remains roughly constant for all transition strengths.
We see that as the transition strength α is increased, the efficiency of fluid kinetic energy production decreases below expectation. For deflagrations, this is associated with reduced wall speeds for the expanding bubbles and reheating of the region in front of the walls, reducing the pressure difference [35,45,46]. The kinetic energy loss leads to a suppression in the gravitational wave power, by a factor which can be as small as O(10 −3 ). This means that current models substantially overestimate gravitational wave production from strong transitions with deflagrations. Detonations are less affected. We consider a system in which the scalar field undergoing the transition is modelled by a real order parameter φ, coupled to a perfect fluid. The model closely follows that used in the previous works [29,30,47], differing by a change in the effective potential and therefore the equation of state. We found that the high-temperature effec-tive potential had numerical instabilities in the presence of large temperature gradients, and we have introduced a simpler bag model equation of state, described below. The new equation of state changes only how the relevant thermodynamic parameters are realised.
Our coupled field-fluid system has a combined energymomentum tensor with v the fluid 3-velocity and γ the associated Lorentz factor. The internal energy and pressure p are and the enthalpy is w = + p.
The effective potential at zero temperature is where V c is chosen such that V 0 (φ b ) = 0, and φ b is the value of φ in the broken phase at T = 0. We denote the potential energy difference between the false and true vacua by We write the thermal effective potential of our bag model as where the function a(φ) models the change in degrees of freedom during the transition. We take the form where a 0 = (π 2 /90)g * with g * the effective number of relativistic degrees of freedom in the symmetric phase. Note that both φ = 0 and φ = φ b are stationary points of the function for all T . For our choice of a(φ) the two minima of V become degenerate at T = T c , as required.
The energy-momentum tensor can be decomposed into field and fluid parts, coupled through a friction term, Ref. [30] used a field and temperature dependent friction parameter η =ηφ 2 /T . Although this models the high temperature physics more accurately [48], strong transitions can reach small temperatures and again the hightemperature approximation breaks down. With small temperatures we also find numerical instabilities and so we have reverted to setting η to a constant. The strength of the phase transition can be parametrised by the trace anomaly difference where ∆V = V (0, T )−V (φ b , T ). The strength parameter is then where T n is the nucleation temperature and r = 3w/4 is the radiation energy density. We assume that the duration of the phase transition is much less than the Hubble time H −1 n , and so neglect the effect of expansion. This is roughly equivalent to the the statement that H n R * 1, where R * is the mean bubble separation. In this regime the contribution of bubble collisions to the gravitational wave signal is negligible.
The mean gravitational wave energy density is where V is the simulation volume, h TT ij is the transverse traceless metric perturbation and the line indicates averaging over a characteristic period of the gravitational waves. We find h TT ij in Fourier space by a standard technique [27,29,49], sourcing only with the fluid, the dominant contribution when α 1 and H n R * 1 [27,29,30]. We express the gravitational wave energy density in terms of the parameter Ω gw = ρ gw /ρ c , with ρ c the critical energy density.
We perform a series of three-dimensional simulations of the field-fluid system. The simulation code is the same as that used in Ref. [30] except for the above changes.
We scan over α for three subsonic deflagrations with asymptotic wall speeds v w = {0.24, 0.44, 0.56}, and two detonations with v w = {0.82, 0.92}. The asymptotic wall speeds, and their fluid profiles, are found with 1D version of the code [29,30,50,51], run with the same parameters to time t = 10000T −1 c . As we increase α, the maximum velocity of the asymptotic fluid profile v p increases (see Fig. 1). For each v w , there is a maximum v p , and hence a maximum strength α max , above which solutions either do not exist (subsonic deflagrations), or change into hybrids. We do not consider hybrids in this work.
The values of η needed for these wall speeds are given in the supplemental material. By comparison, the Standard Model estimate is these slices for three simulations are available at [55]. Selected stills are included in the supplemental material. We measure the RMS fluid 3-velocity v, and its irrotational and rotational parts v and v ⊥ . We also track the enthalpy-weighted RMS four-velocity U f defined as where w the mean enthalpy density. This gives an indication of the magnitude of the shear stress, the source of gravitational waves. A similar quantity U φ can be constructed to track the progress of the phase transition which is proportional to the total area of the phase boundary. We call the time at which U φ reaches its maximum the peak collision time, t pc . Note that t pc ∝ R * /v w .
To check the dependence of our key observables on lattice spacing, we perform simulations with δx T c = {2.0, 1.0, 0.5} for v w = 0.24 and v w = 0.92 and α = 0.5. We find that between δx T c = 1.0 and δx T c = 0.5 U f and Ω gw vary by O(1%). The quantity that is most sensitive to the grid is v 2 ⊥ , which increases by about 10% over the same change of δx.
In Fig. 2 we plot how U f and U φ evolve for a deflagration and a detonation, both with strength α = 0.5. We see that a rotational component of velocity v ⊥ is generated during the bubble collision phase, and that the deflagration generates v ⊥ more efficiently than the deflagration. We also see that, for the deflagration, U φ decreases more slowly than it increases, indicating a slowing down of the phase boundary.
In order to gauge the amount of kinetic energy in the rotational component of velocity, we consider the ratio of the maxima of the mean square 3-velocities v 2 ⊥,max /v 2 max . We plot this for all simulations in Fig. 3. As we increase α for the deflagrations, we see that the proportion of the velocity found in rotational modes increases dramatically, whereas for detonations it stays constant. The deflagrations with smaller wall velocities have a larger proportion of the velocity field in rotational modes. For v w = 0.24, α = 0.34 the ratio v 2 ⊥,max /v 2 max = 0.65, and if we naively extrapolate the trend in the last few simulation points up to α max this increases to 0.95. Fig. 7 of the supplemental material show that the vorticity is generated inside the bubbles, not outside where the fluid shells first interact.
To better understand transfer of energy from the scalar field to the fluid, we plot how U φ and U f change as we increase α for detonations with v w = 0.92 and deflagrations with v w = 0.44 (Fig. 4). When U φ reaches its maximum, the volumes in each phase are approximately equal. As the phase boundary sweeps out the remaining regions of metastable phase, U φ relaxes to zero. It is striking that for deflagrations the relaxation takes longer as we increase α, whereas for detonations the shape of U φ remains unchanged. This means that the phase boundaries in a deflagration must move more slowly in the later stages as the strength of the transition increases.
The reason for the slowing is that the metastable phase is reheated by the fluid shells in front of the bubble walls [35,45,46]. Towards the end of the transition the remaining metastable phase forms into hot droplets (see Fig. 7  in the supplemental material). The higher pressure in the heated droplets opposes their collapse.
For detonations, where the fluid shell develops behind the bubble wall, the shrinking regions of the metastable phase are not reheated, and the disappearance of the droplets is not hindered.
The temperature field in a slice of a simulation of a detonation is shown in Fig. 8 of the supplemental material. Fig. 4 also shows that U f increases with α, as one would expect from the increasing scalar potential energy. However, the maximum is below that expected from a single bubble, which is a good estimate of U f at low α [29,30].
To obtain the single-bubble estimate, 1D simulations of expanding spherical bubbles are performed, and the expected enthalpy-weighted RMS velocity U f,exp is that of the fluid shell around the bubble when the wall reaches a diameter of R * . We then take the ratio with maximum of U f in each simulation, shown in Fig. 5. Note that due to finite volume effects U f oscillates in our simulations, giving an O(10%) uncertainty to this estimate. For all wall speeds we see that the ratio of U f,max to U f,exp decreases as we increase the transition strength. However, for deflagrations we see that the decrease in the kinetic efficiency is more dramatic, and more rapid for slower walls: in the slowest deflagration (v w = 0.24) U f,max /U f,exp reaches 0.3. The decrease is approximately linear; a naive linear extrapolation to the maximum possible strength is indicated by open circles. The loss of kinetic energy is probably a result of the slowing of the walls discussed above, limiting the transfer of energy.
The deficit in kinetic energy can be expected to reduce the gravitational wave signal. In current modelling, the expected gravitational wave density parameter from a flow with U f,exp at time t H −1 n is whereΩ gw has been shown to be a constant of O(10 −2 ) in weak and intermediate transitions.
Here, we takeΩ gw = 10 −2 [30,56]. In Fig. 6 we plot the ratio of Ω gw /t to Ω gw,exp /t, where Ω gw /t is averaged over the final ∆t = 2R * of the simulation. In the most extreme case, v w = 0.24 and α = 0.34, the ratio is 2 × 10 −3 . This is even less than the kinetic energy suppression would suggest, a factor of (U f,max /U f,exp ) 4 8 × 10 −3 . In Fig. 1 we also color code each simulation point according to this ratio. A table of simulation parameters and measured quantities can be found in the supplemental material.
To conclude, we have performed the first 3-dimensional simulations of strong first-order phase transitions, with the strength parameter α up to an order of magnitude larger than the strongest previously studied [30].
A rotational component of velocity v ⊥ is generated during the collision phase. For deflagrations, the ratio v 2 ⊥,max /v 2 max grows rapidly with α, reaching 0.65 for v w = 0.24. For detonations, the ratio is O(10 −2 ) and shows no consistent trend with α.
For stronger phase transitions a smaller proportion of the scalar potential energy is transferred into fluid kinetic energy than is expected from the behaviour of isolated bubbles. For deflagrations, we suppose that the deficit is due to the reheating of the metastable phase slowing the bubble walls during the collision phase. The deficit can be substantial, with U f,max /U f,exp falling to ∼ 0.3 for v w = 0.24 in our simulations, and could fall as low as 0.1 using a naive linear extrapolation to the maximum possible strength at that wall speed.
The gravitational wave intensity is lower than expected, by a factor of order 10 −3 for the strongest deflagration with the lowest wall speed. This can mostly be accounted for by the kinetic energy deficit. Detonations do not suffer such a dramatic suppression, with the smallest suppression factor about 0.2 for v w = 0.92.
Our results have important consequences for gravitational waves from phase transitions. They indicate that the current model [30,56] overestimates the gravitational wave power spectrum for strong transitions, by a factor of a few for detonations, and by an order of magnitude or more for deflagrations. We plan larger simulations aiming to characterise the suppression more precisely.
The authors would like to thank Chiara Caprini, Kari Rummukainen, and Danièle Steer for helpful discussions. Our simulations made use of the resources of the Finnish Centre for Scientific Computing CSC. DC (ORCID ID 0000-0002-7395-7802) is supported by an STFC Studentship. MH (ORCID ID 0000-0002-9307-437X) acknowledges support from the Science and Technology Facilities Council, grant no. ST/P000819/1. DJW (ORCID ID 0000-0001-6986-0517) is supported by an Science and Technology Facilities Council Ernest Rutherford Fellowship, grant no. ST/R003904/1, by the Research Funds of the University of Helsinki, and by the Academy of Finland, grant no. 286769.

SUPPLEMENTAL MATERIAL
In this supplemental material we include various stills taken from movies of our simulations of strong phase transitions in the early universe, which can be seen in Fig. 7 and Fig. 8. The movies these stills have been taken from can be found in a vimeo album [55].
We also provide a table of simulation parameters and measured quantities.