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

TL;DR

This study uses 3D simulations to analyze fluid dynamics and gravitational wave production during strong first-order phase transitions in the early Universe, revealing suppressed gravitational wave signals especially in deflagrations.

Contribution

First 3D simulations of strong first-order phase transitions showing differences in fluid velocity and gravitational wave production between deflagrations and detonations.

Findings

Rotational fluid velocity increases with transition strength in deflagrations.

Kinetic energy transfer efficiency decreases as transition strength increases.

Gravitational wave energy density growth is suppressed in strong deflagrations.

Abstract

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 are less affected.

Tables

Table 1. Table 1: Key simulation parameters and measured quantities used to generate the graphs in this paper. Note that ( Ω gw H n ​ t ) ¯ ¯ subscript Ω gw subscript 𝐻 n 𝑡 \overline{\left(\dfrac{\Omega_{\text{gw}}}{H_{\text{n}}t}\right)} signifies that we average the quantity inside the brackets over the final Δ ​ t = 2 ​ R ∗ Δ 𝑡 2 subscript 𝑅 \Delta t=2R_{*} of the simulation.

$v_{\text{w}}$	$\alpha$	$T_{\text{n}}/T_{\text{c}}$	$t_{\mathrm{fin}}{T}_{\text{c}}$	$\eta /T_{\text{c}}$	${\overline{U}}_{\text{f}\text{,max}}$	${\overline{U}}_{\text{f,exp}}$	${\overline{v}}_{⟂}^2$	${\overline{v}}^2$	$\overline{\left(\frac{{\Omega }_{\text{gw}}}{H_{\text{n}}t}\right)}\frac{1}{H_{\text{n}}{R}_{\ast }}$	$\frac{{\Omega }_{\text{gw,exp}}}{H_{\text{n}}t}\frac{1}{H_{\text{n}}{R}_{\ast }}$
$0.24$	$0.050$	$0.45$	$4.8\times {10}^3$	$1.2$	$0.036$	$0.039$	$6.4\times {10}^{-5}$	$1.3\times {10}^{-3}$	$4.6\times {10}^{-8}$	$1.2\times {10}^{-7}$
$0.24$	$0.073$	$0.41$	$4.8\times {10}^3$	$1.3$	$0.048$	$0.055$	$1.5\times {10}^{-4}$	$2.3\times {10}^{-3}$	$1.4\times {10}^{-7}$	$5.0\times {10}^{-7}$
$0.24$	$0.11$	$0.37$	$4.8\times {10}^3$	$1.5$	$0.063$	$0.082$	$3.7\times {10}^{-4}$	$4.1\times {10}^{-3}$	$4.1\times {10}^{-7}$	$2.4\times {10}^{-6}$
$0.24$	$0.16$	$0.33$	$4.8\times {10}^3$	$1.8$	$0.074$	$0.12$	$9.1\times {10}^{-4}$	$5.7\times {10}^{-3}$	$6.0\times {10}^{-7}$	$9.5\times {10}^{-6}$
$0.24$	$0.23$	$0.30$	$4.8\times {10}^3$	$2.4$	$0.075$	$0.16$	$2.0\times {10}^{-3}$	$5.9\times {10}^{-3}$	$4.6\times {10}^{-7}$	$3.3\times {10}^{-5}$
$0.24$	$0.34$	$0.28$	$4.8\times {10}^3$	$5.1$	$0.066$	$0.22$	$2.9\times {10}^{-3}$	$4.4\times {10}^{-3}$	$2.3\times {10}^{-7}$	$1.2\times {10}^{-4}$
$0.44$	$0.0050$	$0.79$	$2.4\times {10}^3$	$0.34$	$0.0083$	$0.0069$	$1.7\times {10}^{-7}$	$6.9\times {10}^{-5}$	$1.0\times {10}^{-10}$	$1.2\times {10}^{-10}$
$0.44$	$0.016$	$0.59$	$2.4\times {10}^3$	$0.52$	$0.025$	$0.021$	$1.4\times {10}^{-6}$	$6.0\times {10}^{-4}$	$9.1\times {10}^{-9}$	$1.1\times {10}^{-8}$
$0.44$	$0.050$	$0.45$	$2.4\times {10}^3$	$0.66$	$0.066$	$0.061$	$1.6\times {10}^{-5}$	$4.3\times {10}^{-3}$	$5.6\times {10}^{-7}$	$7.2\times {10}^{-7}$
$0.44$	$0.073$	$0.41$	$2.4\times {10}^3$	$0.71$	$0.090$	$0.085$	$4.0\times {10}^{-5}$	$7.8\times {10}^{-3}$	$2.0\times {10}^{-6}$	$2.8\times {10}^{-6}$
$0.44$	$0.11$	$0.37$	$2.4\times {10}^3$	$0.81$	$0.12$	$0.12$	$1.1\times {10}^{-4}$	$1.3\times {10}^{-2}$	$7.0\times {10}^{-6}$	$1.0\times {10}^{-5}$
$0.44$	$0.16$	$0.33$	$2.4\times {10}^3$	$0.94$	$0.15$	$0.16$	$3.0\times {10}^{-4}$	$2.1\times {10}^{-2}$	$2.0\times {10}^{-5}$	$3.7\times {10}^{-5}$
$0.44$	$0.23$	$0.30$	$2.4\times {10}^3$	$1.2$	$0.18$	$0.22$	$1.0\times {10}^{-3}$	$3.2\times {10}^{-2}$	$4.3\times {10}^{-5}$	$1.3\times {10}^{-4}$
$0.44$	$0.34$	$0.28$	$2.4\times {10}^3$	$1.7$	$0.18$	$0.30$	$3.6\times {10}^{-3}$	$3.6\times {10}^{-2}$	$5.2\times {10}^{-5}$	$4.2\times {10}^{-4}$
$0.44$	$0.50$	$0.25$	$2.4\times {10}^3$	$3.5$	$0.19$	$0.39$	$8.0\times {10}^{-3}$	$4.0\times {10}^{-2}$	$5.8\times {10}^{-5}$	$1.2\times {10}^{-3}$
$0.56$	$0.050$	$0.45$	$2.8\times {10}^3$	$0.53$	$0.080$	$0.075$	$1.3\times {10}^{-5}$	$5.6\times {10}^{-3}$	$1.1\times {10}^{-6}$	$1.7\times {10}^{-6}$
$0.56$	$0.073$	$0.41$	$2.8\times {10}^3$	$0.59$	$0.10$	$0.10$	$3.5\times {10}^{-5}$	$9.6\times {10}^{-3}$	$3.6\times {10}^{-6}$	$5.3\times {10}^{-6}$
$0.56$	$0.11$	$0.37$	$2.8\times {10}^3$	$0.67$	$0.14$	$0.14$	$9.2\times {10}^{-5}$	$1.6\times {10}^{-2}$	$1.2\times {10}^{-5}$	$2.0\times {10}^{-5}$
$0.56$	$0.16$	$0.33$	$2.8\times {10}^3$	$0.76$	$0.18$	$0.19$	$2.4\times {10}^{-4}$	$2.7\times {10}^{-2}$	$3.6\times {10}^{-5}$	$6.4\times {10}^{-5}$
$0.56$	$0.23$	$0.30$	$2.8\times {10}^3$	$0.90$	$0.22$	$0.25$	$5.8\times {10}^{-4}$	$4.3\times {10}^{-2}$	$9.3\times {10}^{-5}$	$2.0\times {10}^{-4}$
$0.56$	$0.34$	$0.28$	$2.8\times {10}^3$	$1.2$	$0.27$	$0.33$	$1.5\times {10}^{-3}$	$6.4\times {10}^{-2}$	$2.1\times {10}^{-4}$	$6.2\times {10}^{-4}$
$0.56$	$0.50$	$0.25$	$2.8\times {10}^3$	$1.7$	$0.28$	$0.43$	$5.2\times {10}^{-3}$	$7.8\times {10}^{-2}$	$3.1\times {10}^{-4}$	$1.8\times {10}^{-3}$
$0.56$	$0.67$	$0.23$	$2.8\times {10}^3$	$2.9$	$0.30$	$0.51$	$1.1\times {10}^{-2}$	$9.0\times {10}^{-2}$	$3.0\times {10}^{-4}$	$3.7\times {10}^{-3}$
$0.82$	$0.0050$	$0.79$	$2.8\times {10}^3$	$0.11$	$0.0064$	$0.0066$	$2.3\times {10}^{-7}$	$4.0\times {10}^{-5}$	$4.8\times {10}^{-11}$	$1.0\times {10}^{-10}$
$0.82$	$0.016$	$0.59$	$2.8\times {10}^3$	$0.16$	$0.019$	$0.020$	$2.2\times {10}^{-6}$	$3.6\times {10}^{-4}$	$4.3\times {10}^{-9}$	$9.1\times {10}^{-9}$
$0.82$	$0.050$	$0.45$	$2.8\times {10}^3$	$0.18$	$0.055$	$0.061$	$1.5\times {10}^{-5}$	$2.8\times {10}^{-3}$	$2.9\times {10}^{-7}$	$7.6\times {10}^{-7}$
$0.82$	$0.073$	$0.41$	$2.8\times {10}^3$	$0.19$	$0.076$	$0.088$	$2.6\times {10}^{-5}$	$5.2\times {10}^{-3}$	$1.1\times {10}^{-6}$	$3.2\times {10}^{-6}$
$0.82$	$0.11$	$0.37$	$2.8\times {10}^3$	$0.20$	$0.11$	$0.13$	$4.7\times {10}^{-5}$	$9.4\times {10}^{-3}$	$4.1\times {10}^{-6}$	$1.4\times {10}^{-5}$
$0.82$	$0.16$	$0.33$	$2.8\times {10}^3$	$0.22$	$0.15$	$0.18$	$8.4\times {10}^{-5}$	$1.6\times {10}^{-2}$	$1.2\times {10}^{-5}$	$5.5\times {10}^{-5}$
$0.92$	$0.0050$	$0.79$	$2.4\times {10}^3$	$0.053$	$0.0051$	$0.0051$	$3.2\times {10}^{-7}$	$2.6\times {10}^{-5}$	$2.0\times {10}^{-11}$	$3.6\times {10}^{-11}$
$0.92$	$0.016$	$0.59$	$2.4\times {10}^3$	$0.086$	$0.015$	$0.016$	$3.4\times {10}^{-6}$	$2.4\times {10}^{-4}$	$1.9\times {10}^{-9}$	$3.6\times {10}^{-9}$
$0.92$	$0.050$	$0.45$	$2.4\times {10}^3$	$0.099$	$0.045$	$0.049$	$2.2\times {10}^{-5}$	$1.9\times {10}^{-3}$	$1.5\times {10}^{-7}$	$3.0\times {10}^{-7}$
$0.92$	$0.073$	$0.41$	$2.4\times {10}^3$	$0.10$	$0.064$	$0.070$	$3.6\times {10}^{-5}$	$3.7\times {10}^{-3}$	$6.0\times {10}^{-7}$	$1.3\times {10}^{-6}$
$0.92$	$0.11$	$0.37$	$2.4\times {10}^3$	$0.10$	$0.087$	$0.10$	$5.8\times {10}^{-5}$	$6.9\times {10}^{-3}$	$2.4\times {10}^{-6}$	$6.0\times {10}^{-6}$
$0.92$	$0.16$	$0.33$	$2.4\times {10}^3$	$0.11$	$0.12$	$0.14$	$8.8\times {10}^{-5}$	$1.2\times {10}^{-2}$	$8.4\times {10}^{-6}$	$2.3\times {10}^{-5}$
$0.92$	$0.23$	$0.30$	$2.4\times {10}^3$	$0.11$	$0.16$	$0.20$	$1.4\times {10}^{-4}$	$2.0\times {10}^{-2}$	$2.6\times {10}^{-5}$	$8.0\times {10}^{-5}$
$0.92$	$0.34$	$0.28$	$2.4\times {10}^3$	$0.12$	$0.21$	$0.27$	$2.5\times {10}^{-4}$	$3.2\times {10}^{-2}$	$8.1\times {10}^{-5}$	$2.9\times {10}^{-4}$
$0.92$	$0.50$	$0.25$	$2.4\times {10}^3$	$0.13$	$0.28$	$0.36$	$4.6\times {10}^{-4}$	$4.9\times {10}^{-2}$	$2.2\times {10}^{-4}$	$9.3\times {10}^{-4}$
$0.92$	$0.60$	$0.24$	$2.4\times {10}^3$	$0.14$	$0.32$	$0.42$	$6.3\times {10}^{-4}$	$5.8\times {10}^{-2}$	$3.5\times {10}^{-4}$	$1.6\times {10}^{-3}$
$0.92$	$0.67$	$0.23$	$2.4\times {10}^3$	$0.15$	$0.34$	$0.45$	$7.8\times {10}^{-4}$	$6.5\times {10}^{-2}$	$4.5\times {10}^{-4}$	$2.2\times {10}^{-3}$

Equations

$$T^{\mu \nu }={\partial }^{\mu }\varphi {\partial }^{\nu }\varphi -\frac{1}{2}{g}^{\mu \nu }(\partial \varphi {)}^2+(ϵ+p)U^{\mu }{U}^{\nu }+g^{\mu \nu }p$$

$$ϵ=3a(\varphi )T^4+V_{\text{0}}(\varphi ),p=a(\varphi )T^4-V_{\text{0}}(\varphi ),$$

$$V_{\text{0}}(\varphi )=\frac{1}{2}{M}^2{\varphi }^2-\frac{1}{3}\mu {\varphi }^3+\frac{1}{4}\lambda {\varphi }^4-V_{\text{c}}\text{,}$$

$$V(\varphi ,T)=V_{\text{0}}(\varphi )-T^4\left(a(\varphi )-a_0\right),$$

$$a(\varphi )=a_0-\frac{\Delta V_{\text{0}}}{T_{\text{c}}^4}\left[3{\left(\frac{\varphi }{{\varphi }_{\text{b}}}\right)}^2-2{\left(\frac{\varphi }{{\varphi }_{\text{b}}}\right)}^3\right],$$

$${\partial }_{\mu }{T}_{\varphi }^{\mu \nu }=-{\partial }_{\mu }{T}_{\text{f}}^{\mu \nu }=\eta U^{\mu }{\partial }_{\mu }\varphi {\partial }^{\nu }\varphi .$$

$$\Delta \theta (T)=\frac{1}{4}\frac{d}{dT}\Delta V-\Delta V\text{,}$$

$$\alpha =\Delta \theta (T_{\text{n}})/{ϵ}_{\text{r}}(T_{\text{n}}).$$

$${\rho }_{\text{gw}}=\frac{1}{32\pi G}\frac{1}{\mathcal{V}}{\int }_{\mathcal{V}}{\mathrm{d}}^{3}x\overline{{\dot{h}}_{ij}^{\text{TT}}{\dot{h}}_{ij}^{\text{TT}}}\text{,}$$

$${\overline{U}}_{\text{f}}^2=\frac{1}{\overline{w}\mathcal{V}}{\int }_{\mathcal{V}}{\mathrm{d}}^{3}xw{\gamma }^2{v}^2,$$

$${\overline{U}}_{\varphi }^2=\frac{1}{\overline{w}\mathcal{V}}{\int }_{\mathcal{V}}{\mathrm{d}}^{3}x{\partial }_i\varphi {\partial }_i\varphi ,$$

$${\Omega }_{\text{gw,exp}}=3{\tilde{\Omega }}_{\text{gw}}{\left(\frac{\overline{w}}{\overline{ϵ}}\right)}^2{\overline{U}}_{\text{f,exp}}^4(H_{\text{n}}t)(H_{\text{n}}{R}_{\ast })\text{,}$$

$$[{\partial }_{\mu }{T}^{\mu \nu }{]}_{\mathrm{field}}$$

$$[{\partial }_{\mu }{T}^{\mu \nu }{]}_{\mathrm{fluid}}$$

$${\delta }^{\nu }=\eta U^{\mu }{\partial }_{\mu }\varphi {\partial }^{\nu }\varphi \text{.}$$

$$-\ddot{\varphi }+{\nabla }^2\varphi -\frac{\partial V}{\partial \varphi }=\eta \gamma (\dot{\varphi }+v^i{\partial }_i\varphi ).$$

$$\dot{E}+{\partial }_i(E{v}^i)+p[\dot{\gamma }+\partial (\gamma v^i)]-\frac{\partial V}{\partial \varphi }\gamma (\dot{\varphi }+v^i{\partial }_i\varphi )$$

$$=\eta {\gamma }^2(\dot{\varphi }+v^i{\partial }_i\varphi {)}^2\text{.}$$

$$\dot{Z}+{\partial }_j(Z_i{v}^j)+{\partial }_{i}p+\frac{\partial V}{\partial \varphi }{\partial }_i\varphi =-\eta \gamma (\dot{\varphi }+v^j{\partial }_j\varphi ){\partial }_i\varphi \text{.}$$

$$\square h_{ij}^{TT}=16\pi G{T}_{ij}^{TT},$$

$$\square u_{ij}=16\pi G{T}_{ij}.$$

$${\tilde{h}}_{ij}^{TT}(\mathbf{k},t)={\Lambda }_{ij,lm}(\mathbf{k}){\tilde{u}}_{lm}(\mathbf{k},t)\text{,}$$

$${\Lambda }_{ij,lm}(\mathbf{k}=P_{im}(\mathbf{k})P_{jl}(\mathbf{k})-\frac{1}{2}{P}_{ij}(\mathbf{k})P_{lm}(\mathbf{k})\text{,}$$

$$P_{ij}(\mathbf{k})={\delta }_{ij}-{\hat{k}}_i{\hat{k}}_j\text{.}$$

Full text

††preprint: HIP-2019-15/TH

The Laser Interferometer Space Antenna (LISA), scheduled for launch in 2034, will open the mHz band of the emerging field of gravitational wave astronomy Audley et al. (2017). One of the most exciting goals of LISA is to probe the early universe by searching for gravitational wave signals from a first-order phase transition.

While the Standard Model is a cross-over Kajantie et al. (1996, 1997), there are many extensions with first-order phase transitions. These range from adding a scalar singlet Profumo et al. (2007); Espinosa et al. (2012); Cline and Kainulainen (2013); Profumo et al. (2015); Beniwal et al. (2019) or doublet Kakizaki et al. (2015); Dorsch et al. (2017); Basler et al. (2017), to models with spontaneously broken conformal symmetry Randall and Servant (2007); Konstandin et al. (2010); Konstandin and Servant (2011); von Harling and Servant (2018); Dillon et al. (2018); Megías et al. (2018); Bruggisser et al. (2018). There are also models with phase transitions in hidden sectors Schwaller (2015); Addazi and Marciano (2018); Aoki et al. (2017); Croon et al. (2018); Breitbach et al. (2018); Okada and Seto (2018); Hasegawa et al. (2019). Non-perturbative methods are sometimes necessary to establish the order of the phase transition Gorda et al. (2019); Gould et al. (2019); Kainulainen et al. (2019).

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 $\alpha$, in which case a strong transition has $\alpha\sim 1$. We call $\alpha\gg 1$ ‘very strong’; our results do not access this region.

Substantial progress has been made in understanding gravitational wave production from first-order transitions with weak ($\alpha\sim 10^{-2}$) to intermediate ($\alpha\sim 10^{-1}$) strength using numerical simulations Hindmarsh et al. (2014); Giblin and Mertens (2014); Hindmarsh et al. (2015, 2017), as well as modelling Hindmarsh (2018); Jinno and Takimoto (2017); Konstandin (2018). While the fluid motion is well-described as a linear superposition of sound waves after a weak transition Hindmarsh et al. (2014), rotational modes and turbulence are expected in stronger transitions Witten (1984); Kurki-Suonio (1985), which could substantially affect the gravitational wave signal Kamionkowski et al. (1994); Caprini et al. (2008); Gogoberidze et al. (2007); Caprini et al. (2009a, b); Niksa et al. (2018).

At the same time, investigation of the underlying particle physics models indicates that intermediate to strong transitions are common in conservative extensions of the Standard Model Ellis et al. (2018, 2019), and very strong transitions are possible in models of composite Higgs and nearly conformal potentials Randall and Servant (2007); Konstandin et al. (2010); Konstandin and Servant (2011); von Harling and Servant (2018); Dillon et al. (2018); Megías et al. (2018); Bruggisser et al. (2018). It is also clear that LISA will be most likely to observe transitions where nonlinear effects like shocks and turbulence become important Hindmarsh et al. (2017). Recent work tackling the non-linear regime includes gravitational wave production from magnetohydrodynamic turbulence Pol et al. (2019) and studies of shock collisions using a mixture of 1-dimensional simulations and modelling Jinno et al. (2019).

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, remaining roughly constant for all transition strengths.

As the transition strength $\alpha$ is increased, the efficiency of fluid kinetic energy production decreases below expectation. For deflagrations, this is associated with reduced wall speeds for expanding bubbles and reheating of the region in front of the walls, reducing the pressure difference Kurki-Suonio (1985); Konstandin and No (2011); Mégevand and Ramírez (2018). 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 model the phase transition with a real scalar field $\phi$, coupled to a perfect fluid. We assume that there is no extra physics generating a magnetic field either before or during the phase transition. The model follows that used in previous work Ignatius et al. (1994); Hindmarsh et al. (2015, 2017), differing by a change in the effective potential and therefore the equation of state. Our previous work used the high-temperature expansion of the one-loop thermal effective potential, and we found that in stronger transitions, the total energy could drop below the scalar potential energy, which is unphysical. In this scenario, our algorithm would compute the temperature to be imaginary, causing a crash. Indeed, the high-temperature expansion is known to fail well below $T_{\text{c}}$; for example, the speed of sound diverges and then becomes imaginary. To fix this we have introduced a simpler bag model equation of state, described below. The new equation of state changes only how the relevant thermodynamic parameters $\alpha$ and $v_{\text{w}}$ are realised in terms of the parameters of the potential and field-fluid coupling term. The flows around the expanding bubbles, and hence the gravitational wave spectrum, depend on the underlying theory only through $\alpha$ and the wall speed $v_{\text{w}}$, with the overall frequency scale set by the redshifted mean bubble separation.

Our coupled field-fluid system has energy-momentum tensor

[TABLE]

where $U=\gamma(1,{\textbf{v}})$, with fluid 3-velocity v and associated Lorentz factor $\gamma$. The internal energy $\epsilon$ and pressure $p$ are

[TABLE]

and the enthalpy is $w=\epsilon+p$.

The zero-temperature effective potential is

[TABLE]

where $V_{\text{c}}$ is chosen such that $V_{\text{0}}(\phi_{\text{b}})=0$, and $\phi_{\text{b}}$ is the value of $\phi$ in the broken phase at $T=0$. We denote the potential energy difference between the vacua by $\Delta V_{\text{0}}=V_{\text{0}}(0)-V_{\text{0}}(\phi_{\text{b}})$.

We write the thermal effective potential of our bag model as

[TABLE]

where $a(\phi)$ models the change in degrees of freedom during the transition. We take

[TABLE]

where $a_{0}=(\pi^{2}/90)g_{*}$ with $g_{*}$ the effective number of relativistic degrees of freedom in the symmetric phase. Both $\phi=0$ and $\phi=\phi_{\text{b}}$ are stationary points of the function for all $T$. For our choice of $a(\phi)$ the minima of $V$ become degenerate at $T=T_{\text{c}}$, as required.

The energy-momentum tensor can be decomposed into field and fluid parts, coupled through a friction term,

[TABLE]

Ref. Hindmarsh et al. (2017) used a field- and temperature-dependent friction parameter $\eta=\tilde{\eta}{\phi^{2}}/{T}$. Although this models high temperature physics more accurately Liu et al. (1992), strong transitions can reach small temperatures and again the high-temperature approximation fails. With small temperatures we also find numerical instabilities and so revert to using a constant $\eta$.

The phase transition strength is parametrised by the trace anomaly difference

[TABLE]

where $\Delta V=V(0,T)-V(\phi_{\text{b}},T)$. The strength parameter is then

[TABLE]

where $T_{\text{n}}$ is the nucleation temperature and $\epsilon_{\text{r}}=3w/4$ the radiation energy density.

We assume that the duration of the phase transition is much less than the Hubble time $H_{\text{n}}^{-1}$, and neglect the effect of expansion. This is comparable to the statement that $H_{\text{n}}R_{*}\ll 1$, where $R_{*}$ is the mean bubble separation. In this regime the contribution of bubble collisions to the gravitational wave signal is negligible. To neglect expansion the final simulation time $t_{\mathrm{fin}}$ must also be much smaller than $H_{\text{n}}^{-1}$. For all our simulations $t_{\mathrm{fin}}\leq 10\,R_{*}$.

The mean gravitational wave energy density is

[TABLE]

where ${\mathcal{V}}$ is the simulation volume, $h^{\text{TT}}_{ij}$ is the transverse traceless metric perturbation and the line indicates averaging over a characteristic period of the gravitational waves. We find $h^{\text{TT}}_{ij}$ in Fourier space by a standard technique Garcia-Bellido et al. (2008); Hindmarsh et al. (2014, 2015), sourced only by the fluid, the dominant contribution when $\alpha\lesssim 1$ and $H_{\text{n}}R_{*}\ll 1$ Hindmarsh et al. (2014, 2015, 2017).

We express the gravitational wave energy density in terms of the parameter $\Omega_{\text{gw}}={\rho_{\text{gw}}}/{\rho_{\text{c}}}$, with $\rho_{\text{c}}$ the critical energy density. Our assumptions on $\alpha$ and $H_{\text{n}}R_{*}$ ensure that $\Omega_{\text{gw}}\ll 1$ at all times. They also ensure that the gravitational backreaction is negligible compared to the pressure forces, as the wavelength of the density perturbations $\sim R_{*}$ is much less than the Jeans length $\sim c_{\text{s}}/H_{\text{n}}\sqrt{\delta}$, where $\delta$ is the energy density contrast averaged over the wavelength being considered111For the parameter space we consider, we determined that the variation of the energy density is at most a factor of twenty in the asymptotic fluid profile that develops around an expanding bubble; $\delta$ will be less than this due to averaging over a given wavelength.. We leave a deeper analysis of gravitational backreaction in the case $H_{\text{n}}R_{*}\sim c_{\text{s}}/\sqrt{\delta}$ to a later study.

We perform a series of three-dimensional simulations of the field-fluid system. The simulation code is the same as used in Ref. Hindmarsh et al. (2017) except for the above changes.

We scan over $\alpha$ for three subsonic deflagrations with asymptotic wall speeds $v_{\text{w}}=\{0.24,\,0.44,\,0.56\}$, and two detonations with $v_{\text{w}}=\{0.82,\,0.92\}$. The asymptotic wall speeds, and their fluid profiles, are found with a spherically symmetric version of the code Hindmarsh et al. (2015, 2017); Kurki-Suonio and Laine (1996, 1995), run with the same parameters until $t=10000T_{\text{c}}^{-1}$. As we increase $\alpha$, the maximum velocity of the asymptotic fluid profile $v_{\mathrm{p}}$ increases. For each $v_{\text{w}}$, there is a maximum $v_{\mathrm{p}}$, and hence a maximum strength $\alpha_{\mathrm{max}}$, above which solutions either do not exist (subsonic deflagrations), or change into hybrids. We do not consider hybrids here.

The values of $\eta$ needed for these wall speeds are given in the supplemental material. By comparison, the Standard Model estimate is $\eta\simeq 3\,\phi_{\text{b}}^{2}/T_{\text{c}}$  John and Schmidt (2001); Liu et al. (1992); Moore and Prokopec (1995).

All simulations have the number of bubbles $N_{\text{b}}=8$, lattice spacing $\delta x=1.0\,T_{\text{c}}^{-1}$, timestep $\delta t=0.2\,T_{\text{c}}^{-1}$, and $L^{3}=960^{3}$ lattice sites, giving a mean bubble separation $R_{*}=L\,\delta x/N_{\text{b}}^{1/3}=480\,T_{\text{c}}^{-1}$. All bubbles are nucleated simultaneously with a gaussian profile at the same locations at the start of each simulation. The initial profile of the bubbles is insignificant as they approach the same asymptotic profile.

We fix $g_{*}=106.75$, $M^{2}=0.0427\,T_{\text{c}}^{2}$, $\mu=0.168\,T_{\text{c}}$ and $\lambda=0.0732$, in turn fixing $\phi_{\text{b}}=2.0\,T_{\text{c}}$. This sets the relative change in degrees of freedom to $[a(\phi_{\text{b}})-a_{0}]/a_{0}=5.9\times 10^{-3}$. To change the transition strength we vary $T_{\text{n}}$.

We output slices of the temperature $T$, fluid speed $v$ and vorticity magnitude $\left|\nabla\times\textbf{v}\right|$. Movies created from these slices are available at Cutting (2019). Selected stills are included in the supplemental material.

We measure the RMS fluid 3-velocity $\overline{v}$, and its irrotational and rotational parts $\overline{v}_{\parallel}$ and $\overline{v}_{\perp}$. We also track the enthalpy-weighted RMS four-velocity $\overline{U}_{\text{f}}$ defined as

[TABLE]

where $\overline{w}$ the mean enthalpy density. This gives an indication of the magnitude of the shear stress, the source of gravitational waves.

A similar quantity $\overline{U}_{\phi}$ can be constructed to track the progress of the phase transition

[TABLE]

proportional to the total area of the phase boundary. We call the time when $\overline{U}_{\phi}$ reaches its maximum the peak collision time, $t_{\text{pc}}$. Note that $t_{\text{pc}}\propto R_{*}/v_{\text{w}}$. To see how these global quantities evolve during a detonation and a deflagration see Fig. 10 in the supplemental material.

To check the dependence of our key observables on lattice spacing, we perform simulations with the same physical volume and various lattice spacings $\delta x\,T_{\text{c}}=\{2.0,1.5,1.0,0.75,0.5\}$ for $v_{\text{w}}=0.24$ and $v_{\text{w}}=0.92$ and $\alpha=0.5$. We find that $\overline{v}_{\perp\text{,max}}^{2}$, $\overline{U}_{\text{f\text{,max}}}$, and $\Omega_{\text{gw}}$ converge with lattice spacing. We perform a quadratic fit with $\delta x$ for each quantity, finding that $\overline{U}_{\text{f\text{,max}}}$ and $\Omega_{\text{gw}}$ differ from the continuum limit by $\mathrm{O}(1\%)$. The quantity that is most sensitive to the grid is $\overline{v}_{\perp\text{,max}}^{2}$ which we underestimate from the continuum limit by $\simeq 25\%$. We also test convergence of key observables with timestep, finding in all cases that convergence is better for $\delta t$ than for $\delta x$. It is important to check how close the colliding bubbles are to their asymptotic profile. We find that spherically symmetric bubbles with diameter $R_{\text{c}}$ have at worst $\overline{U}_{\text{f}}$ within $20\%$ of the asymptotic $\overline{U}_{\text{f}}$. In the supplemental material we show our convergence test results and the deviation from asymptotic $\overline{U}_{\text{f}}$ for all $v_{\text{w}}$ and $\alpha$ considered here.

From our simulations we see that a rotational component of velocity is generated during the bubble collision phase. In order to gauge the relative amount of kinetic energy in the rotational component of velocity, we consider the ratio of the maxima of mean square 3-velocities $\overline{v}_{\perp\text{,max}}^{2}/\overline{v}_{\text{max}}^{2}$. We plot this in Fig. 1. As we increase $\alpha$ 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 in rotational modes. For $v_{\text{w}}=0.24$, $\alpha=0.34$ the ratio $\overline{v}_{\perp\text{,max}}^{2}/\overline{v}_{\text{max}}^{2}=0.65$, and if we naively extrapolate the trend in the last few points up to $\alpha_{\mathrm{max}}$ this increases to $0.95$.

Fig. 11 of the supplemental material shows 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 $\overline{U}_{\phi}$ and $\overline{U}_{\text{f}}$ change as we increase $\alpha$ for detonations with $v_{\text{w}}=0.92$ and deflagrations with $v_{\text{w}}=0.44$ (Fig. 2). When $\overline{U}_{\phi}$ reaches its maximum, the volumes in each phase are approximately equal. As the phase boundary sweeps out the remaining regions of metastable phase, $\overline{U}_{\phi}$ relaxes to zero. It is striking that for deflagrations the relaxation takes longer as we increase $\alpha$, whereas for detonations the shape of $\overline{U}_{\phi}$ remains unchanged. The phase boundaries in a deflagration must therefore move more slowly in the later stages, as the transition strength increases.

The reason for the slowing is that the metastable phase is reheated by the fluid shells in front of the bubble walls Kurki-Suonio (1985); Konstandin and No (2011); Mégevand and Ramírez (2018). Towards the end of the transition the remaining metastable phase forms into hot droplets (see Fig. 11 in the supplemental material). The higher pressure inside the droplets opposes their collapse.

For detonations, where the fluid shell develops behind the bubble wall, shrinking regions of the metastable phase are not reheated (see Fig. 12 in the supplemental material).

Fig. 2 also shows that $\overline{U}_{\text{f}}$ increases with $\alpha$, as one expects from the increasing scalar potential energy. However, the maximum is below that expected from a single bubble, which is a good estimate of $\overline{U}_{\text{f}}$ at low $\alpha$ Hindmarsh et al. (2015, 2017).

To obtain the single-bubble estimate, simulations of expanding spherical bubbles are performed, and the expected enthalpy-weighted RMS velocity $\overline{U}_{\text{f,exp}}$ is that of the fluid shell when the wall reaches a diameter of $R_{*}$. We then take the ratio with the maximum of $\overline{U}_{\text{f}}$ in each simulation, shown in Fig. 3. Note that due to finite volume effects $\overline{U}_{\text{f}}$ oscillates in our simulations, giving an O(10%) uncertainty to this estimate.

For all wall speeds, the ratio of $\overline{U}_{\text{f\text{,max}}}$ to $\overline{U}_{\text{f,exp}}$ decreases as we increase the transition strength. However, for deflagrations the decrease in the kinetic efficiency is more dramatic, and more rapid for slower walls: in the slowest deflagration ($v_{\text{w}}=0.24$), $\overline{U}_{\text{f\text{,max}}}/\overline{U}_{\text{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 discussed above, limiting the transfer of energy.

The deficit in kinetic energy can be expected to reduce the gravitational wave signal. In current modelling Caprini et al. (2016); Hindmarsh et al. (2017), the expected gravitational wave density parameter from a flow with $\overline{U}_{\text{f,exp}}$ at time $t\ll H_{\text{n}}^{-1}$ is

[TABLE]

where $\tilde{\Omega}_{\text{gw}}$ has been shown to be a constant of $\mathrm{O}(10^{-2})$ in weak and intermediate transitions. Here, we take $\tilde{\Omega}_{\text{gw}}=10^{-2}$. In Fig. 4 we plot the ratio of $\Omega_{\text{gw}}/t$ to $\Omega_{\text{gw,exp}}/t$, where $\Omega_{\text{gw}}/t$ is averaged over the final $\Delta t=2R_{*}$ of the simulation. In the most extreme case, $v_{\text{w}}=0.24$ and $\alpha=0.34$, the ratio is $2\times 10^{-3}$. This is even less than the kinetic energy suppression suggests, a factor of $(\overline{U}_{\text{f\text{,max}}}/\overline{U}_{\text{f,exp}})^{4}\simeq 8\times 10^{-3}$.

A table of simulation parameters and measured quantities can be found in the supplemental material.

We have performed the first 3-dimensional simulations of strong first-order phase transitions, with the strength parameter $\alpha$ up to an order of magnitude larger than those previously studied Hindmarsh et al. (2017).

A rotational component of velocity $\overline{v}_{\perp}$ is generated during the collision phase. For deflagrations, the ratio $\overline{v}_{\perp\text{,max}}^{2}/\overline{v}_{\text{max}}^{2}$ grows rapidly with $\alpha$, reaching $0.65$ for $v_{\text{w}}=0.24$. For detonations, the ratio is O($10^{-2}$)—showing no consistent trend with $\alpha$.

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 reheating of the metastable phase slowing the bubble walls during the collision phase. The deficit can be substantial, with $\overline{U}_{\text{f\text{,max}}}/\overline{U}_{\text{f,exp}}$ falling to $\sim 0.3$ for $v_{\text{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_{\text{w}}=0.92$.

Our results have important consequences for gravitational waves from phase transitions. They indicate that the current model Caprini et al. (2016); Hindmarsh et al. (2017) 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 estimate222We use PTPlot v1.01, http://www.ptplot.org/ptplot/Caprini:2019egz, to compute and plot signal to noise ratio (SNR) curves from first-order phase transitions for LISA. We choose $g_{*}=106.75$ and an optimistic $T_{\text{n}}=100\,\mathrm{GeV}$. The resulting plot shows that an SNR of 10 requires $\overline{U}_{\text{f}}$ of at least 0.07 for all $H_{\text{n}}R_{*}$. that to obtain a signal to noise ratio of 10 $\overline{U}_{\text{f}}$ must be at least 0.07. Therefore the kinetic energy suppression we observe will probably render transitions with $v_{\text{w}}=0.24$ unobservable except for within a very small region of parameter space. Faster walls suffer less suppression, though the observable parameter space is still reduced.

We plan larger simulations to characterise more precisely the suppression, and its effect on observability.

I Supplemental Material

I.1 Field and fluid equations of motion

In order to obtain the equations of motion for our field and fluid system we focus on the coupling between the field and fluid parts of our energy momentum tensor. The current of the energy-momentum tensor can be split into field and fluid parts and coupled through a dissipative friction term,

[TABLE]

We can then write this coupling term as

[TABLE]

From these two equations we can extract the equation of motion for our system. By taking Eq (13) and dividing through by $\delta^{\nu}\phi$ we obtain

[TABLE]

We find the equation of motion for the fluid energy density $E=\gamma\epsilon$ by contracting Eq (14) with $U_{\nu}$ giving

[TABLE]

Finally we obtain an expression for the fluid momentum density $Z_{i}=\gamma(\epsilon+p)U_{i}$ by considering the spatial components of Eq (14),

[TABLE]

I.2 Gravitational waves

To obtain the gravitational wave energy density we must first calculate the transverse traceless perturbations in the metric, $h^{TT}_{ij}$. We operate in linearised gravity and therefore the equation of motion for $h^{TT}_{ij}$ is

[TABLE]

where $T^{TT}_{ij}$ is the transverse traceless projection of the energy-momentum tensor.

Due to the numerical cost of computing the transverse traceless components of the energy-momentum tensor, it is useful to instead track an auxiliary tensor $u_{ij}$ Garcia-Bellido et al. (2008) which evolves according to

[TABLE]

Then to obtain $h^{TT}_{ij}$ from $u_{ij}$ we apply the transverse traceless projector in wave space,

[TABLE]

where

[TABLE]

and

[TABLE]

This method then allows us to only need to perform the necessary Fourier transforms and projections to calculate the gravitational wave energy density at regular intervals rather than every timestep.

I.3 Resolution convergence

To ensure the validity of our simulations we performed a series of lattice resolution checks. To do this we repeated two simulations with $v_{\text{w}}=0.44$ and $v_{\text{w}}=0.92$ and $\alpha=0.5$ for a variety of different lattice spacings and timesteps while keeping the total physical volume and duration of the simulations fixed. We plot the convergence of several key quantities with $\delta x$ in Fig. 5 through Fig. 7. We also plot a quadratic fit for the convergence of each quantity with $\delta x$. We can see that all quantities converge. The worst convergence is for $\overline{v}_{\perp\text{,max}}^{2}$ which for $\delta x=1.0\,T_{\text{c}}^{-1}$ we underestimate by up to $25\%$ from the extrapolation to the continuum limit. We also performed tests for convergence of our simulations with $\delta t$. For $\delta t=0.2T_{\text{c}}^{-1}$ the error from our simulations is within $\sim 1\%$ from the continuum limit for $\overline{\left(\Omega_{\text{gw}}/H_{\text{n}}t\right)}(1/H_{\text{n}}R_{\text{c}})$ and $\overline{U}_{\text{f\text{,max}}}$ and $\sim 5\%$ for $\overline{v}_{\perp\text{,max}}^{2}$.

I.4 Convergence to asymptotic fluid flow

In addition to testing convergence with lattice spacing, we also check how close the fluid shells around colliding bubbles in our simulation are to the final asymptotic profiles. To do this we perform spherically symmetric 1D simulations of isolated bubbles and calculate $\overline{U}_{\text{f,exp}}$ from the fluid shell when the bubble has diameter $R_{\text{c}}$. We then compare this to $\overline{U}_{\text{f,exp}}$ calculated from the fluid shell at $t=10000T^{-1}_{\mathrm{c}}$, i.e when the diameter is $\gg R_{\text{c}}$ and the profile has reached its asymptotic solution. We plot the ratio of these two quantities for all $\alpha$ and $v_{\text{w}}$ in Fig. 8. We can see that the bubbles colliding with the diameter of the average bubble seperation are within $20\%$ of the asymptotic $\overline{U}_{\text{f}}$ for all simulations. We believe this to be sufficient for this study, and save a further investigation on the convergence with increasing $R_{\text{c}}$ for a future work.

I.5 Parameter space

In order to understand the regions of parameter space mapped out by our simulations, it can be illuminating to plot the asymptotic maximum fluid flow velocity $v_{\mathrm{p}}$ against the wall velocity $v_{\text{w}}$ for each simulation point. We do this in Fig. 9. Plotting the parameter space in this manner separates subsonic-deflagrations, supersonic-deflagrations, and detonations. Stronger phase transitions with the same wall velocity have a larger value of $v_{\mathrm{p}}$. Transitions with $v_{\mathrm{p}}>v_{\text{w}}$ are forbidden as this would mean that in the wall frame fluid was flowing out from the bubble. We additionally colour each simulation point by the suppression factor in gravitational waves found in our study.

I.6 Evolution of global quantities

In Fig. 10 we plot how $\overline{U}_{\text{f}}$ and $\overline{U}_{\phi}$ evolve for a deflagration and a detonation, both with strength $\alpha=0.5$. We see that a rotational component of velocity $\overline{v}_{\perp}$ is generated during the bubble collision phase, and that the deflagration generates $\overline{v}_{\perp}$ more efficiently than the detonation. We also see that, for the deflagration, $\overline{U}_{\phi}$ decreases more slowly than it increases, indicating a slowing down of the phase boundary.

I.7 Simulation slice stills

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. 11 and Fig. 12. The movies these stills have been taken from can be found in a Vimeo album Cutting (2019).

I.8 Simulation parameters and measurements