Shape of the acoustic gravitational wave power spectrum from a first order phase transition

TL;DR

This study uses large-scale simulations to analyze the shape of gravitational wave spectra from early universe phase transitions, confirming some theoretical predictions and identifying new features related to sound shell thickness.

Contribution

It provides the first detailed numerical comparison of the acoustic gravitational wave spectrum with the sound shell model, revealing new spectral features and parameter dependencies.

Findings

The $k^{-3}$ behavior appears clearly for detonations.

Deflagrations show steeper high-$k$ decay and additional features.

A simple broken power law fit can describe the spectrum for certain wall speeds.

Abstract

We present results from large-scale numerical simulations of a first order thermal phase transition in the early universe, in order to explore the shape of the acoustic gravitational wave and the velocity power spectra. We compare the results with the predictions of the recently proposed sound shell model. For the gravitational wave power spectrum, we find that the predicted $k^{-3}$ behaviour, where $k$ is the wavenumber, emerges clearly for detonations. The power spectra from deflagrations show similar features, but exhibit a steeper high-$k$ decay and an extra feature not accounted for in the model. There are two independent length scales: the mean bubble separation and the thickness of the sound shell around the expanding bubble of the low temperature phase. It is the sound shell thickness which sets the position of the peak of the power spectrum. The low wavenumber behaviour of the velocity power spectrum is consistent with a causal $k^{3}$, except for the thinnest sound shell, where it is steeper. We present parameters for a simple broken power law fit to the gravitational wave power spectrum for wall speeds well away from the speed of sound where this form can be usefully applied. We examine the prospects for the detection, showing that a LISA-like mission has the sensitivity to detect a gravitational wave signal from sound waves with an RMS fluid velocity of about $0.05c$, produced from bubbles with a mean separation of about $10^{-2}$ of the Hubble radius. The shape of the gravitational wave power spectrum depends on the bubble wall speed, and it may be possible to estimate the wall speed, and constrain other phase transition parameters, with an accurate measurement of a stochastic gravitational wave background.

Tables

Table 1. Table 1: Input parameters and derived equilibrium and non-equilibrium quantities for our simulations. Our parameters are the effective number of relativistic degrees of freedom g ∗ subscript 𝑔 g_{*} , scalar potential parameters [see Eq. ( 6 )] and nucleation temperature T n subscript 𝑇 n T_{\text{n}} . From these we obtain the latent heat ℒ ℒ {\cal L} , phase boundary tension σ 𝜎 \sigma and the thickness ℓ ℓ \ell . For studying phase transitions, it is useful to also compute the equilibrium value of the scalar field at the nucleation temperature ϕ b ​ ( T n ) subscript italic-ϕ b subscript 𝑇 n \phi_{\text{b}}(T_{\text{n}}) , transition strength parameters α w subscript 𝛼 𝑤 \alpha_{w} and α θ subscript 𝛼 𝜃 \alpha_{\theta} [see Eqs. ( 8 ) and ( 10 )], the Chapman-Jouguet speed for detonations v CJ subscript 𝑣 CJ v_{\text{CJ}} , and the critical bubble radii R c subscript 𝑅 c R_{\text{c}} . Finally, we use Eq. ( 20 ) to compute a value for the Hubble constant H n subscript 𝐻 n H_{\text{n}} .

Parameter	Weak	Intermediate
$g_{\ast }$	$106.75$	$106.75$
$T_0/T_{\text{c}}$	$1/\sqrt{2}$	$1/\sqrt{2}$
$\gamma$	$2/9$	$4/9$
$A$	$0.1990$	$0.1990$
$\lambda$	$0.0792$	$0.0396$
$T_{\text{n}}/T_{\text{c}}$	$0.86$	$0.80$
$ℒ/T_c^4$	$0.7013$	$5.6102$
$\sigma /T_c^3$	$0.1558$	$0.8816$
$\ell T_{\text{c}}$	$3$	$2.1213$
${\varphi }_{\text{b}}\left(T_{\text{n}}\right)/T_{\text{c}}$	$1.7838$	$3.5810$
${\alpha }_w\left(T_{\text{n}}\right)$	$0.010$	$0.084$
${\alpha }_{\theta }\left(T_{\text{n}}\right)$	$0.0046$	$0.050$
$v_{\text{CJ}}$	$0.63$	$0.73$
$R_{\text{c}}{T}_{\text{c}}$	$8.1$	$4.3$
$H_{\text{n}}/\sqrt{G}{T}_{\text{c}}^2$	$12.686$	$10.978$

Table 2. Table 2: Simulation parameters η ~ ~ 𝜂 \tilde{\eta} (field-fluid coupling), N b subscript 𝑁 b N_{\text{b}} (number of bubbles nucleated) and lattice spacing d ​ x 𝑑 𝑥 dx , with the corresponding bubble wall speed v w subscript 𝑣 w v_{\text{w}} , mean bubble separation R ∗ subscript 𝑅 R_{*} and sound shell width Δ ​ R ∗ = R ∗ ​ | v w − c s | / c s Δ subscript 𝑅 subscript 𝑅 subscript 𝑣 w subscript 𝑐 s subscript 𝑐 s \Delta R_{*}=R_{*}|v_{\text{w}}-c_{\text{s}}|/c_{\text{s}} . The potential parameters and derived quantities for the “weak” and “intermediate” phase transitions are given in Table  1 .

Type	$\tilde{\eta }$	$N_{\text{b}}$	$dx{T}_{\text{c}}$	$v_{\text{w}}$	$R_{\ast }{T}_{\text{c}}$	$\Delta R_{\ast }{T}_{\text{c}}$
Weak	0.19	84	2	$0.92$	1918	714.4
	0.35	84	2	$0.80$	1918	533.8
	0.51	84	2	$0.68$	1918	289.5
	0.59	11	1	$0.56$	1889	58.51
	0.93	84	2	$0.44$	1918	598.7
	0.51	5376	2	$0.68$	480	72.38
	0.59	5376	2	$0.56$	480	14.86
	0.93	5376	2	$0.44$	480	149.7
Int.	0.17	84	2	0.92	1918	714.4
	0.40	11	1	0.72	1889	374.2
	0.62	84	2	0.44	1918	598.7

Table 3. Table 3: Wall speed v w subscript 𝑣 w v_{\text{w}} , average bubble separation R ∗ subscript 𝑅 R_{*} , with peak bubble collision time t pc subscript 𝑡 pc t_{\text{pc}} , the maximum fluid RMS velocity U ¯ f max superscript subscript ¯ 𝑈 f max \overline{U}_{\text{f}}^{\text{max}} , the maximum contribution of transverse fluid motion U ¯ f , ⟂ max superscript subscript ¯ 𝑈 f perpendicular-to max \overline{U}_{\text{f},\perp}^{\text{max}} , and two estimates for U ¯ f subscript ¯ 𝑈 f \overline{U}_{\text{f}} based on 1D fluid profiles and a fitting formula given in Ref.  Espinosa  et al. ( 2010 ) .

Type	$v_{\text{w}}$	$R_{\ast }{T}_{\text{c}}$	$t_{\text{pc}}{T}_{\text{c}}$	${10}^3{\overline{U}}_{\text{f}}^{\text{max}}$	${10}^3{\overline{U}}_{\text{f},⟂}^{\text{max}}$	${10}^3{\overline{U}}_{\text{f}}^{1D}$	${10}^3{\overline{U}}_{\text{f}}^{\mathrm{Esp}}$
Weak	0.92	1918	1210	4.60	0.0833	5.31	5.32
	0.80	1918	1380	5.75	0.0665	6.39	6.50
	0.68	1918	1630	8.65	0.116	9.17	10.0
	0.56	1889	1860	13.8	0.190	14.3	14.7
	0.44	1918	2520	7.51	0.0775	7.70	7.76
	0.68	480	430	8.74	0.252	9.17	10.0
	0.56	480	480	11.7	0.498	14.3	14.7
	0.44	480	660	6.99	0.131	7.70	7.76
Int.	0.92	1918	1180	43.7	0.869	51.6	53.7
	0.72	1889	1480	65.0	1.97	72.8	95.0
	0.44	1918	2650	54.5	2.77	51.7	67.7

Table 4. Table 4: Wall speed v w subscript 𝑣 w v_{\text{w}} and mean bubble separation R ∗ subscript 𝑅 R_{*} , with the resulting integral scale of the fluid flow ξ f end superscript subscript 𝜉 f end \xi_{\text{f}}^{\text{end}} , the integral scale of the gravitational wave power ξ gw end superscript subscript 𝜉 gw end \xi_{\text{gw}}^{\text{end}} , and the dimensionless gravitational wave amplitude parameters Ω ~ gw ξ f superscript subscript ~ Ω gw subscript 𝜉 f \tilde{\Omega}_{\text{gw}}^{\xi_{\text{f}}} and Ω ~ gw R ∗ superscript subscript ~ Ω gw subscript 𝑅 \tilde{\Omega}_{\text{gw}}^{R_{*}} .

Type	$v_{\text{w}}$	$R_{\ast }{T}_{\text{c}}$	${\xi }_{\text{f}}^{\text{end}}{T}_{\text{c}}$	${\xi }_{\text{gw}}^{\text{end}}{T}_{\text{c}}$	${10}^2{\tilde{\Omega }}_{\text{gw}}^{{\xi }_{\text{f}}}$	${10}^2{\tilde{\Omega }}_{\text{gw}}^{R_{\ast }}$
Weak	0.92	1918	1490	1620	1.5	1.2
	0.80	1918	1290	1600	2.2	1.4
	0.68	1918	888	1410	1.5	0.62
	0.56	1889	530	865	1.4	0.32
	0.44	1918	1450	1750	1.4	1.1
	0.68	480	268	323	1.5	0.88
	0.56	480	233	416	1.6	0.86
	0.44	480	416	493	1.5	1.3
Int.	0.92	1918	1530	1780	2.6	2.0
	0.72	1889	1100	1180	3.3	1.8
	0.44	1918	1980	2090	1.6	1.7

Table 5. Table 5: The fit parameters of the ansatz ( 35 ). ( k ​ R ∗ ) max subscript 𝑘 subscript 𝑅 max (kR_{*})_{\rm max} is the location of the maximum of the power spectrum, and A 𝐴 A its amplitude, along with the estimate from Eq. ( 39 ).

Type	$v_{\text{w}}$	$R_{\ast }$	${\left(k{R}_{\ast }\right)}_{\max }$	$A$	$A_{\text{est}}$
Weak	0.92	1918	8.6	$1.4\times {10}^{-11}$	$1.9\times {10}^{-11}$
	0.80	1918	10.4	$3.1\times {10}^{-11}$	$5.9\times {10}^{-11}$
	0.68	1918	18.3	$8.1\times {10}^{-11}$	$14\times {10}^{-11}$
	0.44	1918	9.9	$8.2\times {10}^{-11}$	$12\times {10}^{-11}$
Int.	0.92	1918	8.5	$1.6\times {10}^{-7}$	$2.7\times {10}^{-7}$
	0.72	1889	16.1	$3.7\times {10}^{-7}$	$13\times {10}^{-7}$
	0.44	1918	6.9	$4.3\times {10}^{-7}$	$5.3\times {10}^{-7}$

Equations

$${\tau }_{ij}^{\varphi }={\partial }_i\varphi {\partial }_j\varphi ,{\tau }_{ij}^{\text{f}}=W^{2}w{V}_i{V}_j,$$

$${\overline{U}}_{\text{f}}^2=\frac{1}{\bar{w}\mathcal{V}}{\int }_{\mathcal{V}}{d}^{3}x{\tau }_{ii}^{\text{f}},$$

$${\overline{U}}_{\varphi }^2=\frac{1}{\bar{w}\mathcal{V}}{\int }_{\mathcal{V}}{d}^{3}x{\tau }_{ii}^{\varphi }.$$

$$p(T,\varphi )$$

$$ϵ(T,\varphi )$$

$$V(\varphi ,T)=\frac{1}{2}\gamma (T^2-T_0^2){\varphi }^2-\frac{1}{3}AT{\varphi }^3+\frac{1}{4}\lambda {\varphi }^4.$$

$$\mathcal{L}(T)=w(T,0)-w(T,{\varphi }_{\text{b}})$$

$${\alpha }_w=\frac{\mathcal{L}(T)}{{ϵ}_{\text{r}}(T)}.$$

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

$${\alpha }_{\theta }=\frac{\Delta \theta (T)}{{ϵ}_{\text{r}}(T)}.$$

$$\begin{gathered} {\lambda }_{ij,kl}(\text{k})⟨{\tau }_{\text{f}}^{ij}(\text{k},t_1){\tau }_{\text{f}}^{kl}({\text{k}}^{\prime },t_2)⟩ \\ ={\Pi }^2(k,t_1,t_2)(2\pi {)}^3\delta (\text{k}+{\text{k}}^{\prime }), \end{gathered}$$

$${\lambda }_{ij,kl}(\text{k})=P_{ik}(\text{k})P_{jl}(\text{k})-\frac{1}{2}{P}_{ij}(\text{k})P_{kl}(\text{k})$$

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

$${\hat{q}}^i{\hat{q}}^j⟨{\tilde{v}}_{{\text{q}}_1}^i(t_1){\tilde{v}}_{{\text{q}}_2}^{\ast j}(t_2)⟩=G(q,t_1,t_2)(2\pi {)}^3\delta ({\text{q}}_1-{\text{q}}_2).$$

$$\frac{d{V}^2}{d\ln (q)}=\frac{q^3}{2{\pi }^2}\left(G(q,t,t)+G^{\perp }(q,t,t)\right).$$

$$\square u_{ij}=(16\pi G){\tau }_{ij}^{\text{f}}.$$

$${\rho }_{\text{gw}}=\frac{1}{32\pi G}⟨{\dot{h}}_{ij}{\dot{h}}_{ij}⟩,$$

$$\frac{d{\Omega }_{\text{gw}}}{d\ln (k)}=\frac{1}{12H^2}\frac{k^3}{2{\pi }^2}{P}_{\dot{h}}(\text{k},t),$$

$$⟨{\dot{h}}_{ij}(\text{k},t){\dot{h}}_{ij}({\text{k}}^{\prime },t)⟩=P_{\dot{h}}(\text{k},t)(2\pi {)}^3\delta (\text{k}+{\text{k}}^{\prime }).$$

$$H^2=\frac{8\pi G}{3}ϵ(0,T).$$

$$\begin{gathered} \frac{d{\Omega }_{\text{gw}}(k)}{d\ln (k)} \\ =3{\Gamma }^2{\overline{U}}_{\text{f}}^4(H_{\text{n}}t)(H_{\text{n}}{L}_{\text{f}})\frac{(k{L}_{\text{f}}{)}^3}{2{\pi }^2}{\tilde{P}}_{\text{gw}}(k{L}_{\text{f}}), \end{gathered}$$

$${\tau }_{\text{sh}}\sim L_{\text{f}}/{\overline{U}}_{\text{f}}.$$

$${\Omega }_{\text{gw}}^{\text{ac}}=3{\Gamma }^2{\overline{U}}_{\text{f}}^4(H_{\text{n}}{L}_{\text{f}}){\tilde{\Omega }}_{\text{gw}},$$

$${\tilde{\Omega }}_{\text{gw}}=\frac{1}{2{\pi }^2}{\int }_0^{\infty }dx{x}^2{\tilde{P}}_{\text{gw}}(x)$$

$${\Omega }_{\text{gw}}^{\text{tu}}\sim {\Gamma }^2{\overline{U}}_{\text{f}}^4(L_{\text{f}}/{\tau }_{\text{sh}}),$$

$${\Omega }_{\text{gw}}^{\text{ac}}\sim {\Omega }_{\text{gw}}^{\text{tu}}\sim {\Gamma }^2{\overline{U}}_{\text{f}}^3(L_{\text{f}}{H}_{\text{n}}{)}^2.$$

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

$$R_{\ast }=(\mathcal{V}/N_{\text{b}}{)}^{\frac{1}{3}}.$$

$$dx\ll \ell \ll \Delta R_{\ast }\lesssim R_{\ast }\ll L_{\text{f}}.$$

$${\left({\overline{U}}_{\text{f}}^{1D}\right)}^2=\frac{3}{v_{\text{w}}^3}\int d\xi {\xi }^2{W}^2{v}^2,$$

$${\kappa }_{\text{v}}{\alpha }_{\theta }=\frac{1}{\overline{ϵ}\mathcal{V}}\int d^{3}x{\tau }_{\text{f}}^{ii}.$$

Full text

††preprint: HIP-2017-02/TH

I Introduction

The first direct observation of gravitational waves by the Laser Interferometer Gravitational Wave Observatory (LIGO) in 2015 has opened a new and unexplored window to the cosmos Abbott et al. (2016a, b). Even more excitingly, while the original detection was related to an astrophysical process, the merger of two black holes, gravitational waves will also allow us to directly probe processes in the very early universe, such as inflation Bartolo et al. (2016), topological defects Sanidas et al. (2012, 2013); Aasi et al. (2014); Blanco-Pillado et al. (2014) and first order phase transitions Caprini et al. (2016). To study gravitational wave signals from these and other sources the Laser Interferometer Space Antenna (LISA) is due to launch about a decade from now Audley et al. (2017). The LISA Pathfinder mission has recently demonstrated the technological feasibility of such a mission Armano et al. (2016).

LISA will have particular sensitivity in the millihertz frequency range, making it an ideal instrument to observe gravitational wave signals from phase transitions in the electroweak era, corresponding to roughly 10 picoseconds after the big bang. First order phase transitions proceed via the nucleation, expansion, collision and merger of bubbles of the low temperature (broken) phase. They can source gravitational waves in a number of ways. Firstly, gravitational radiation is produced by the collisions of the bubble walls, where the scalar order parameter changes from the symmetric to the broken phase. If the phase transition occurred in vacuum, this would be the only source of gravitational waves. The energy momentum tensor of this source is well approximated by the envelope of a configuration of infinitely thin shells Kosowsky et al. (1992a, b); Kamionkowski et al. (1994); this is known as the envelope approximation Huber and Konstandin (2008). It leads to a characteristic high-wavenumber (UV) fall-off of the spectrum proportional to $k^{-1}$, where $k$ is the wave number, which has been recently confirmed by both numerical simulations  Weir (2016) and analytic modelling Jinno and Takimoto (2017).

If the scalar bubbles expand in a hot plasma, as expected to be the case in the early universe, friction from the plasma will slow down the walls, which after initial acceleration will expand with a constant speed $v_{w}$. In this case most energy released by the transition will be transferred from the scalar field into the plasma. Only a tiny fraction, on the order of microphysics scale to bubble radius at collision, will remain in the scalar field. So the scalar field source is completely negligible for a thermal cosmological phase transition111According to Ref. Bodeker and Moore (2009) for very strong thermal transitions, friction from the plasma will not stop the bubbles from accelerating towards the speed of light (“runaway bubbles”). However, it is expected that additional friction from higher order corrections related to particle production will modify this result Bodeker and Moore (2017). So runaway bubbles are expected to correspond to very fast standard detonations.. The energy transferred to the medium can either go to heat or fluid motion. Numerical simulations show that the energy momentum tensor of the fluid after bubble collisions corresponds to an ensemble of sound waves. These sound waves in turn are an efficient source of gravitational radiation Hogan (1986); Hindmarsh et al. (2014, 2015). An analytical model of the velocity perturbations in this acoustic phase of the transition based on a picture of the acoustic phase fluctuations as overlapping shells of sound waves has recently been proposed Hindmarsh (2016), in which the UV fall-off of the spectrum is roughly $k^{-3}$, so distinctly different from the transition in the vacuum case. This means that an observation of such a gravitational wave signal with LISA could allow to distinguish between the two cases, giving valuable information on the nature of the transition. For very strong transitions it is expected that the acoustic phase turns over into a turbulent stage Kosowsky et al. (2002); Nicolis (2004); Caprini and Durrer (2006); Gogoberidze et al. (2007); Caprini et al. (2009a); Kahniashvili et al. (2010, 2012); Kisslinger and Kahniashvili (2015). This turbulence will continue to produce gravitational radiation until it decays.

How do the different possible components of gravitational radiation from a thermal first order phase transition compare? Collisions of the bubble walls contribute on the order of $R_{\text{c}}/R_{*}$ which for electroweak bubbles is about $10^{-14}$ and therefore completely negligible.

The relative weight of the contributions from the acoustic phase and the subsequent turbulent phase crucially depend on the strength of the transition. For very strong transitions, the plasma will quickly enter the turbulent phase, which will then have a noticeable if not dominating impact on the resulting gravitational wave signal.

Sufficiently weak transitions will not become turbulent before gravitational wave production is effectively switched off by cosmic expansion. The spectrum will then be dominated by gravitational radiation from sound waves. This will be discussed in more detail in Section II.3.

Magnetic fields may be present at the electroweak phase transition, generated either earlier in the history of the Universe Durrer and Neronov (2013), or by charge separation at the bubble walls if CP violation is present Baym et al. (1996). In this case, any turbulent flow will redistribute energy between the fluid and the magnetic field towards equipartition Kahniashvili et al. (2013). Magnetohydrodynamic turbulence is expected to have its own characteristic gravitational wave signal Caprini and Durrer (2006); Kahniashvili et al. (2010), although there is still significant uncertainty about the shape of the power spectrum.

The peak frequencies of these different comments are typically quite similar. For example, according to Ref. Caprini et al. (2016) the peak frequency of gravitational waves obtained from magnetohydrodynamics (MHD) is about a factor of 1.5 higher than that from sound waves, but given our limited understanding of turbulence this may change as a result of future research. Also this difference may well depend on the parameters of the phase transition, such as the wall velocity.

In this paper we present results from extensive numerical simulations, building on earlier work reported in Refs. Hindmarsh et al. (2014, 2015). We model the system by a scalar order parameter field coupled to a relativistic fluid by means of a phenomenological friction term. The resulting Klein-Gordon equation coupled to relativistic hydrodynamics is solved on a lattice. We use it to study the acoustic phase at an unprecedented level of accuracy.

We show detailed velocity and gravitational wave power spectra, for both deflagrations and detonations, and compare them to the predictions of the sound shell model. The UV power laws agree with the model in the case of detonations, and the prediction of an intermediate $k^{1}$ power law for wall speeds close to the speed of sound is also corroborated. We establish that there are two length scales in the power spectrum: the mean bubble separation, and the width of the sound pulse around the expanding bubble wall, the “shell” of the sound shell model. The UV power law for deflagrations is steeper than the model prediction, with an interesting break or knee.

We show that the gravitational wave power spectrum for wall speeds well away from the speed of sound can be modelled with a broken power law and an amplitude proportional to the fourth power of the RMS fluid velocity and to the ratio of the fluid flow length scale to the Hubble length. We use the model to forecast the sensitivity of LISA Audley et al. (2017) to acoustically generated gravitational waves.

Characterising the fluid flow length scale by the mean bubble separation, we find that the peak sensitivity is to transitions with a mean bubble separation of order $10^{-2}$ of the Hubble length at a transition with critical temperature $10^{2}\;\textrm{GeV}$. Transitions generating an RMS fluid velocity of about $0.05$ (in natural units) give rise to acoustic gravitational waves with signal-to-noise ratio of about 10.

We also estimate the timescale on which the acoustic waves become shocked due to non-linear evolution, which would cause the velocity and gravitational wave power spectra to deviate from their acoustic form. A significant part of the parameter space generating observable gravitational waves is likely to feature shocks and eventually turbulence, for which further simulation is required to establish an accurate power spectrum.

When interfaced with a particular microphysics realisation of a phase transition to provide the main input parameters, latent heat, bubble size and wall velocity, our results allow accurate estimations of the resulting gravitational wave signal. Conversely, the differences in the shapes we observe point the way towards estimating the wall speed and constraining combinations of other phase transition parameters from accurate observations of a primordial gravitational wave power spectrum.

In the following section we recap the physics of the acoustic generation of gravitational waves after a first-order thermal phase transition; in Section III we discuss our numerical methods, highlighting aspects of our approach which differ from Refs. Hindmarsh et al. (2014, 2015); our results for the fluid velocity power spectrum can be found in Section IV and for gravitational waves in Section V. We then compare these results to the power law ansatz used for the LISA Cosmology Working Group report (Ref. Caprini et al. (2016)) in Section VI. Our conclusions are in Section VII.

II Acoustic gravitational waves

The source of the gravitational waves is shear stress in the system, induced by the nucleation, explosive growth and merger of bubbles of the Higgs phase. These perturbations take the form of compression and rarefaction waves laid down around the growing bubbles - that is, the sound of the Higgs explosions.

II.1 Thermodynamics

The sources of shear stress are the order parameter $\phi$ and the relativistic fluid to which it is coupled. Because we need only the transverse-traceless part of the energy-momentum tensor, it is sufficient to consider as a source tensor $\tau_{ij}=\tau^{\phi}_{ij}+\tau^{\text{f}}_{ij}$, which is decomposed into fluid and field pieces according to

[TABLE]

where $w=\epsilon+p$ is the enthalpy density, $\epsilon$ is the energy density, $p$ is the pressure, $V_{i}$ is the fluid 3-velocity, and $W$ is the corresponding Lorentz factor. Unless the transition is strongly supercooled, most of the available energy of the transition goes into thermal and kinetic energy of the fluid; the scalar contribution is negligible.

It is useful to describe the overall amplitude of the fluid shear stress by a root mean square (RMS) four-velocity $\overline{U}_{\text{f}}$ defined through

[TABLE]

where ${\mathcal{V}}$ is the averaging volume and $\bar{w}$ is the volume averaged enthalpy density. One can define a similar quantity for the scalar field

[TABLE]

Although these are not quite the magnitudes of the transverse traceless part of the shear stress, they are easy to compute, and do have a direct connection to the gravitational wave amplitude for random fields, as we shall see in Eq. (21).

The fluid energy density and pressure have a contribution from the scalar order parameter of the phase transition $\phi$, through its effective potential $V(\phi,T)$,

[TABLE]

where $g_{*}$ is the effective number of relativistic degrees of freedom.

Following Enqvist et al. (1992); Ignatius et al. (1994), we use a simple quartic form for the potential:

[TABLE]

The detailed form is not important: its function is to supply a metastable state with a latent heat

[TABLE]

where $\phi_{\text{b}}$ is the equilibrium value of the field in the symmetry-broken phase at temperature $T$. The strength of the transition can be parametrised by the ratio of the latent heat to the total radiation density in the high temperature symmetric phase222Note that the radiation density has no unique definition in a system with a scalar order parameter in a thermal bath: here we follow Ref. Espinosa et al. (2010) and define it as $\epsilon_{\text{r}}(T)=3w(T,0)/4$.

[TABLE]

A commonly used alternative is the difference in the trace anomaly divided by a conventional factor of 4, or

[TABLE]

where $\Delta V(T)=V(0,T)-V(\phi_{\text{b}},T)$. Expressed relative to the energy density in the symmetric phase

[TABLE]

Other important parameters are the surface tension and width of the phase boundary $\sigma$ and $\ell$, which can be computed straightforwardly from the parameters of the potential Enqvist et al. (1992); Hindmarsh et al. (2015). The transition takes place at the nucleation temperature $T_{\text{n}}$, at which the radius of the critical bubble is denoted $R_{\text{c}}$.

II.2 Shear stress and velocity correlations

As mentioned above, the dominant source of shear stress is the fluid, unless $\alpha\gg 1$ by the measures of the phase transition strength outlined above, and the scalar field is so weakly coupled with the fluid that the walls continue to accelerate until collision. Our simulations explore the more generic situation, \alpha\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}1.

One can characterise the fluid source by the unequal time correlator (UETC) of the shear stress $\Pi^{2}$ Caprini et al. (2009b); Figueroa et al. (2013), defined by projecting out the spatially transverse and traceless part of the energy-momentum tensor

[TABLE]

where

[TABLE]

and

[TABLE]

The shear stresses are the result of the sound waves, which can be characterised by the longitudinal part of the velocity unequal time correlator $G(q,t_{1},t_{2})$. This is defined from the Fourier transform of the velocity field $\tilde{v}^{i}_{\textbf{q}_{1}}$ through

[TABLE]

The transverse part $G^{\perp}$ can be defined analogously. We will be interested in the velocity power spectrum,

[TABLE]

In our simulations, the longitudinal part is always much greater than the transverse part (see Table 3), reflecting the dominance of sound waves in the fluid perturbations.

II.3 Gravitational waves

The transverse traceless metric perturbation $h_{ij}$ is extracted by projection from an auxiliary tensor $u_{ij}$ satisfying the equation

[TABLE]

The energy density in gravitational waves is

[TABLE]

where $h_{ij}$ is the transverse-traceless projection of $u_{ij}$. It is most useful to consider the gravitational wave power spectrum relative to the critical density, defined as

[TABLE]

where $H$ is the Hubble parameter and $P_{\dot{h}}$ is the power spectral density, defined from the two-point correlation

[TABLE]

In our simulations, the expansion of the universe is scaled out using the scale invariance of the relativistic fluid equations Brandenburg et al. (1996); Hindmarsh et al. (2015). Nonetheless, one can still formally compute the Hubble rate from the Friedmann equation

[TABLE]

It was shown in Ref. Hindmarsh et al. (2015) that, after the acoustic source has been on for time $t$, the dimensionless gravitational wave power spectrum takes the form

[TABLE]

where $\Gamma=1+\bar{p}/\bar{\epsilon}$ is the adiabatic index, $\overline{U}_{\text{f}}$ is the RMS fluid velocity, $H_{\text{n}}$ is the Hubble rate at the bubble nucleation temperature $T_{\text{n}}$, $L_{\text{f}}$ is the characteristic length of the fluid flow, and $\tilde{P}_{\text{gw}}$ is a dimensionless spectral density for the gravitational waves.

It was also shown that, provided that turbulence does not develop within a Hubble time, the effective time for which the acoustic source operates is precisely the Hubble time, so that $H_{\text{n}}t\to 1$.

The turbulence timescale, both for appearance and decay, is the shock appearance or eddy turn-over time Landau and Lifshitz (1987); Pen and Turok (2016)

[TABLE]

The maximum duration of all our simulations is much less that $\tau_{\text{sh}}$, and so no turbulence develops. Our results therefore apply to flows for which $\overline{U}_{\text{f}}\ll L_{\text{f}}H_{\text{n}}$.

With this assumption, the total gravitational wave energy density from the acoustic phase is

[TABLE]

where

[TABLE]

is a dimensionless parameter, quantifying the efficiency with which shear stress is converted to gravitational waves. A significant result in Hindmarsh et al. (2015) was that this parameter is approximately independent of the length scale and RMS velocity of the fluid flow.

The formula (21) and (23) are derived from general considerations of the velocity correlation function, so one can also estimate the relative amplitude of gravitational waves produced during the turbulent phase. Eq. (23) arises from a flow with RMS velocity $\overline{U}_{\text{f}}$ with characteristic length scale $L_{\text{f}}$ starting when the Hubble parameter is $H_{\text{n}}$. The turbulent phase starts when the Hubble parameter is of order $\tau_{\text{sh}}^{-1}$, with similar RMS velocity and characteristic length scale. Hence once can estimate that the density parameter of gravitational waves produced during the turbulent phase is

[TABLE]

a factor $(H_{*}\tau_{\text{sh}})^{-1}$ smaller than from the acoustic phase. Therefore, for transitions with $\overline{U}_{\text{f}}\ll L_{\text{f}}H_{\text{n}}$, the gravitational wave signal from the turbulent phase can be neglected. One reaches the same conclusion with a more careful derivation based on Eq. (A9) of Ref. Hindmarsh et al. (2015).

In the case of a strong phase transition, turbulence develops in less than a Hubble time, the lifetime of both the acoustic and turbulent phases is of order $\tau_{\text{sh}}$. Hence, using $t\sim\tau_{\text{sh}}\sim L_{\text{f}}/\overline{U}_{\text{f}}$ in (21) and integrating over wavenumber, one finds

[TABLE]

Hence for strong phase transitions the acoustic and turbulent signals should have similar magnitudes.

III Methods

The system is a set of coupled partial differential equations governing the evolution of the scalar field $\phi$ and the relativistic ideal fluid with 4-velocity $U^{\mu}$. We use the techniques previously described in Hindmarsh et al. (2015) (see also Refs. Kurki-Suonio and Laine (1995, 1996), and the textbooks Wilson and Matthews (2003); Rezzolla and Zanotti (2013)).

The field and fluid parts of the system are coupled together through a dissipative term that turns field stress-energy $T^{\mu\nu}_{\phi}$ into fluid stress-energy $T^{\mu\nu}_{\text{f}}$ such that,

[TABLE]

where $\eta$ is in general a function of $\phi$ and $T$ with mass dimension 1. In previous work, including our own Hindmarsh et al. (2014, 2015), this was taken to be a constant. In this work we take $\eta=\tilde{\eta}{\phi^{2}}/{T}$, where $\tilde{\eta}$ is a dimensionless parameter, which is better motivated by the underlying physics Moore and Prokopec (1995); Huber and Sopena (2013). The fluid velocity around the expanding bubbles for a given wall speed and phase transition strength is minimally affected by the change, as it is determined purely by hydrodynamics, except right at the bubble wall Steinhardt (1982); Kurki-Suonio and Laine (1996); Espinosa et al. (2010).

The parameters chosen in the numerical simulations are given in Table 1. All dimensionful quantities are expressed in terms of the critical temperature $T_{\text{c}}$, defined from $\Delta V(T_{\text{c}})=0$.

We take the effective number of relativistic degrees of freedom to be the Standard Model’s high temperature value, although its exact value is not important333From Eq. (8) and the potential Eq. (6), one can see that the phase transition strength $\alpha_{w}$ can be kept constant if $g_{*}$, $\gamma$, $A$ and $\lambda$ are all scaled by some constant $C$. This will change the correlation length $\ell$ and hence the bubble wall thickness by a factor $\sqrt{C}$ but will otherwise have little impact on the position or amplitude of the resulting gravitational wave power spectrum, so long as Eq. (29) holds.. The value of Newton’s constant $G$ is arbitrary, as we will compute quantities which are independent of $G$.

We simulate at two different transition strengths, which we label “weak” ($\alpha\sim 10^{-2}$) and “intermediate” ($\alpha\sim 10^{-1}$). These have the same phase transition strengths as our “weak” and “intermediate” simulations in Refs. Hindmarsh et al. (2014, 2015), but the different form of the field-fluid coupling term, and the changed value of the relativistic degrees of freedom $g_{*}$ mean that the simulations are not identical. The differences, however, are minor and do not affect the results.

For both weak and intermediate categories, we simulate at a variety of wall velocities giving both deflagrations and detonations, including some which move at close to the Chapman-Jouguet speed $v_{\text{CJ}}$, defined as the wall speed at which the exit velocity of the fluid in the wall frame is the speed of sound Landau and Lifshitz (1987). For a deflagration, $v_{\text{CJ}}=c_{\text{s}}$, while for a detonation it depends on the strength of the transition, but is always greater than $c_{\text{s}}$ (see e.g. Ref. Espinosa et al. (2010)). For a weak transition, $v_{\text{CJ}}\simeq c_{\text{s}}(1+\sqrt{2\alpha_{\theta}})$. Note that the values of the transition strength parameter in Table 1 do not take into account the small increase in the temperature of the fluid near the moving bubble wall, and hence the stated Chapman-Jouguet speed is slightly higher than the true minimum speed of a detonation.

Each bubble is nucleated by inserting a scalar field configuration with Gaussian profile, as described in Ref. Hindmarsh et al. (2015). This profile is slightly larger than the critical bubble radius $R_{\text{c}}$. The bubble expands and perturbs the fluid, which evolves towards the scaling solution (c.f. Fig. 2).

Our bubbles are nucleated simultaneously, rather than with a physical nucleation rate. It is possible to rescale the results to yield the power spectrum from a more realistic nucleation rate. This was demonstrated in Ref. Weir (2016) and is discussed further in Section VI. We simulate with a variety of bubble numbers $N_{\text{b}}$, which controls the average bubble separation

[TABLE]

After time $t$, each bubble is surrounded by a sound shell Hindmarsh (2016) of approximate thickness $\Delta R=\Delta v_{\text{w}}t$, where $\Delta v_{\text{w}}=|v_{\text{w}}-c_{\text{s}}|$.

For near-Jouguet bubble walls the sound shell is typically quite thin relative to the bubble separation. In the sound shell model Hindmarsh (2016), the velocity power spectrum peaks at $k\Delta R_{*}\sim 1$, where $\Delta R_{*}$ is the average bubble shell thickness at collision.

The bubble wall speeds and bubble separations are chosen to explore the dependence of the power spectra on $R_{*}$ and $\Delta R_{*}$, and also to compensate for the limited dynamic range of the simulations. The values of $R_{*}$ and $\Delta R_{*}$ are listed in Table 2.

In principle we need sufficiently large lattices to explore

[TABLE]

In practice, the bubble wall does not need high resolution, as the detailed dynamics of the scalar field are not important beyond the transfer of energy to the fluid. For the near-Jouguet transitions, we found it useful to explore the wavenumber range $kR_{*}\ll 1$ and $k\Delta R_{*}\gg 1$ separately by adjusting the number of bubbles. Simulations with large $R_{*}$ have the advantage that the velocity field is closer to its asymptotic self-similar form.

Velocity and gravitational wave power spectra are computed on cubic periodic lattices with $N=4200$ points per side. The lattice spacing $dx$, bubble number $N_{\text{b}}$, and bubble wall speed $v_{\text{w}}$, along with the field-fluid coupling required to obtain this speed, are listed in Table 2.

We generally run with $dx=2/T_{\text{c}}$ which has been established to work well for single bubble self-collisions for a weak deflagration and $R_{*}\approx 384/T_{\text{c}}$ Hindmarsh et al. (2015).

Wall velocities $v_{\text{w}}$ depend on discretisation effects, and it is difficult to determine the final asymptotic value from numerical simulations of limited duration. For concreteness, our quoted $v_{\text{w}}$ values are determined from spherically symmetric simulations of a single bubble with $dx\,T_{\text{c}}=0.2$ and are measured at time $5000/T_{\text{c}}$. For coarser lattice spacings the actual wall velocity will be slightly smaller, but the difference in the cases studied here is at most 3% (this effect is more pronounced for faster wall velocities Cutting (2017)). In our simulations, the bubbles collide before the asymptotic profile is reached (see Fig. 2), and this is generally more significant than discretisation effects, particularly for simulations with smaller mean bubble separations.

For simulations with $v_{\text{w}}$ close to $c_{\text{s}}$ the lattice discretisation effects on the fluid profile were more significant than for other cases. In the weak, $v_{\text{w}}=0.59$, $R_{*}\approx 1900/T_{\text{c}}$ case, we ran simulations at both $dx=1/T_{\text{c}}$ and $dx=2/T_{\text{c}}$ as a check against lattice effects in our final results. Agreement in the velocity power spectra is excellent until $kR_{*}\approx 100$, deteriorating to an error of about $50\%$ at $kR_{*}\approx 200$. The discrepancy is more pronounced in the gravitational wave power spectrum because it convolves the power at different wavelengths Hindmarsh (2016).

IV Results: fluid velocity

In Fig. 1 we plot $\overline{U}_{\text{f}}$ and $\overline{U}_{\phi}$ against time, showing the development, completion, and aftermath of the phase transition. We divide the transition into three phases Hindmarsh et al. (2015): the expansion phase before any bubble collisions take place, the collision phase, and the acoustic phase. These can be traced in the Figures, using the fact that $\overline{U}_{\phi}^{2}$ is proportional to the surface area of the phase boundaries. In the expansion phase, $\overline{U}_{\phi}$ grows linearly with time. In the collision phase, $\overline{U}_{\phi}$ falls below the initial linear behaviour, peaks, and drops to zero, marking the start of the acoustic phase. The RMS velocity $\overline{U}_{\text{f}}$ also grows with $\overline{U}_{\phi}$ in the expansion and collision phase, and levels off in the acoustic phase. We take the peak of $\overline{U}_{\phi}$ to mark the collision time of the bubbles $t_{\text{pc}}$.

The RMS fluid velocity is generally constant in the acoustic phase, as noted in Ref. Hindmarsh et al. (2015), although some reduction can be seen in the transitions with bubbles expanding at near the Chapman-Jouguet speed, where the fluid velocity profile is narrower and peaks at higher values (see deflagrations with $v_{\text{w}}=0.56$ and detonations with $v_{\text{w}}\simeq 0.7$ in Fig. 2). It is possible that this reduction represents the beginning of turbulent transport of energy to the lattice dissipation scale; the time for shocks and turbulence to appear is in these cases $\tau_{\text{sh}}\sim 10^{5}/T_{\text{c}}$, only a factor of order 10 longer than the simulation time.

Fig. 2 also shows how the fluid profiles evolve as the bubbles expand. The profile takes some time to settle to its asymptotic self-similar form, and we plot both this, and the profile at the peak collision time $t_{\text{pc}}$. The difference is particularly noticeable for small bubble separations where collisions happen much earlier (see Fig. 2b); the fluid profiles are smooth and do not have the characteristic sharp edges at the wall position. This affects the high wavenumber behaviour of the velocity power spectra, as we shall see.

In Table 3 we list the maximum RMS fluid velocity $\overline{U}_{\text{f}}^{\text{max}}$, along with the transverse component $\overline{U}_{\text{f},\perp}^{\text{max}}$. We also give theoretical values for the mean square fluid velocity estimated in two ways. Firstly, $\overline{U}_{\text{f}}^{\mathrm{1D}}$ is obtained by integrating the numerical 1D fluid profiles out to $t=7000/T_{\text{c}}$, according to (see Refs. Espinosa et al. (2010); Hindmarsh et al. (2015))

[TABLE]

where $\xi=r/t$ is the scaled radius, $v(\xi)$ is the radial fluid velocity and $W(\xi)$ the associated fluid gamma factor.

The second estimate is $\overline{U}_{\text{f}}^{\mathrm{Esp}}=\sqrt{\frac{3}{4}\kappa_{\text{v}}\alpha}$, where the function $\kappa_{\text{v}}(v_{\text{w}},\alpha)$ is given in the Appendix of Espinosa et al. Espinosa et al. (2010), using $v_{\text{w}}$ extracted from 1D simulations at $t=7000/T_{\text{c}}$, and $\alpha=\alpha_{\theta}(T_{\text{n}})$ from Table 1. Note that $\kappa_{\text{v}}$ is defined from the trace of the spatial part of the energy-momentum tensor as

[TABLE]

We see that the rotational component of the velocity field is very small, consistent with the interpretation of the fluid flow as sound waves, and with the linearity of the flow. We also see that the Espinosa et al. Espinosa et al. (2010) fitting formula for the mean square fluid velocity around a single bubble gives a good estimate of the mean square velocity of the 3D flow. Finally, the time at which the scalar gradient energy peaks is approximately $t_{\text{pc}}\simeq 0.6R_{*}/v_{\text{w}}$ in the simulations where $R_{*}$ is large. This is consistent with interpreting this peak as a bubble collision time.

In Figs. LABEL:f:VelPSDet, 4 and LABEL:f:VelPSDefJou we plot velocity power spectra (15) for a set of simulations. We divide the data according to the bubble wall speed $v_{\text{w}}$, separating out the special case of deflagrations moving at close to sound speed.

In Fig. LABEL:f:VelPSDet we show detonations with wall speeds $v_{\text{w}}=0.92$, $0.80$ and $v_{\text{w}}\simeq 0.7$ for transitions of weak and intermediate strength. The intermediate strength transition at $v_{\text{w}}=0.72$ is run at a higher resolution ($dx=1/T_{\text{c}}$) to resolve the higher velocity gradients, so that the mean bubble separation ($1889/T_{\text{c}}$) is close to the box size of the simulation ($2100/T_{\text{c}}$), hence there is less dynamic range on the long-wavelength side of the peak in the power spectrum.

The general form is a broken power law, with a domed peak at $kR_{*}=\textrm{O}(10)$. The shape is similar between the weak and intermediate cases at the same velocity; the intermediate strength transitions have higher amplitude, as more energy is transferred into kinetic energy. This is in accord with the sound shell model Hindmarsh (2016), where the velocity power spectra are of a universal shape for a given wall velocity and non-relativistic fluid flows with negligible shocks.

The width of the dome is larger for the detonations with wall speeds closer to the speed of sound, with the peak displaced to the right relative to the fast detonation. This is also consistent with the sound shell model prediction that the peak position in the power spectrum is determined by the inverse width of the sound shell.

The power-law to the right of the dome is close to the $k^{-1}$ predicted by the sound shell model, particularly in the detonations, where there is a clear separation between the peak scale and the wall width scale. The long-wavelength power law index is not so clear, as there are few bins and there are fewer $k$ vectors in each bin, but is consistent with the predicted $k^{3}$.

In Fig. 4 we show velocity power spectra from deflagrations with wall speed $v_{\text{w}}=0.44$. The power law to the right of the dome appears steeper than the $k^{-3}$ prediction, and there is a knee in the power spectrum at higher $k$ neither of which is in accord with the sound shell model. These features need further investigation.

We recall that in the sound shell model, there are two scales in the velocity field: the mean bubble separation $R_{*}$, and the sound shell width $\Delta R_{*}=R_{*}\Delta v_{\text{w}}/v_{\text{w}}$. If the scales are well separated (i.e. if the bubble wall is moving close to the speed of sound), the long-wavelength $k^{3}$ power law is predicted to turn into a $k^{1}$ power law at $kR_{*}=\textrm{O}(1)$, and finally to a $k^{-1}$ power law at $k\Delta R_{*}=\textrm{O}(1)$.

The clearest scale separation should be found in transitions where the wall speed is closest to the sound speed. In Fig. LABEL:f:VelPSDefJou we show the power spectra from a deflagration with speed $v_{\text{w}}=0.56$, very close to $c_{\text{s}}=1/\sqrt{3}\simeq 0.577$.

We do not have the dynamic range to resolve all three wavelength ranges simultaneously, but by altering the number of bubbles we can try to resolve two ranges at a time. At the top, for the largest $R_{*}$, we see that the dome around the peak has broadened into a slowly rising plateau, consistent with a $k^{1}$ behaviour. At higher wavenumber the plateau drops off, although we do not have enough range to confirm a $k^{-1}$ behaviour.

For larger numbers of bubbles (centre, bottom in Fig. LABEL:f:VelPSDefJou) the long-distance behaviour emerges. With $N_{\text{b}}=5376$ ($R_{*}=480/T_{\text{c}}$) one can clear see a steep power law, even steeper than $k^{3}$, in the range $1\lesssim kR_{*}\lesssim 3$. A possible reason for the discrepancy with the generic sound shell model prediction at low $kR_{*}$ is that the coefficient of the $k^{3}$ term is proportional to $\Delta v_{\text{w}}^{2}\simeq 4\times 10^{-4}$ Hindmarsh (2016), and so the next order term in a series expansion in $kR_{*}$ can dominate for $kR_{*}=\mathrm{O}(1)$. The short-distance behaviour at these large bubble numbers is not reliable, as there is insufficient distinction between the sound shell width $\Delta R_{*}\simeq 15/T_{\text{c}}$ and the bubble wall width $\ell\simeq 3/T_{\text{c}}$. Indeed, in Fig. 2 one can see that the fluid velocity profiles are far from their asymptotic forms.

V Results: gravitational waves

In Table 4 we show global quantities computed from the gravitational wave power spectrum

[TABLE]

where $z=kR_{*}$. These are the integral scale of the gravitational waves

[TABLE]

calculated at the end of the simulations, and the dimensionless gravitational wave amplitude parameter defined in (24), with the fluid length scale $L_{\text{f}}$ is taken to be either the integral scale of the velocity field $\xi_{\text{f}}$ (defined analogously to Eq. (33)) or the mean bubble separation $R_{*}$. With either of these length scales, $\tilde{\Omega}_{\text{gw}}$ is approximately constant and of order $10^{-2}$ for weak transitions444Note that there is a numerical error in the computation of $\tilde{\Omega}_{\text{gw}}$ in Hindmarsh et al. (2015) which resulted in a value which was a factor $(2\pi)^{3}/32\pi$ too high.. The integral scales of both the fluid and the gravitational waves are significantly smaller for the near-Jouguet transitions, showing the influence of the sound shell thickness.

In Figs. LABEL:f:GWPSDet to LABEL:f:GWPSDefJou we plot power spectra of the fractional energy density in gravitational waves. The spectra are divided by the mean bubble separation $R_{*}$ and the time $t$ in units of the Hubble distance and time, and plotted against wavenumber in units of the inverse bubble separation, for ease of comparison with Eq. (21). Taking the fluid flow length scale to be $R_{*}$, we have

[TABLE]

Note that by dividing by time, gravitational wave power generated in the collision phase will decrease, while acoustic phase gravitational wave power will asymptote to a constant. Note also that, by dividing by $H_{\text{n}}^{2}$, we arrive at a quantity which is independent of $G$.

One can see that, at late times, the shape of the power spectrum appears to change little, and is settling down to a characteristic shape. We would expect power laws to be less clear in the gravitational wave power spectrum than in the velocity power spectrum, as the former is a convolution of the latter over a range $\Delta k=\pm c_{\text{s}}k$ at wavenumber $k$ Hindmarsh et al. (2015).

Where power laws are established over a sufficient range, we expect that a velocity power spectrum going as $k^{n}$ should produce a gravitational wave power spectrum going as $k^{2n-1}$. The sound shell model Hindmarsh (2016) predicts $n=-1$ for $R_{*}/\Delta R_{*}\ll kR_{*}\ll R_{*}/\ell$, $n=1$ for $1\ll kR_{*}\ll R_{*}/\Delta R_{*}$ (if the scales $R_{*}$ and $\Delta R_{*}$ are well separated), and $n=3$ for $kR_{*}\ll 1$.

In Fig. LABEL:f:GWPSDet we show the gravitational wave power spectra from detonations ($v_{\text{w}}=0.92$, $0.80$, and $v_{\text{w}}\simeq 0.7$) for transitions of weak and intermediate strength, arising from the velocity fields with power spectra in Fig. LABEL:f:VelPSDet. Again, the general form is a broken power law, with a domed peak at $kR_{*}=\textrm{O}(10)$. As with the velocity power spectrum, the shape is similar between the weak and intermediate cases at the same wall velocity; the intermediate strength transitions have higher amplitude, resulting from the higher RMS velocity. We see that the broader dome in the velocity power spectrum translates to a broader feature in the gravitational wave spectrum, although the features around the peak are not well resolved, particularly in the intermediate strength transition.

The power-law to the right of the dome is close to the $k^{-3}$ predicted by the sound shell model, particularly in the the detonations. The predicted $k^{5}$ power law at long-wavelengths is difficult to discern, as it is buried under a feature established during the collision phase.

In Fig. 7 we show gravitational wave power spectra from a deflagrations with $v_{\text{w}}=0.44$, for weak and intermediate strength transitions. The power law at high $k$ is steeper than $k^{-3}$, inconsistent with the sound shell model, but consistent with velocity power spectra steeper than $k^{-1}$.

In Fig. LABEL:f:GWPSDefJou we show the gravitational wave power spectra from a deflagration with speed $v_{\text{w}}=0.56$, where the sound shell is very thin, and there are clearly two scales in the power spectra. Again, the low-$k$ behaviour is hidden behind the gravitational waves from the collision phase due to the limited duration of the simulation, although there is a suggestion of a steepening below $kR_{*}\sim 5$ in the case with the maximum long-wavelength resolution, $N_{\text{b}}=5376$. It is clear that the peak is at around $kR_{*}\simeq 50$, which is understandable in terms of the scale $k\Delta R_{*}\simeq 2$.

VI Modelling the power spectrum for LISA

Our simulations nucleate all bubbles simultaneously, whereas in a real thermal phase transition away from metastability, the nucleation rate rises exponentially as $p(t)=p_{0}\exp[\beta(t-t_{\text{f}})]$ after the temperature drops below the critical temperature, where $\beta$ is the transition rate parameter, and $t_{\text{f}}$ is the time at which the volume fraction of the symmetric phase is $1/e$ Enqvist et al. (1992). This means that there are a few larger, bubbles earlier in the transitions, and more smaller bubbles as the transition ends. Numerical experiments with the gravitational wave power spectrum in the envelope approximation Weir (2016) show that the principal effect of instantaneous nucleation for a given mean bubble separation is to increase the peak frequency by a factor of approximately 1.7, and decrease the amplitude by a factor of about 3. The shape of the power spectrum is not significantly changed. This rescaling can be regarded as a rescaling of the relationship between the mean bubble separation $R_{*}$ and the transition rate parameter $\beta$.

We shall assume that the same is true for the velocity field generated by the colliding bubbles, i.e. that the shape is not significantly changed by using a realistic nucleation history, and that the principal effect is to change the proportionality constant in the standard relation $R_{*}=(8\pi)^{\frac{1}{3}}v_{\text{w}}/\beta$ Enqvist et al. (1992). With this assumption we can directly use our measured power spectra as models for the gravitational wave power spectrum from a phase transition.

In the LISA Cosmology Working Group report Caprini et al. (2016) the acoustic gravitational wave power spectrum was modelled using a broken power law function:

[TABLE]

where

[TABLE]

and

[TABLE]

The dimensionless parameters $A$ and $(kR_{*})_{\rm max}$ determine the magnitude and the location of the maximum of the power spectrum, respectively. The form of the function is motivated by the results from hydrodynamical gravitational wave production simulations in Ref. Hindmarsh et al. (2015). The power spectrum of the ansatz at small $k$ is $\propto k^{3}$, turning over to $\propto k^{-4}$ at large $k$.

The detailed structure of the ansatz determines the width of the dome between the small-$k$ and large-$k$ regions; indeed, the width of the dome can be adjusted by generalising the ansatz to a form

[TABLE]

In order to enable direct comparison with the working group results we fit the measured power spectrum using the original ansatz (35).

The resulting fits are shown as dashed lines in Fig. LABEL:f:GWPSDet for detonations and in Fig. 7 for deflagrations, and fit parameters are listed in Table 5. The functions are fitted only in the neighbourhood of the domes, but at least for the weak detonations the fitted functions describe the overall behaviour of the data surprisingly well.

By integrating (35) and comparing to Eq. (23), we can derive an estimate for the amplitude parameter from our simulations555There is a missing factor of 3 in this expression, see the erratum at the end of this paper.,

[TABLE]

The estimates, computed from the values of $\overline{U}_{\text{f}}$ in Table 3 and $\tilde{\Omega}_{\text{gw}}^{R_{*}}$ in Table 4, are shown in the last column of Table 5.

It can be seen that the amplitude estimates based on numerical integration of the scaled gravitational wave power spectrum $\tilde{\mathcal{P}}_{\text{gw}}$ are generally higher than those derived from the fit: this is because the numerical power spectrum exaggerates the low-$k$ part of the power spectrum derived from the collision phase. In the case of the intermediate strength transition at $v_{\text{w}}=0.72$, the dome is less apparent, and the fitting formula under-predicts the gravitational wave spectrum at high $k$. We recommend using the fitting formula for wall velocities away from the Chapman-Jouguet speed by about $|v_{\text{w}}-v_{\text{CJ}}|\gtrsim 0.1$.

For the near-Jouguet deflagration ($v_{\text{w}}=0.56$) the dome in the power spectrum becomes very broad, as can be seen in Fig. LABEL:f:GWPSDefJou. In this case the fit ansatz (35) cannot describe the behaviour well. It is possible to construct more complicated fit functions which can capture the structure at intermediate scales, but the limited numerical data makes it difficult to see universal features. We leave the detailed analysis of the Jouguet case for further analysis.

The peak angular frequency in units of the mean bubble separation $z_{\text{p}}=(kR_{*})_{\rm max}$ is generally around 10, except near the Chapman-Jouguet speed where it is larger. Qualitatively this agrees with the estimate for the peak frequency made in Ref. Caprini et al. (2016) $f_{\text{p}}\simeq 1.2\beta/v_{\text{w}}$, with the phase transition rate parameter $\beta\simeq 3v_{\text{w}}/R_{*}$. For a more precise estimate we must take into account the fact that nucleating bubbles simultaneously changes the effective transition rate for a given $R_{*}$ Weir (2016), as outlined above. A power spectrum peaking at $(kR_{*})=z_{\text{p}}$ in our simulations corresponds to a true peak frequency

[TABLE]

where $S\simeq 2$ is a factor which takes into account the overestimate of the frequency for a given $R_{*}$.

In order to calculate the observed frequency of waves emitted with frequency $f_{\text{p}}$ (40) at time $t_{\text{n}}$ when the Hubble parameter was $H_{\text{n}}$, we note that the peak frequency today can be written

[TABLE]

Here

[TABLE]

is the Hubble rate at the nucleation temperature, redshifted to today, assuming that the dominant source of energy density is radiation. Hence

[TABLE]

To obtain the amplitude of the gravitational wave power spectra today, the power spectrum (21) and total power (23) must be multiplied by a factor $F_{\text{gw},0}=\Omega_{\gamma,0}(h_{0}/h_{*})^{\frac{1}{3}}(h_{0}/2)$, where $\Omega_{\gamma,0}$ is the density parameter of photons today, $h_{0}$ is the effective number of relativistic degrees of freedom contributing to the entropy today, and $h_{*}\simeq g_{*}$ is the corresponding number at the time of gravitational wave generation.

Using the Planck best-fit value $H_{0}=67.8\pm 0.9\,\textrm{km}\,\textrm{s}^{-1}\,\textrm{Mpc}^{-1}$ Ade et al. (2016), and the Far Infrared Absolute Spectrophotometer (FIRAS) temperature for the Cosmic Microwave Background $T_{\gamma,0}=2.725\pm 0.002\,\textrm{K}$ Mather et al. (1999), we have

[TABLE]

Our final expression for the acoustic gravitational wave power spectrum today is666The numerical prefactor in this expression contains some errors; see the erratum at the end of this paper.

[TABLE]

In Fig. 8 we show the signal-to-noise ratio expected at a LISA-like gravitational wave observatory with 6 laser links of arm length $2$ Gm and a 5-year mission duration Klein et al. (2016) (intermediate between configurations C1 and C2 of the LISA Cosmology Working Group report Caprini et al. (2016)). The SNR is shown as contours in the $(\overline{U}_{\text{f}},H_{\text{n}}R_{*})$ plane, for a phase transition with a reference temperature of $T_{\text{n}}=100\,\textrm{GeV}$, and taking999This should read $\tilde{\Omega}_{\text{gw}}=1.2\times 10^{-2}$. See the erratum. $\tilde{\Omega}_{\text{gw}}=1.2\times 10^{-1}$. This value reproduces the peak amplitude for the intermediate strength transition with $v_{\text{w}}=0.92$. It under-predicts the power spectrum for the intermediate strength transition at other wall speeds we simulated, and is therefore a conservative estimate.

We also show contours of $R_{*}H_{\text{n}}/\overline{U}_{\text{f}}$, which gives the timescale for the appearance of shocks relative to the Hubble time. We do not expect the acoustic gravitational wave power spectrum to be an accurate description for $R_{*}H_{\text{n}}/\overline{U}_{\text{f}}\ll 1$.

From the figure, we can conclude that, for a phase transition at 100 GeV, LISA’s peak sensitivity will be at a mean bubble separation of about a hundredth the Hubble length, at which scale an intermediate strength deflagration or detonation ($\overline{U}_{\text{f}}\simeq 5\times 10^{-2}$) should produce a SNR of around 10, taken as the detection threshold in Caprini et al. (2016). However, this is also in the region where one cannot safely assume that the fluid has not developed shocks. Further simulations are required for a more accurate determination of LISA’s ability to explore the parameter space of phase transitions.

VII Conclusions

We have performed the largest numerical simulations to date of first order phase transitions in the early universe, and computed the resulting fluid velocity and gravitational wave spectra for a range of bubble wall speeds $v_{\text{w}}$, mean bubble separations $R_{*}$, and transition strengths $\alpha_{w}$. The power spectra are more tightly pinned down than in our previous campaign of simulations Hindmarsh et al. (2015).

We observe a gravitational wave power spectrum with a rising power law, a broad dome at $kR_{*}=\textrm{O}(10)$, and a decreasing power law at higher $k$. The dome widens to a slowly rising plateau for bubble wall speeds close to the speed of sound.

In the case of detonations the spectra exhibit a $k^{-3}$ power law in good agreement with the sound shell model Hindmarsh (2016), but the high-$k$ power law for deflagrations appears slightly steeper.

We establish that there are two length scales in the gravitational wave power spectrum: besides the mean bubble separation $R_{*}$ there is also the mean sound shell thickness $\Delta R_{*}=R_{*}\Delta v_{\text{w}}/c_{\text{s}}$. For a thin sound shell, where $\Delta R_{*}\ll R_{*}$, the slowly rising plateau has a power law consistent with the $k^{1}$ predicted by the sound shell model.

We do not have sufficient computational volume to determine the gravitational wave power spectrum at $kR_{*}\lesssim\textrm{O}(1)$, but the velocity power spectra for generic wall speeds are consistent with $k^{3}$ there, which should lead to $k^{5}$ in the gravitational waves. In the special case of a just-subsonic deflagration, where the sound shell is very thin, the velocity power spectrum is steeper than $k^{3}$ for $kR_{*}\lesssim\textrm{O}(1)$. This is not in contradiction with the sound shell model, which only gives $k^{3}$ for $kR_{*}\lesssim\textrm{O}(|v_{\text{w}}-c_{\text{s}}|)$.

Our simulations are still not large enough to properly resolve all the different wavenumber regimes simultaneously, and they reveal an interesting knee in the velocity power spectra for deflagrations at high $k$. This feature appears to be established at about the peak bubble collision time $t_{\text{pc}}$, where the area of the phase boundary is at its maximum. Yet larger simulations are therefore needed to resolve the full power spectrum, and a separate simulation campaign to investigate the deflagration power spectra.

We also need to run the simulations for longer, with larger fluid velocities, in order to investigate the transition to turbulence. Flows with larger fluid velocities become turbulent earlier, and are likely to be important for the gravitational wave signals observable by LISA.

However, for transitions which are weak enough that turbulence does not develop, and have bubble wall speeds not too close to the speed of sound, we are confident in the form of the power spectrum, and offer the fitting formula (35) for the gravitational wave power spectrum. This can safely be used for RMS fluid velocities up to $\overline{U}_{\text{f}}\lesssim R_{*}H_{\text{n}}$, or $\kappa_{\text{v}}\alpha_{\theta}\lesssim(H_{\text{n}}/\beta)^{2}$ in terms of the kinetic energy conversion efficiency $\kappa_{\text{v}}$ and the transition strength parameter $\alpha_{\theta}$ (see Eq. (31)).

Using the fitting formula, we computed the signal-to-noise ratio at a LISA-like mission with 6 laser links, arm length 2 Gm, and a duration of 5 years. At a transition temperature of $T_{\text{n}}=100\;\textrm{GeV}$, we find that such a mission has greatest sensitivity to transitions with mean bubble separation of about one hundredth of the Hubble length. At this bubble separation, LISA will be able to detect the signal from phase transitions down to latent heat to energy ratio $\alpha_{w}\simeq 0.1$, i.e. intermediate strength in our terminology.

We also checked the timescale for the development of shocks and turbulence, finding that a significant fraction of the detectable parameter space is in a region where one cannot safely assume that the fluid flow remains linear. Further simulations are required to examine the transition to turbulence, necessary for a more accurate determination of LISA’s ability to explore the parameter space of phase transitions.

There is a range of extensions of the standard model (SM), which could produce a gravitational wave signal in reach of LISA Caprini et al. (2016). Most of these models are extensions of the SM Higgs sector by additional scalar fields, principally SU(2) singlets or doublets, or modifications of the SM Higgs potential itself. Furthermore, LISA could also probe phase transitions in the TeV range, e.g. possible confinement transitions in strong coupling completions of the electroweak theory.

A strong phase transition in these models can either come from 1) a modified zero temperature Higgs potential, as is often the case in singlet extensions, or 2) additional thermal contributions to the Higgs potential, or 3) a combination of both, as e.g. in the Two-Higgs-Doublet model (2HDM). Models of type 2) and 3) will show a stronger dependence of the nucleation rate on temperature (via the energy of the critical bubble), which results in a smaller value of $\beta/H$. So models of type 1) will have larger bubbles at fixed $\alpha$ compared to types 2) and 3). As a result models of type 1) may show a stronger gravitational wave signal without reaching the state of turbulence.

Most of the models will show an observable gravitational wave signal signal only when the walls propagate as fast detonations, but there may be exceptions as has been shown e.g. for the 2HDM Dorsch et al. (2016).

More generally, the fact that the shape of the gravitational wave power spectrum depends on the bubble wall speed means that it could become possible to determine $v_{\text{w}}$ from a detection of the stochastic gravitational wave background with sufficiently high signal-to-noise ratio.

Delineating the ability of future space-based gravitational wave detectors to constrain the full set of phase transition parameters $\alpha$, $\beta$, $v_{\text{w}}$ and $T_{\text{n}}$ of an electroweak-scale first order transition, thereby opening a new quantitative window onto physics beyond the Standard Model, is now an important goal.

Acknowledgements.

The LISA Cosmology Working Group has provided an important forum for discussing our work and we thank its members as well as its coordinators, Chiara Caprini and Germano Nardini. We owe a huge debt also to the late Pierre Binétruy for his work on cosmological sources of gravitational waves, his work with LISA, and in particular his encouragement of our work. We acknowledge PRACE for awarding us access to resource HAZEL HEN based in Germany at the High Performance Computing Center Stuttgart (HLRS). Our work was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”. MH (ORCID ID 0000-0002-9307-437X) and SH acknowledge support from the Science and Technology Facilities Council (grant number ST/L000504/1). KR (ORCID ID 0000-0003-2266-4716) is supported by the Academy of Finland grant 267286. DJW (ORCID ID 0000-0001-6986-0517) was supported by the People Programme (Marie Skłodowska-Curie actions) of the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement number PIEF-GA-2013-629425.

Erratum

There is a factor of 3 missing from the right hand side of Eq. 39, which should read:

[TABLE]

Without this additional factor of 3 this expression is not consistent with Eqs. 21 and 23, nor with Ref. Hindmarsh et al. (2015). However, the values of $A_{\text{est}}$ in Table 5 are correct and include this factor.

This factor of 3 is, in turn, missing from Eq. 45, which should read (note that original equation also truncated the numerical coefficient as 0.68 rather than 0.687):

[TABLE]

We also wish to highlight that this formula is written without the customary factors of the reduced Hubble constant $h^{2}$ on each side. This itself does not affect the validity of the equation and indeed allows us to select a value of $h$ for comparison with a specific precomputed gravitational wave energy density sensitivity curve.

Note also that in the paragraph below Eq. 45 the exponent in the value for $\tilde{\Omega}_{\text{gw}}$ is incorrect. The expression should read $\tilde{\Omega}_{\text{gw}}=1.2\times 10^{-2}$.

The SNR contours in Fig. 8 were produced using Eq. 45. The above factor of 3 (mistakenly absent from Eq. 45) was included but the factor of 0.68[7] (present in Eq. 45) was excluded, and there was an additional erroneous factor of 2. This omission means the correct predicted SNR is a factor of $0.687\times 0.5\approx 0.34$ lower than in the original plot.

Figure 1 below is a corrected version of Fig. 8 that takes account of the numerical errors explained above.

However, the mission profile of LISA has changed substantially since this paper was published, with a longer expected mission duration and a modified sensitivity curve. This will further change the predicted SNR curves. The plots in the following referenced paper are free of the errors discussed in the erratum, and use the most up-to-date mission profile101010The parameter space can also be explored at ptplot.org..

Acknowledgements.

We are grateful to Aleksandr Azatov, Francesco Sgarlata and Daniele Barducci for bringing this to our attention.

Reference

C. Caprini et al., “Detecting gravitational waves from cosmological phase transitions with LISA: an update,” JCAP 2003, 024 (2020), [arXiv:1910.13125 [astro-ph.CO]].