\section{Introduction}

Gravitational waves (GWs) promise a new and exciting window to the cosmos. Two ground-based interferometer experiments, LIGO and VIRGO,  are about to restart operations with greatly increased sensitivity \cite{Harry:2010zz,Accadia:2009zz}, and will be joined in a few years by KAGRA \cite{Somiya:2011np}.
Once working at their design sensitivity, they are expected to quickly detect gravitational wave signals from binary neutron stars \cite{Abadie:2010cf}. 
Gravitational waves also offer a unique way to learn about the early universe. A range of phenomena, such as inflation, topological defects and phase transitions may lead to observable gravitational wave signals across a wide range of frequencies (for a review see \cite{Binetruy:2012ze}). 
There are a number of proposals to realise a gravitational wave detector in space, in the first place eLISA \cite{Seoane:2013qna}, which is scheduled for launch in 2034. Space-based detectors have much longer arm lengths than ground based ones, and have maximum sensitivity in a frequency range which is relevant for a first order phase transition at the electroweak scale.

Given these exciting observational prospects, we revisit the generation of gravitational waves in first order thermal phase transitions in the early universe. 
We have in mind
an electroweak-scale phase transition, but nothing in our formalism is specific to electroweak scale physics.  In the Standard Model the electroweak transition is known to be a cross-over \cite{Kajantie:1995kf,Kajantie:1996mn,Gurtler:1997hr,Csikor:1998eu,DOnofrio:2014kta}, which does not lead to a gravitational wave signal.  However, a strong first order phase transition is possible in various extensions of the Standard Model~\cite{Carena:1996wj,Delepine:1996vn,Laine:1998qk,Grojean:2004xa,Huber:2000mg,Huber:2006wf,Dorsch:2013wja}.

We reduce the original physical system to a model consisting of a scalar order parameter field coupled to an ideal fluid. The parameters of the model can in principle be fixed by matching to the thermodynamical quantities of the original theory.  We perform very large scale numerical simulations  to determine the fluid and gravitational wave power spectra. The ultimate goal is to understand what information on the phase transition can be extracted from the future observation of a gravitational wave signal. 

Since the early nineties there have been a number of studies of gravitational waves from phase transitions. In Refs.~\cite{Turner:1990rc,Kosowsky:1991ua,Kosowsky:1992rz,Kosowsky:1992vn}, the case of a scalar field only, i.e.\ a vacuum transition without fluid, was considered, motivated by models of inflation terminated by a first order transition.  Vacuum transitions during inflation and with a fluid were considered in Ref.\ \cite{Chialva:2010jt}.

In a vacuum transition, all the energy released goes into the bubble wall, which as a result is accelerated to the speed of light.
After solving numerically the field equations for the collision of two scalar bubbles \cite{Hawking:1982ga,Kosowsky:1991ua}, it was realised that the energy-momentum tensor sourcing the gravitational wave production can be approximated by the envelope of the colliding bubbles \cite{Kosowsky:1992rz,Kosowsky:1992vn}. This ``envelope approximation'' models a configuration of expanding bubbles by the overlap of a corresponding set of infinitely thin shells. The envelope disappears once the transition is completed and gravitational wave production stops. It is found that the gravitational wave spectrum peaks at a frequency determined by the average bubble size at collision. In the UV, the spectrum falls as a power law, subsequently shown to be $k^{-1}$ \cite{Huber:2008hg}, where $k$ is the wave number. 
Numerical studies in Ref.~\cite{Child:2012qg} did not have the dynamic range to clearly confirm this behaviour, but the larger simulations done for the present work show some supporting evidence.

The case of a thermal phase transition, where the scalar field is coupled to a fluid, is more complicated. The nucleated bubbles will show accelerated expansion until the pressure inside is balanced by friction caused by the plasma. The bubbles then expand with a constant velocity. This is because the energy released by the transition grows with the volume of the bubble, i.e. $\sim R^3$, while the energy transferred to the scalar bubble wall only grows with the bubble surface, i.e. $\sim R^2$, where $R$ denotes the bubble radius. Hence only a tiny fraction of the released energy, on the order of the  ratio of the initial to final bubble radius, stays in the scalar field. In the case of a first order electroweak scale thermal phase transition this ratio is about $10^{-4}M_W/M_{\rm Pl}~\sim10^{-13}$. Therefore gravitational wave production in thermal phase transitions is completely dominated by the fluid.\footnote{An exception may be the case where the bubble wall ``runs away'', i.e. friction is not sufficient to prevent the wall from approaching the speed of light~\cite{Bodeker:2009qy}, similar to a vacuum transition. Then both the scalar and the fluid could contribute sizeably to the generated gravitational wave signal.} The energy which is released into the fluid mostly goes into reheating the plasma. A small and calculable fraction \cite{Kamionkowski:1993fg,Espinosa:2010hh} goes into bulk motion of the fluid and can source gravitational waves. 

Having established the fluid as the main source of gravitational waves, the question of the production mechanism arises. Several mechanisms have been suggested and studied in the literature. In the simplest approach, one assumes that the fluid put into motion by the scalar wall can still be treated as a thin shell and the energy momentum tensor sourcing gravitational wave production can again be approximated by the shell overlaps~\cite{Kamionkowski:1993fg}.
In this case gravitational wave production finishes with the completion of the phase transition, and a characteristic prediction is the $k^{-1}$ UV power law of the spectrum.
Another possibility is that the collision of bubbles induces turbulent motion of the fluid \cite{Kamionkowski:1993fg}. The resulting eddies generate gravitational waves even after the transition is completed~\cite{Kamionkowski:1993fg,Caprini:2006jb,Gogoberidze:2007an,Caprini:2009yp}. Various UV power laws of the the gravitational wave spectrum have been suggested in this context, such as $k^{-3.5}$~\cite{Dolgov:2002ra} and $k^{-8/3}$~\cite{Caprini:2006jb}.

To shed light onto these competing scenarios, we recently performed large scale numerical simulations of a thermal phase transition of a scalar field plus fluid system~\cite{Hindmarsh:2013xza}. We found no indications that fluid turbulence was an important source of gravitational radiation. Instead sound waves are generated by the explosive bubble growth, which propagate through the plasma until long after the transition is completed. In our simulations these sound waves are the dominant source of gravitational waves.  After the phase transition, the fluid energy-momentum tensor  clearly does not show the form assumed in the envelope approximation. The nearly linear behaviour of sound waves is very different to the highly nonlinear behaviour of the scalar field. 

Other numerical simulations of the generation of gravitational waves by the coupled field-fluid system, using an explicit update algorithm for the fluid, have been described recently in Ref.~\cite{Giblin:2014qia}. 
The generation of gravitational waves through sound in QCD and electroweak phase transitions was also recently studied in Ref.~\cite{Kalaydzhyan:2014wca}, with special focus on the effect of possible non-linear sound dispersion relations, which were argued to lead to an inverse acoustic cascade. In Ref.~\cite{Ghiglieri:2015nfa}, generation of gravitational waves in the Standard Model in the absence of such a cascade was discussed.

In the present work we simulate at larger volumes and larger average bubble separations than in \cite{Hindmarsh:2013xza}, for the same range of bubble wall speeds and phase transition strengths. We widen the dynamic range even more by nucleating all bubbles at the same time. 
We confirm that the gravitational wave density parameter is proportional to the fourth power of the mean square fluid velocity, the ratio of lifetime of the source to the Hubble time, and the ratio of length scale of the source to the Hubble length. 
We measure the length scale of the source, approximately the average bubble separation in \cite{Hindmarsh:2013xza}, directly from the fluid flow. With this improvement,  the proportionality constant for the gravitational wave density parameter varies much less between phase transitions with different strengths and bubble wall speeds. Our measurements show that it is $0.8 \pm 0.1$, where the uncertainly is the root mean square fluctuation between simulations.


We show that the resulting gravitational wave spectrum exhibits UV power laws which are clearly steeper than the $k^{-1}$ predicted by the envelope approximation. In the case of deflagrations (where the bubble walls are subsonic), we are reasonably confident that the power law is $k^{-3}$. For detonations we do not have sufficient dynamic range to be certain of the power law index. 

We compare the acoustic gravitational waves with the standard prediction from the envelope approximation. 
We argue that the envelope approximation is based on an incorrect picture of the dynamics of the fluid, in which the fluid perturbations are destroyed by bubble collisions in the same way as the bubble walls. Instead, they pass through one another, and keep oscillating, resulting in a gravitational wave source whose effective lifetime is the
Hubble time.
The true gravitational wave energy density is therefore a factor $\beta/\Hc$ higher, where $\Hc$ is the Hubble rate at the phase transition, 
and $\beta^{-1}$ is the duration of the phase transition. 
For a thermal electroweak-scale phase transition, the gravitational wave signal is larger than hitherto believed by at least two orders of magnitude.

\section{First order phase transitions in cosmology}


\subsection{Hydrodynamics}

We describe the phase transition using the cosmic fluid -- order parameter field model~\cite{Ignatius:1993qn,KurkiSuonio:1995vy}, which we summarise here.  The model contains a classical scalar field $\phi$ (effective order parameter), which is coupled to ideal fluid hydrodynamics.  The variables describing the local state of the matter are local temperature $T$, fluid 4-velocity $U^\mu$ and the scalar order parameter field $\phi$.  The first order dynamics are obtained by introducing a temperature dependent effective potential $V(\phi,T)$.  
Following \cite{Enqvist:1991xw,Ignatius:1993qn}, we use a simple $\phi^4$ form for the potential:
\begin{equation}
\label{e:ScaPot}
V(\phi, T) = \frac{1}{2} \gamma (T^2-T_0^2) \phi^2 - \frac{1}{3} A T \phi^3 + \frac{1}{4}\lambda\phi^4.
\end{equation}
The detailed form of the potential is not important, as long as it allows for a first order phase transition of sufficient strength.  A first order transition occurs if $2A^2 < 9\lambda\gamma$.

The equation of state of the coupled scalar field and fluid system is
\begin{align}
\epsilon(T,\phi) &= 3 a T^4 + V(\phi,T) - T\frac{\partial V}{\partial T},\\
  p(T,\phi) &= a T^4 - V(\phi,T)
\end{align}
where $a=(\pi^2/90)g_*$, and $g_*$ is the effective number of degrees of freedom.
The latent heat density (usually just called the latent heat) is
\begin{equation}
{\cal L}(T) = w(T,0) - w(T,\phi_\text{b})
\end{equation}
where $w = \epsilon + p$ is the enthalpy density, and $\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 characterised by the ratio of the latent heat to the total radiation density in the high temperature symmetric phase,
\begin{equation}
\strengthPar{T} = \frac{{\cal L}(T)}{3aT^4}.
\label{e:StrParA}
\end{equation}
The total energy-momentum tensor of the system can be written as
\begin{equation}
\label{eq:tmunu}
T^{\mu\nu} = \partial^\mu \phi \partial^\nu \phi - {\textstyle \frac12} g^{\mu\nu} (\partial\phi)^2
+ \left[\epsilon + p \right] U^\mu U^\nu + g^{\mu\nu} p
\end{equation}
where the metric convention is (-+++).
The energy-momentum tensor is conserved, $\partial_\mu T^{\mu\nu} = 0$.
The interaction between field gradients and the fluid is introduced by
splitting the conserved current nonuniquely into field and fluid parts, 
which are then coupled together through a dissipative term proportional to field gradient:
\begin{align}
  [\partial_\mu T^{\mu\nu}]_{\rm field} &=
    (\partial_\mu\partial^\mu\phi) \partial^\nu\phi -
    \frac{\partial V}{\partial\phi} \partial^\nu \phi = 
    \delta^\nu 
    \label{t-munu-1}\\
    [\partial_\mu T^{\mu\nu}]_{\rm fluid} &=
   \partial_\mu[(\epsilon+p)U^\mu U^\nu] - \partial^\nu p +
    \frac{\partial V}{\partial \phi} \partial^\nu \phi = 
    -\delta^\nu,
    \label{t-munu-2}
\end{align}
where the coupling term is 
\begin{equation}
\label{e:CouTer}
\delta^\nu = \eta U^\mu \partial_\mu\phi \partial^\nu \phi
\end{equation}
with $\eta$
an adjustable friction parameter \cite{Ignatius:1993qn}. 
Equations analogous to Eqs.~(\ref{t-munu-1})-(\ref{t-munu-2}) can, at least in principle, be derived from field theory (see e.g.\ \cite{Moore:1995si,Konstandin:2014zta}), but
the simplified model here is adequate for parametrising the entropy
production~\cite{KurkiSuonio:1996rk}.

From Eqs.~(\ref{t-munu-1}) and~(\ref{t-munu-2}) we can derive the equations of
motion in a form suitable for numerical simulation.  For the field we obtain
\begin{equation}
- \ddot{\phi} + \nabla^2 \phi - \frac{\partial V}{\partial \phi} = \eta W (\dot{\phi} + V^i \partial_i \phi)\,,
\label{eqfield}
\end{equation}
where $W$ is the relativistic $\gamma$-factor and $V^i$ is the fluid 3-velocity,
$U^i = WV^i$.  For the fluid energy density $E=W\epsilon$, contracting $[\partial_\mu T^{\mu\nu}]_{\rm fluid}$ with $U_\nu$ gives
\begin{multline}
\dot{E} + \partial_i (E V^i) + p [\dot{W} + \partial_i (W V^i)]  - \frac{\partial V}{\partial \phi} W (\dot{\phi} + V^i \partial_i \phi) \\ = \eta W^2 (\dot{\phi} + V^i \partial_i \phi)^2.
\label{eqE}
\end{multline}
Finally, the equations of motion for the fluid momentum density $Z_i = W(\epsilon+p) U_i$ are
\begin{equation}
\dot{Z}_i + \partial_j(Z_i V^j) + \partial_i p + \frac{\partial V}{\partial \phi} \partial_i \phi  = -\eta W (\dot{\phi} + V^j \partial_j \phi)\partial_i \phi.
\label{eqZ}
\end{equation}
The implementation of Eqs.~(\ref{eqfield})-(\ref{eqZ}) on a discrete lattice is described in Section \ref{sect:numerical}.

The parameters of the potential in Eq.~(\ref{e:ScaPot}) are related to thermodynamic quantities at the phase transition: the critical temperature $T_c$, latent heat ${\cal L}(\Tc)$, surface tension $\sigma$ and the broken phase correlation length (which is also of order the bubble wall thickness) $\ell$~\cite{Enqvist:1991xw}
\begin{align}
  T^2_c &= \frac{T_0^2}{1-2A^2/(9\lambda\gamma)} \\
  {\cal L} &= \frac{A^2\gamma}{\lambda^2} T_0^2 T_c^2 \\
  \sigma &= \frac{2\sqrt{2}}{81}\frac{A^3}{\lambda^{5/2}} T_c^3 \\
  \ell^2 &= \frac{9\lambda}{2A^2} \frac{1}{T_c^{2}}.
  \label{e:PhaTraPar}
\end{align}
Due to supercooling, the phase transition (bubble nucleation) starts at temperature $T_N$, where $T_0 < T_N < T_c$.  We are mostly interested in the large supercooling (LSC) case, where $T_N$ is typically somewhere in the middle between $T_0$ and $T_c$.  
However, we emphasise that our focus in this work is not the nucleation of critical bubbles, which in a given microscopic theory is a thermal field theory problem and can be studied in perturbation theory 
or with numerical simulations~\cite{Moore:2000jw}.  In our simulations both the density of the initial bubbles 
and the nucleation temperature $T_N$ are set by hand. 


\subsection{Bubble nucleation}
\label{ss:BubNuc}

The phase transition proceeds by the nucleation and growth of bubbles of the broken phase~\cite{Linde:1978px,Linde:1981zj}.
Bubble nucleation occurs at an exponentially growing rate per unit volume below the critical temperature~\cite{Hogan:1984hx,Enqvist:1991xw},
\begin{equation}
p(t) \simeq \Gamma_0  e^{-S(\tN) + \be(t - \tN)},
\label{e:TunRat}
\end{equation}
where $-\be$ is the time derivative of the action of the critical bubble $S(t)$, and $\Gamma_0$ is a dimensional prefactor of order $\alpha_W^5 \Tc^4$ \cite{Moore:2000jw}, where $\alpha_W \approx 1/30$. The nucleation time $\tN$ can be defined to be the time at which the nucleation rate reaches one bubble per Hubble volume per Hubble time, or $p(\tN) = H^4(\tN)$. 

The tunnelling rate parameter $\be$ not only sets the timescale of the transition, but also 
the average separation between bubbles once the transition has completed, $\Rbc$. 
Having defined $\Rbc$ to be the inverse cube root of the number density of bubbles, 
it can be shown that \cite{Enqvist:1991xw}
\begin{equation}
\Rbc = (8\pi)^{\frac13} \frac{\vw}{\be}.
\label{e:RbcBet}
\end{equation}
Strictly, Eq.~(\ref{e:RbcBet}) applies only for detonations.  For deflagrations, one should take into account the suppression of the tunnelling rate ahead of the bubble wall, where the fluid is heated by the release of latent heat.  In this case we would expect $\Rbc \sim  {\cs}/{\be}$.

The important ratio ${\be}/{\Hc}$ (the transition rate relative to the Hubble rate) follows from simple considerations of the temperature of the transition \cite{Hogan:1984hx}.
One can straightforwardly argue that 
\begin{equation}
S(\tN) \sim 4\ln (m_{\text{P}}/\TN),
\end{equation}
and that for tunnelling in a thermal effective potential (\ref{e:ScaPot})
\begin{equation}
\frac{\be}{\Hc} \simeq \frac{2 S(\tN)}{(1 - \TN/\Tc)}.
\end{equation}
Hence, for a thermal electroweak-scale transition, the critical bubble action must be $\mathrm{O}(10^2)$, and the ratio $\be/\Hc$ must be at least $\mathrm{O}(10^2)$.

A detailed non-perturbative evaluation of the bubble nucleation rate in the standard model electroweak theory is presented in Ref.~\cite{Moore:2000jw}, using an unphysically small Higgs mass in order to ensure a first order phase transition.  In this case the critical bubble action was found to be $\approx 90$, and $\be/\Hc \approx 2\times 10^4$.
These are expected to be generic numbers for any first order thermal electroweak-scale transition.


\section{Theory of GW generation}
\label{s:TheGWGen}
\subsection{GW power spectrum definition}

A gravitational wave is a propagating mode of the transverse and traceless part of the metric perturbation, $h_{ij}$.
We are interested in calculating the gravitational wave energy density power spectrum, where the gravitational wave energy-momentum tensor is 
\begin{equation}
T_{\mu\nu}^\text{GW} = \frac{1}{32\pi G}\left< \partial_\mu h_{ij} \partial_\nu h_{ij} \right>.
\end{equation}
To this end, we define the spectral density of the time derivative of the metric perturbation $\SpecDen{\dot h}(\textbf{k},t)$ by 
\begin{equation}
\vev{\dot{h}_{ij}(\textbf{k},t) \dot{h}_{ij}(\textbf{k}',t) } = \SpecDen{\dot h}(\textbf{k},t) \debar3{\textbf{k}+\textbf{k}'}.
\end{equation}
The gravitational wave energy density power spectrum is then
\begin{equation}
\frac{d \rho_\text{GW}}{d \ln(k)} = \frac{1}{32\pi G} \frac{k^3}{2\pi^2} \SpecDen{\dot h}(\textbf{k},t).
\label{e:GWPowSpeDef}
\end{equation}

\subsection{GW power spectrum from fluid and field}

The source of gravitational waves is the transverse traceless part of the spatial components of the energy-momentum tensor.  
Given that we will be removing the trace anyway, it suffices to consider a source tensor
$\tau_{ij}  = \tau^\phi_{ij} + \tau^\text{f}_{ij}$, which is decomposed into fluid and field pieces according to 
\begin{equation}
\label{e:TauDef}
\tau^\phi_{ij} = \partial_i \phi \partial_j \phi, \quad\tau^\text{f}_{ij} = W^2 (\epsilon + p)V_i V_j.
\end{equation}
The physical metric perturbations are recovered in momentum space by 
applying the projector onto transverse, traceless symmetric rank 2 tensors:
\begin{equation}
\label{e:ProDef}
\lambda_{ij,kl}(\textbf{k}) = P_{ik}(\textbf{k}) P_{jl}(\textbf{k}) - \frac12 P_{ij}(\textbf{k})P_{kl}(\textbf{k}) 
\end{equation}
with 
\begin{equation}
P_{ij}(\textbf{k})  = \delta_{ij} - \hat{k}_i \hat{k}_j. 
\end{equation}
The particular solution for the gravitational wave is therefore
\begin{equation}
h_{ij} (\mathbf{k},t) = (16\pi G)\lambda_{ij,kl}(\textbf{k}) \int_0^t dt' \frac{\sin[k(t-t')]}{k} \tau_{kl}(\textbf{k},t'),
\label{e:GWsol}
\end{equation}
where we have assumed that the source vanishes for $t'<0$.


Using the fact that the fluid shear stress dominates the spatial parts of the energy-momentum tensor, we write 
\begin{widetext}
\begin{eqnarray}
\vev{\dot{h}^{ij}_\textbf{k}(t) \dot{h}^{ij}_{\textbf{k}'}(t) } &=&  (16\pi G)^2 \int^t_0 dt_1 dt_2 \cos[ k(t - t_1)] \cos[ k(t - t_2)] 
\lambda_{ij,kl}(\textbf{k})\vev{\tau_\text{f}^{ij}(\textbf{k},t_1)\tau_\text{f}^{kl}(\textbf{k}',t_2)}.
\end{eqnarray}
Introducing the unequal time correlator (UETC) of the fluid shear stress $\Pi^2$ \cite{Caprini:2009fx,Figueroa:2012kw} through 
\begin{equation}
\lambda_{ij,kl}(\textbf{k})\vev{\tau_\text{f}^{ij}(\textbf{k},t_1)\tau_\text{f}^{kl}(\textbf{k}',t_2)} = \Pi^2(k,t_1,t_2) \debar3{\textbf{k} + \textbf{k}'}
\end{equation}
and averaging over a period $T$, much longer than the periods of the gravitational waves of interest, we can 
write
\begin{equation}
\SpecDen{\dot h}(k,t) = (16\pi G)^2 \int^t_0 dt_1 dt_2 \frac{\cos[ k(t_1 - t_2)]}{2} \Pi^2(k,t_1,t_2).
\end{equation}
\end{widetext}
On dimensional grounds, we can write the UETC as 
\begin{equation}
\label{e:UETCmod}
\Pi^2(k,t_1,t_2) \simeq [(\bar\epsilon+\bar p)\fluidV^2]^2L_\text{f}^3\tilde\Pi^2
\end{equation}
where 
$\bar\epsilon$ and $\bar{p}$ are the spatially averaged energy density and pressure;
$\fluidV$ is the root mean square fluid velocity, 
defined through
\begin{eqnarray}
\label{e:fluidVdef}
(\bar\epsilon+\bar p) \fluidV^2 &=&  \frac{1}{{\mathcal V}}\int_{{\mathcal V}} d^3x\tau^\text{f}_{ii}, 
\end{eqnarray}
where ${\mathcal V}$ is the averaging volume;
$L_\text{f}$ is a characteristic length scale in the velocity field; and $\tilde\Pi^2$ is a dimensionless function of $k$, $t_1$ and $t_2$.
In Ref.~\cite{Hindmarsh:2013xza} we estimated that $L_\text{f}$ would be the mean bubble separation, but we will not make that assumption yet. We will see that we can understand the numerical results better if we extract the scale directly from the fluid velocity field in the simulations.

We also assumed that the UETC would be a function of $t_1 - t_2$ for times between the nucleation time $\tN$ and the lifetime of the velocity perturbations $\tau_\text{v}$, and that there is no separate timescale in the function $\tilde\Pi^2$, apart from that generated from $L_\text{f}$, the speed of sound $\cs$, and the speed of light.
With these assumptions, we can write the dimensionless UETC as a function of $kL_\text{f}$ and $z = k(t_1 - t_2)$, 
and the spectral density of $\dot h$ becomes
\begin{widetext}
\begin{equation}
\SpecDen{\dot h}(k,t) = \left[16\pi G(\bar\epsilon+\bar p)\fluidV^2\right]^2  t k^{-1} L_\text{f}^3 \int dz \frac{\cos(z)}{2} \tilde\Pi^2(kL_\text{f},z).
\label{e:SpeDenExpA}
\end{equation}
Note that one could follow through the same arguments for the scalar field, which would contribute in exactly an analogous manner
\begin{equation}
\SpecDen{\dot h}^\phi(k,t) = \left[16\pi G(\bar\epsilon+\bar p)\fieldV^2\right]^2  t k^{-1} L_\phi^3\int dz \frac{\cos(z)}{2} \tilde\Pi^2_\phi(kL_\phi,z),
\end{equation}
\end{widetext}
where 
\begin{eqnarray}
\label{e:fieldVdef}
(\bar\epsilon+\bar p) \fieldV^2 &=&  \frac{1}{{\mathcal V}}\int_{{\mathcal V}} d^3x\tau^\phi_{ii}, 
\end{eqnarray}
$L_\phi$ is a characteristic scale in the scalar field configuration, and $ \tilde\Pi^2_\phi$ is the dimensionless unequal time correlator of the scalar field shear stress tensor.   However, as explained in the introduction, the field contribution is negligible in most phase transitions.

Hence, putting together (\ref{e:GWPowSpeDef}) and (\ref{e:SpeDenExpA}), we may write
the gravitational wave energy density power spectrum as 
\begin{equation}
\frac{d \rho_\text{GW}}{d \ln(k)} = 8 \pi G\left[(\bar\epsilon+\bar p)\fluidV^2\right]^2  t L_\text{f} \frac{(kL_\text{f})^3}{2\pi^2} \tilde P_{\text{GW}}(kL_\text{f}),
\label{e:GWEneDenPowSpe}
\end{equation}
where
\begin{equation}
\tilde P_{\text{GW}}(kL_\text{f})  = \frac{1}{kL_\text{f}}\int dz \frac{\cos(z)}{2} \tilde\Pi^2(kL_\text{f},z),
\label{e:NoDimSpecDen}
\end{equation}
is a dimensionless spectral density for the gravitational waves. 

The gravitational wave power spectrum at time $t$ can then be written
\begin{eqnarray}
\rho_\text{GW} &=& (\bar\epsilon+\bar p)^2 \fluidV^4 (tL_\text{f}) (8\pi G \tilde\Omega_\text{GW}),
\label{e:RhoGWdot}
\end{eqnarray}
where 
\begin{equation}
\tilde\Omega_\text{GW} = \int \frac{dk}{k}\frac{(kL_\text{f})^3}{2\pi^2} \tilde P_{\text{GW}}(kL_\text{f})
\end{equation}
is a dimensionless number.  
We see that the gravitational wave power spectrum grows linearly with time, for as long as the velocity perturbations are active, with a slope which depends on the square of the enthalpy density, the fourth power of the mean square fluid velocity, the fluid length scale, and a dimensionless number describing the fluid flow $\tilde\Omega_\text{GW}$.

In principle, the value of $\tilde\Omega_\text{GW}$ depends on the parameters of the phase transition in dimensionless combinations, which we can expect to include the bubble wall speed $\vw$ and the latent heat relative to the total energy density $\strengthPar{T}$.  
In Fig.\ 2 (bottom) of \cite{Hindmarsh:2013xza}, we plotted $\rho_\text{GW}/\left[(\bar\epsilon+\bar p)^2 \fluidV^4 L_\text{f}\right]$ against time. 
Noting that we have $G=1$, the slope of the graph is $8\pi \tilde\Omega_\text{GW}$. 
We found that $\tilde\Omega_\text{GW}$ was approximately constant, varying by no more than a factor 2, when 
we took the fluid scale $L_\text{f}$ to be the mean bubble separation $\Rbc = \sqrt[3]{V/N_\text{b}}$. 
Hence most of the dependence of the gravitational radiation energy density on the phase transition parameters is accounted for by the explicit factors in Eq.~(\ref{e:RhoGWdot}).


\subsection{Integral scale}
\label{ss:IntSca}

The question of which scale to take for $L_\text{f}$ affects the value of $\tilde\Omega_\text{GW}$, and hence its variation between simulations.  
As mentioned above, in Ref.~\cite{Hindmarsh:2013xza} we took the scale to be $\Rbc$, the average bubble separation.
However, one could equally estimate the length scale from the velocity field itself, and to this end we can use the following quantity  (sometimes referred to as the integral scale) 
\begin{equation}
\xi_\text{f} =   \frac{1}{\vev{V^2}}\int \frac{d^3k}{(2\pi)^3} |k|^{-1} \SpecDen{V}(k),
\end{equation}
where $\vev{V^2}$ is the RMS velocity.
We will see that when the scale $L_\text{f}$ is chosen to be the integral scale, the variation in the parameter $\tilde\Omega_\text{GW}$ is reduced to about 10\%.  
This emergence of $\tilde\Omega_\text{GW}$ as a quasi-universal constant for first order phase transitions with $\strengthPar{\TN} \lesssim 0.1$ is an important result.

One can also define an integral scale $\xi_\text{GW}$ for the gravitational wave energy density from its spectral density $\tilde P_{\text{GW}}$.  We will also confirm that the integral scale of the gravitational radiation is related to the integral scale of the velocity field, as one would expect.

\begin{widetext}
\subsection{Dimensionless GW power spectrum parameter $\tilde\Omega_\text{GW}$ }

It is often useful to express the gravitational wave power spectrum as a fraction of the critical density, $\rho_\text{c} = 3H^2/8\pi G$. Hence we are led to consider 
a dimensionless gravitational wave power spectrum 
\begin{equation}
\frac{d \Omega_\text{GW}(k,t)}{d \ln(k)} =  \left[16\pi G(\bar\epsilon+\bar p)\fluidV^2\right]^2   \frac{t L_\text{f}}{\Hc^2} \frac{(kL_\text{f})^3}{24\pi^2} \tilde P_{\text{GW}}(kL_\text{f}),
\end{equation}
where $\Hc$ is the Hubble parameter at the time the bubbles are nucleated. Noting that the critical density is the energy density $\bar\epsilon$, and denoting the lifetime of the source by $\tau_\text{v}$, we find that the dimensionless gravitational wave power spectrum during the radiation era can be written   
\begin{eqnarray}
\label{e:GWPowSpe}
\frac{d \Omega_\text{GW}(k)}{d \ln(k)} &=&  3(1+w)^2 \fluidV^4  (\Hc \tau_\text{v}) (\HcL_\text{f}) \frac{(kL_\text{f})^3}{2\pi^2} \tilde P_{\text{GW}}(kL_\text{f}),
\end{eqnarray}
\end{widetext}
where $w = \bar{p}/\bar{\epsilon}$ is the equation of state parameter.
Integrating over wavenumber, we see that the total relative energy density is
\begin{equation}
\Omega_\text{GW} = 3(1+w)^2 \fluidV^4 (\Hc \tau_\text{v}) (\HcL_\text{f}) \tilde\Omega_\text{GW}.
\label{e:OmgwEqn}
\end{equation}

\subsection{Source lifetime}
\label{ss:LifSouWav}

It is clearly important for the calculation of the gravitational wave energy density to calculate the lifetime of the source, the shear stress caused by the sound waves. 
We show in Appendix \ref{s:GWExpUni} that in an expanding universe, the shear stresses decay and decorrelate in such a way to make $\tau_\text{v}$ precisely equal to the Hubble time.
The shear stresses also decay due to the viscosity of the fluid at a scale-dependent rate. We should therefore estimate on which scales viscous damping time is smaller than $\tau_\text{v}$.

For linear non-relativistic flows induced by sound waves (i.e.\ for velocity fields $V^i_\parallel$ which are purely longitudinal), viscosity adds a term of the form 
\begin{equation}
\left(\frac{4}{3}\etaS + \zetaB\right) \nabla^2 V^i_\parallel
\end{equation}
to the left hand side of Eq.\ (\ref{eqZ}), where $\etaS$ is the shear viscosity and $\zetaB$ is the bulk viscosity.  For a plasma of relativistic particles in a gauge theory, the bulk viscosity is negligible compared to the shear viscosity~\cite{Arnold:2006fz}, and the shear viscosity can be estimated as  
\begin{equation}
\etaS \sim T^3/e^4 \ln(1/e),
\end{equation}
where $e$ is the electromagnetic gauge coupling~\cite{Arnold:2000dr}.
Hence velocity perturbations of wavenumber $k$ are damped as $\exp(-4 \etaS k^2 t/3)$, and the lifetime due to viscous damping of sound waves with wavelength $R$ is  
\begin{equation}
\Tvisc(R) \sim R^2\epsilon/\etaS \sim e^4 \ln(1/e)R^2 T.
\end{equation}
Therefore, at a transition with temperature just below the critical temperature $\Tc$, the viscous damping lifetime exceeds the Hubble time $\Hc^{-1}$ for all scales 
\begin{equation}
R \gg \frac{\vw}{\Hc}\left(\frac{\sqrt{a}\Tc}{m_{\text{P}} e^4}\right) \sim 10^{-11} \frac{\vw}{\Hc}\left(\frac{\Tc}{100\; \text{GeV}}\right),
\end{equation}
where we have neglected the logarithm of the gauge coupling.

We will see in the next section that the scale of the fluid perturbations is set by the average separation of the nucleating bubbles $\Rbc$, and that the bubble separation at an electroweak-scale phase transition with any interesting degree of supercooling will satisfy this inequality. 

Hence for a first order transition at the electroweak scale -- or even a few orders of magnitude above -- the lifetime of the 
source of the gravitational waves
 is the Hubble time,
\begin{equation}
\tau_\text{v} = \Hc^{-1} \ll \Tvisc(\Rbc).
\end{equation}

\subsection{Comparison to envelope approximation}
\label{ss:ComEnvApp}
In the envelope approximation, the relative energy density in gravitational waves is given by~\cite{Kamionkowski:1993fg,Huber:2008hg}
\begin{equation}
\label{e:EnvAppFor}
\Omega_\text{GW}^\text{ea} \simeq \frac{0.11 \vw^3}{0.42+\vw^2} \left( \frac{\Hc}{\be}\right)^2 \frac{\kappa^2 \al^2}{(\al+1)^2},
\end{equation}
where $\al$ is the ratio between the ``vacuum'' energy (defined below) and the radiation energy density in the symmetric phase, $\kappa$ is the efficiency with which vacuum energy is converted to kinetic energy, and $\be$ is the nucleation rate parameter also defined above.

The vacuum energy $V_0$ is defined in Ref.~\cite{Kamionkowski:1993fg} from the trace anomaly, 
\begin{equation}
\theta = \epsilon - 3p,
\end{equation} 
as a quarter of the difference between the symmetric and broken phases: 
\begin{equation}
V_0 = \frac14(\theta_\text{s} - \theta_\text{b}).
\end{equation}
In our convention, the trace anomaly vanishes in the symmetric phase, and in the broken phase is
\begin{equation}
\label{e:VacEne}
\theta_{\text{b}} =  -T \frac{d}{dT}V(\phi_\text{b},T) + 4V(\phi_\text{b},T),
\end{equation}
where $\phi_\text{b}$ is the value of $\phi$ in equilibrium in the broken phase at temperature $T$. 
In the conventions of \cite{Kamionkowski:1993fg}, the trace anomaly vanishes in the broken phase, and in the symmetric phase is equal and opposite to (\ref{e:VacEne}). 
Hence for our thermal potential (\ref{e:ScaPot}) the parameter $\al$ is 
\begin{equation}
\al = \frac{V_0}{3aT^4} = \frac{1}{3aT^4}\left( \frac{1}{4}T \frac{d}{dT}V(\phi_\text{b},T) - V(\phi_\text{b},T) \right).
\label{e:StrParB}
\end{equation}
The efficiency parameter is defined from the average fluid kinetic energy density (\ref{e:fluidVdef}) as 
\begin{equation}
\kappa = \frac{1}{V_0} \frac{1}{{\mathcal V}} \int d^3 x \tau^\text{f}_{ii}.
\end{equation}
Therefore
\begin{equation}
\label{e:KappaEquiv}
(1+w) \fluidV^2 = \frac{\kappa\al}{1+\al}. 
\end{equation}
The factor of $1+\al$ in the denominator of the right hand side comes from the fact that we are dividing by the average total energy density in the symmetric phase, which is 
$3aT^4 + V_0$ in the conventions of \cite{Kamionkowski:1993fg}. 
 
Note that $\kappa\al$ is conventionally estimated analytically from the radial fluid velocity around an isolated expanding bubble
$v(r,t)$, where $r$ is the distance from the centre of the bubble, and $t$ is the time since nucleation~\cite{Kamionkowski:1993fg,Espinosa:2010hh}.
At large times, the radial fluid velocity is a function of a scaling variable $\xi = r/t$, rather than $r$ and $t$ separately. 
The ratio of the kinetic energy density to the total energy density can then be estimated as  
\begin{equation}
\kappa\al = \frac{3}{\vw^3 \epsilon} \int d\xi \xi^2 (\epsilon + p) W^2 v^2(\xi).
\end{equation}
We will compare this estimate to the numerically obtained $(1+w) \fluidV^2$ in the results section, finding good agreement.
 
In order to compare our expression for the gravitational wave energy density (\ref{e:OmgwEqn}) with the envelope approximation formula (\ref{e:EnvAppFor}), we estimate the fluid flow scale $L_\text{f}$ as the bubble separation scale $\Rbc$, which in turn is related to the nucleation rate parameter by Eq.\ (\ref{e:RbcBet}).

Hence the ratio between the gravitational wave energy density generated acoustically and in the envelope approximation is\footnote{Note that for a deflagration, if tunnelling is suppressed behind the shock wave, the ratio is boosted by a factor $\sim \cs/\vw$ -- see the discussion after Eq.\ (\ref{e:RbcBet}).}
\begin{equation}
\frac{\Omega_\text{GW}}{\Omega_\text{GW}^\text{ea}} \simeq \frac{3(8\pi)^{\frac13}\tilde\Omega_\text{GW}}{0.11\vw^2(0.42+\vw^2)}\left({\be}\tau_\text{v}\right).
\label{e:OmGWRat}
\end{equation}
Given that the ratio (\ref{e:OmGWRat}) is smallest for $\vw=1$, and that the lifetime of the sound waves is approximately $\Hc^{-1}$ (see Section \ref{ss:LifSouWav}), we can estimate that 
\begin{equation}
\frac{\Omega_\text{GW}}{\Omega_\text{GW}^\text{ea}} \gtrsim 60 \tilde\Omega_\text{GW} \frac{\be}{\Hc}.
\end{equation}
We will see from our numerical simulations that $\tilde\Omega_\text{GW} \sim 0.04$.
The ratio ${\be}/{\Hc}$ was discussed in Section~\ref{ss:BubNuc}, and shown to be at least $\mathrm{O}(10^2)$, and possibly significantly greater if there is only small supercooling.
We conclude that the energy density in acoustically generated gravitational waves is at least two orders of magnitude greater than the envelope approximation suggests.


\section{Numerical simulations}
\label{sect:numerical}

\subsection{Methods}

Our numerical methods are a development of those first used in this context to study the case of isolated bubbles in Ref.~\cite{KurkiSuonio:1995vy}. In that paper a spherically symmetric bubble was assumed. Here we extend those simulations to a full 3+1-dimensional simulation volume. In addition, we couple the linearised stress-energy tensor to perturbations around a flat metric, to measure the gravitational wave power produced by the simulation.

\subsubsection{Coupled field-fluid system}

\begin{figure}[t]
\begin{centering}
\includegraphics{space.eps}\\
\vspace{0.5cm}
\includegraphics{timespace.eps}\\
\caption{\label{fig:layout} Layout of quantities simulated. The positions of quantities related to simulating an ideal relativistic fluid are standard~\cite{WilsonMatthews}. Because the field and fluid are coupled together, it is important that the scalar field $\phi$ and its conjugate momentum $\pi$ are correctly centred. We take $\phi$ to reside in zones (like pressure, temperature, etc.), so that no centring is required to compute, for example, the equation of state.}
\end{centering}
\end{figure}


The coupled hydro-scalar equations, outlined above, can be treated quite easily using standard numerical techniques. The scalar field is evolved with the leapfrog (Verlet) algorithm, while standard operator splitting methods are used for the fluid~\cite{WilsonMatthews}. These are equivalent to numerically integrating the equations of motion given above.

Although the full details of how to implement relativistic hydrodynamics is beyond the scope of this paper, it is instructive to consider how the quantities are laid out on the lattice both in the spatial and temporal directions (see Fig.~\ref{fig:layout}). Furthermore, for good energy conservation it is essential that the discretised version of the damping term couple the field and fluid quantities at equal times during the simulation.

We have tested the results of our simulations against changing timestep (as well as the lattice spacing); see the following section.

As our simulations do not run for sufficiently long to develop strong shocks (indeed, we choose our lattice spacing parameters such that the fluid velocity profile is always resolved by several $\delta x$), the simulations presented in this paper do not involve any artificial viscosity. The importance of an artificial viscosity term was previously studied using 1+1-dimensional simulations of two colliding bubbles.

\subsubsection{Metric perturbations}

Our principal observables are the energy density and power spectrum of the gravitational waves. 
The goal of our simulations is to compute  the power per unit logarithmic frequency interval in gravitational waves $d\rho_\text{GW}(k)/d\ln k$, 
and the total energy density $\rho_\text{GW}$.

Perturbations of the metric are sourced by transverse-traceless part of the stress-energy tensor $\Pi_{ij}$
\begin{equation}
\ddot{h}_{ij} - \nabla^2 h_{ij} = 16 \pi G \Pi_{ij}.
\end{equation}
Obtaining $\Pi_{ij}$ from $T_{\mu\nu}$ involves a projection in momentum space. Therefore, evolving $h_{ij}$ (whether in momentum space or position space)
would involve Fourier transforms at each timestep.
As we go to large volumes, the execution time of fast Fourier transforms (FFTs) scales as $\mathrm{O}(N \log N)$, while there are few optimised FFT codes that offer domain decomposition in more than one dimension. 
It is therefore vital that the number of steps requiring Fourier transforms be minimised, to yield a scalable simulation.

Our approach is to evolve the unprojected equation of motion in real space~\cite{GarciaBellido:2007af}
\begin{equation}
\label{eq:unprojected}
\ddot{u}_{ij} - \nabla^2 u_{ij} = 16 \pi G (\tau^\phi_{ij} + \tau^\text{f}_{ij}),
\end{equation}
where $u_{ij}$ is an auxiliary tensor and 
the sources are defined in Eq.~(\ref{e:TauDef}).
Only when we wish to recover the metric perturbations $h_{ij}$ do we Fourier transform $u_{ij}$ and project out the transverse-traceless components through
\begin{equation}
h_{ij} (\mathbf{k}) = \lambda_{ij,lm} (\hat{\mathbf{k}}) u_{lm} (t,\mathbf{k}),
\end{equation}
where
the projector is defined in Eq.~(\ref{e:ProDef}).
We evolve Eq.~(\ref{eq:unprojected}) using a leapfrog algorithm in a similar manner to the scalar field.

Note that we choose the units of the code such that the critical temperature $\Tc = 1$ and the gravitational constant $G=1$.


\subsection{Tests}

Our basic tests principally involve varying the lattice spacing and timestep independently, on simulations of a single bubble colliding with itself in a small periodic box. These allow us to test that the simulations perform accurately between length scales $1/\Rbc$ and $1/\ell$. Longer distances do not need to be tested, and in any case $\Rbc$ is set by the box size $L$ in these tests.

\subsubsection{Changing the lattice spacing}

We performed tests on the effect of changing the lattice spacing using the self-collision of a single bubble in a cubic box, in relatively modest volumes (with parameters given in the following section). We considered $\delta x = 0.5/\Tc$, $\delta x = 1/\Tc$, $\delta x = 2/\Tc$ and $\delta x = 4/\Tc$. We notice no significant difference between these choices until $\delta x = 4/\Tc$.

\begin{figure}[t]
\begin{centering}
\includegraphics[width=0.4\textwidth,clip=true]{singlebubble.eps}\\
\caption{\label{fig:singlebubble} Single bubble test simulation, with the correlation length $\ell$ shown as an indication of the wall width. Only the fluid source is shown here; discretisation errors for the field source are the same or smaller. There is good agreement between $1/L$ and $1/\ell$, as desired.}
\end{centering}
\end{figure}

It is worth mentioning that, even for an isolated bubble which would (in continuum) have vanishing quadrupole moment and hence not source gravitational waves, the lattice discretisation breaks the spherical symmetry and results in a small amount of gravitational wave production. 
This power goes to zero as $(\delta x)^4$ for both the field and the fluid sources. After collision, however, agreement is very good with relative differences of at most $7\%$ for $k \lesssim 1/\ell$ between $\delta x=1/\Tc$ and $\delta x=2/\Tc$ for the fluid source. Furthermore, at higher momenta there are only $\mathrm{O}(1)$ differences between these two choices, consistent across seven orders of magnitude. This is surprisingly good given the relatively coarse wall width and the complicated microphysics. Similarly, discrepancies between $\delta x = 0.5/\Tc$ and $\delta x = 1/\Tc$ at late times were at worst $2 \%$ for $k \lesssim 1/\ell$; see Figure~\ref{fig:singlebubble}. Discretisation errors were always less severe for the field source than for the fluid source.

In summary, we note no significant sensitivity to lattice spacing so long as it is kept well below the scalar field wall width. 

While our previous work used simulations with $\delta x = 1/\Tc$, we use a lattice spacing of $\delta x = 2/\Tc$ in the present paper. The inferred discrepancies are demonstrably smaller than $10\%$, and the doubling of the accessible dynamic range that this allows is very useful.

\subsubsection{Changing the timestep}

With $\delta x=2/\Tc$ having been chosen, we varied the timestep to explore the effect of inaccuracies in our evolution algorithm. There is agreement at the $1\%$ level or better for all $k \lesssim 1/\ell$ and $5\%$ or better up to $k \sim 0.5$ (all points plotted on Fig.~\ref{fig:singlebubble}) as we varied $\delta t$ between $0.2/\Tc$, $0.1/\Tc$ and $0.05/\Tc$, in the same single-bubble tests for the fluid power spectrum described above. We use $\delta t=0.1$ for the remainder of the paper, although we could probably have achieved acceptable results with $\delta t=0.2$.

In the present paper our simulation durations are typically the same order of magnitude as one light-crossing time, and rather less than one sound-crossing time. This means, in particular, that the production of gravitational radiation by acoustic waves (or by scalar radiation, which is in any case heavily damped) is not likely to be affected by signals propagating around-the-lattice.

\subsection{Parameter choices}

We use the same parameters for the potential as in our previous paper. The exact values of these parameters are not particularly important: it is the latent heat and the wall velocity which mainly determine the gravitational wave power spectrum. Our aim in the present paper is to develop the ideas underlying our previous paper as well as the spherical studies carried out earlier, and so we work with the same parameters as before.

No attempt is made in the present paper to look at strong fluid flows or fast `runaway' bubble walls. We leave these harder topics for future work and instead seek to comprehensively explain the generation of gravitational waves by more gentle phase transitions ($\alpha_{T_\mathrm{N}} \lesssim 0.1$).

We discuss bubble nucleation further in the next section but note that we nucleate all of our bubbles simultaneously in the present work.

Given $\delta x =2/\Tc$ and a simulation size of $2400^3$ points, our physical simulation volume is $(4800/\Tc)^3$ for all the results presented in this paper.

\begin{table}

\begin{tabular}{l | c | c | c |}
 & Weak & Weak (scaled)  & Intermediate \\
\hline
$T_0/\Tc$ & $1/\sqrt{2}$ & $1/\sqrt{2}$ & $1/\sqrt{2}$ \\
$\gamma$ & $1/18$ & $4/18$ & $2/18$ \\
$A$ & $\sqrt{10}/72$ & $\sqrt{10}/9$ & $\sqrt{10}/72$ \\
$\lambda$ & $10/648$ & $160/648$ & $5/648$ \\
${\cal L}/T_c^4$      & $9/40$ & $9/40$ & $9/5$ \\
$\sigma/T_c^3$ & $1/10$ & $1/20$ & $4\sqrt{2}/10$ \\
$\ell \Tc$ & $6$ & $3$ & $6/\sqrt{2}$ \\
$\TN/\Tc$ & $0.86$ & $0.86$ & $0.8$\\
$\al_{\TN}$ & $0.010$ & $0.010$ & $0.084$ \\
$\al$ & $0.0046$ & $0.0046$ & $0.050$ \\
$\Rc \Tc$ & $ 16$ & $8.1$ & $8.6$ \\
\hline

\end{tabular}

\caption{Scalar potential parameters (\ref{e:ScaPot}), nucleation temperature $\TN$, 
phase transition parameters (\ref{e:PhaTraPar}), transition strength parameters 
(\ref{e:StrParA}) and (\ref{e:StrParB}), and critical bubble radii (\ref{eq:critbub})  for our simulations.}

\label{t:SimParsPot} 
\end{table}


\subsection{Initial conditions}

\label{sec:initconds}

At the start of our simulation, we nucleate a controllable number of bubbles, which was usually $\mathrm{O}(1000)$ (yielding bubbles of average collision radius slightly larger than in Ref.~\cite{Hindmarsh:2013xza}), but was as small as 37 or in one case as large as 32558. These have a Gaussian scalar field profile. This profile is initially at rest, meaning that the conjugate momentum, and also the fluid velocity are zero in the vicinity of the bubble. We ensure that all the initially nucleated bubbles are well separated at the start of the simulation. For runs with the same number of bubbles but different wall velocities, all bubbles are nucleated at the same positions, but from testing we found that even 37 bubbles was enough to remove any discernible dependence on the initial bubble positions.

The critical bubble radius can be computed from the surface tension $\sigma$ and the difference in potential energy at $\TN$ from the thin-wall formula (noting that the potential energy in the symmetric phase is zero)
\begin{equation}
\label{eq:critbub}
\Rc = \frac{2\sigma}{- V(\phi_b,\TN)}.
\end{equation}
Values of $\Rc$ for our simulations are shown in Table~\ref{t:SimParsPot}. Rather than find the critical bubble profile exactly, we use a spherically symmetric Gaussian field profile
\begin{equation}
\phi(r) = \phi_\text{b} \exp(-r^2/2\Rc).
\end{equation}
This is rather broad, and therefore sufficiently large compared to the true critical bubble profile to ensure that the bubbles reliably expand despite lattice effects.

The bubbles are sufficiently large that they immediately start growing, driven by the pressure difference between the interior and the exterior.  The  scalar field quickly settles into a kink-like configuration, interpolating between the metastable and stable minima over a distance of order $\ell$, the correlation length of the scalar field (see Table~\ref{t:SimParsPot} for the values the correlation length takes).  For the scalar field dynamics to be valid we must have a lattice spacing that resolves the wall width (see previous section), which places an upper limit on the physical simulation volume possible for a given amount of computer memory.

In this paper, the bubbles are nucleated simultaneously. 
Nucleating at a single time helps to ensure clear scale separation in the limited dynamic range available to our numerical simulations, although it does produce oscillatory patterns in the resulting power spectrum (we cover the case of unequal nucleation times in Appendix~\ref{a:BubNucTim}). We could in principle recover the power spectrum produced by bubbles of all different sizes by a linear 
superposition of the resulting power spectra, weighted by the bubble size distribution.

Once nucleated, the bubbles grow, and the fluid approaches a characteristic radial velocity distribution, which is a function of $\xi=r/t$, where $r$ is the distance from the centre of the bubble, and $t$ is the time since nucleation (see Fig.~\ref{fig:profilecomparison}). The form of this function depends  on the bubble wall velocity~\cite{KurkiSuonio:1995vy}, and we will refer to it as the scaling profile. The rate of approach to the fluid scaling profile is generally much slower than the relaxation of the scalar field. 
In a true electroweak phase transition the bubble size at collision is many orders of magnitude larger than the bubble size at nucleation, giving a lot of time for the radial velocity distribution to reach its asymptotic profile. 

In our numerical simulations, the ratio of the bubble size at collision, $\Rbc$ to the bubble size at nucleation ($\approx \Rc$) is at most 90 for $N_\text{b}=37$, $\ell=16$ and as small as 9.4 for $N_\text{b}=32558$, $\ell=16$ (our simulation parameters are outlined in Tables \ref{t:SimParsPot} and \ref{t:SimParsRuns}).
We should therefore be alert to the fact that the fluid has definitely not settled down to its final scaling profile in a collision.  One can test for the effect of a non-scaling fluid profile by repeating simulations with fewer bubbles, so that there is a longer time before collision. We have carried out simulations such that $\Rb$ varies by around a factor of three in the two sets of simulations for which we present plots, and by as much as a factor of ten in our full set of simulations for this paper.

\begin{figure}
\begin{centering}
\includegraphics[width=0.4\textwidth,clip=true]{profilecomp.eps}

\caption{\label{fig:profilecomparison} Comparison of radial fluid velocity profiles for simulations at approximate collision times for $N_\text{b}=1000$ ($t=500/\Tc$; gray), $N_\text{b}=37$ ($t=1000/\Tc$; black) and at late times (dashed red). }

\end{centering}
\end{figure}


\begin{table}
\begin{tabular}{cccrcccc}
    Type      & $\eta/T_c$ & $\vw$ & $N_\text{b}$ & $\fluidVmax$ & $\fluidVmaxperp$ & $\xi_\text{f, end} T_c$ & $8\pi\tilde\Omega_\text{GW}$ \\
\hline
  Weak        &  0.06  &  0.83 &   988 &      0.0052  &         0.00037  &     351              &     0.88       \\
              &        &       &   125 &      0.0052  &         0.00028  &     649              &     0.84       \\
\cline{2-8}
              &  0.1   &  0.68 &   988 &      0.0084  &         0.00036  &     244              &     0.73       \\
              &        &       &   125 &      0.0082  &         0.00026  &     451              &     0.71       \\
              &        &       &    37 &      0.0080  &         0.00021  &     644              &     0.60       \\
\cline{2-8}
              &  0.121 &  0.59 &   988 &      0.0116  &         0.00052  &     182              &     0.69       \\
\cline{2-8}
              &  0.15  &  0.54 &   988 &      0.0102  &         0.00037  &     230              &     0.54       \\
              &        &       &    37 &      0.0120  &         0.00025  &     428              &     0.80       \\
\cline{2-8}
              &  0.2   &  0.44 & 32558 &      0.0059  &         0.00047  &     136              &     0.97       \\
              &        &       &   988 &      0.0073  &         0.00031  &     368              &     0.70       \\
              &        &       &   125 &      0.0075  &         0.00023  &     613              &     0.86       \\
              &        &       &    37 &      0.0078  &         0.00019  &     942              &     0.70       \\
\cline{2-8}
              &  0.4   &  0.24 &   988 &      0.0036  &         0.00049  &     756              &     0.86       \\
\hline
Wk. (sc.)     &  0.4   &  0.44 &   988 &      0.0075  &         0.00029  &     365              &     0.81       \\
\hline
Interm.       &  0.4   &  0.44 &   988 &      0.0595  &         0.00328  &     485              &     1.04       \\
\hline
\end{tabular}
\caption{\label{t:SimParsRuns} Simulation parameters $\eta$
  (field-fluid coupling), $N_\text{b}$ (number of bubbles nucleated), with
  the resulting bubble wall speed $\vw$, the maximum fluid RMS
  velocity $\fluidVmax$, the maximum contribution of transverse fluid
  motion $\fluidVmaxperp$, the integral scale of the fluid
  $\xi_\text{f, end}$, and the scaled slope parameter for the
  growth of the gravitational wave energy density $\tilde\Omega_\text{GW}$. The
  potential parameters and derived quantities for each type --
  ``weak'', ``weak scaled'' and ``intermediate'' -- are given in
  Table~\ref{t:SimParsPot}.}
\end{table}


\subsection{Scaling}

A cosmological first order phase transition is a multiscale problem, with length scales varying from the microscopic ($1/T$, bubble wall thickness) up to the Hubble scale, a range spanning 
17 orders of magnitude at the electroweak scale.
The typical bubble sizes  at collision time are somewhere between these scales, depending on the metastability of the high temperature phase.  It is of course impossible to include all of these scales in a single numerical simulation, where scale hierarchies only of order $10^2$ are achievable.  To obtain a stable numerical description of the bubble wall, the wall thickness has to span a few lattice units (denoted by $\delta x$).  In order to have collisions within the simulation volume, this restricts the bubble separation to be unphysically small.

However, this restriction can be relaxed, at least partly:
we expect that the dynamics of the bubble growth, collisions and the subsequent generation of gravitational waves are mostly determined by the ``bulk'' thermodynamics ($\epsilon$, $p$, latent heat ${\cal L}$) and the friction parameter $\eta$, but not by microscopic details of the bubble wall (surface tension $\sigma$, wall thickness $\ell$).  Dimensionally, it is clear that the contribution from quantities proportional to bubble volume (e.g. latent heat) will dominate over quantities proportional to the area of the bubbles when the bubble radius is large enough.

This motivates us to search for a way to modify the equations of motion (\ref{eqfield})--(\ref{eqZ}) so that we could simulate bubbles which are significantly larger than the microscopic length scale, while preserving the bulk thermodynamics of bubble expansion while possibly sacrificing the properties of the bubble wall.  Indeed, this can be achieved with the following simple rescaling of the parameters and fields:
\begin{equation}
  \begin{array}{l}
  \gamma \rightarrow r^2 \gamma, ~~
  A \rightarrow r^3 A, ~~
  \lambda \rightarrow r^4 \lambda, ~~
  \eta \rightarrow r \eta,  \\
  \phi(x) \rightarrow r^{-1}\phi(rx), ~~
   V^i(x) \rightarrow V^i(rx),\\
   E(x) \rightarrow E(rx),~~~~~
   Z_i(x) \rightarrow Z_i(rx),
  \end{array}
  \label{scaling}
\end{equation}
Here $r$ is a dimensionless scaling factor, and $x=(\textbf{x},t)$.
Clearly, the equations of motion (\ref{eqfield})--(\ref{eqZ}) remain valid.  The crucial feature of the scaling is that the potential remains constant, $V(\phi,T) = V_{\rm scaled}(r^{-1}\phi,T)$, indicating that the bulk quantities $T_c$, ${\cal L}$ and also $\epsilon(T,\phi)$ and $p(T,\phi)$ remain invariant, as desired.  However, the surface tension and wall thickness scale as $\sigma\rightarrow r^{-1}\sigma$ and $\ell\rightarrow r^{-1}\ell$.  In effect the scaling stretches the field configuration by a factor of $r^{-1}$ in spatial and temporal directions.

We note that in spite of the non-trivial scaling of the friction parameter $\eta$, the total frictional force imparted on the moving bubble wall does not change: it is obtained by integrating the $\eta$-terms in Eqs.~(\ref{eqfield})--(\ref{eqZ}) over the bubble wall thickness, which is scaled by a factor of $r^{-1}$.

What does the rescaling gain us?   It is straightforward to see that the lattice implementation of the equations of motion (\ref{eqfield})--(\ref{eqZ}) do not change (in lattice units) under scaling (\ref{scaling}), provided that the lattice spacing is also scaled as $\delta x \rightarrow r^{-1} \, \delta x$.  This implies that a single lattice simulation exactly corresponds to a whole family of results, given by the scaling with $r$.  All of them have the same bulk thermodynamical properties.  Thus, provided that the detailed bubble wall properties are not important for bubble collisions and gravitational wave generation, we can take a simulation run where bubbles have been nucleated at specific locations, and rescale it to 
the desired physical bubble separation scale.

We can test the assumption that the surface properties are not important by comparing results from simulations which differ only in surface tension and wall thickness.  This can be achieved by applying the scaling (\ref{scaling}) to the parameters of the theory, but leaving the lattice spacing constant.  A set of parameters for unscaled ($r=1$) and scaled ($r=2$) runs are shown in Table \ref{t:SimParsPot}.  The surface tension and the bubble wall thickness have been halved in the scaled simulation.  The results of the test are shown in Table \ref{t:SimParsRuns}; here the scaled run is done with $\eta/\Tc=0.4$ and 988 bubbles, which can be compared with unscaled $\eta/\Tc=0.2$, 988 bubble results.  In both simulations the bubbles are nucleated at identical times and locations.  The numerical results match well within uncertainties of the measurements, supporting our assumption that the surface properties of the scalar field profile are unimportant.

\section{Numerical results}
\label{s:NumRes}

\begin{figure*}[t]


\begin{centering}
\includegraphics[height=0.29\textwidth]{slice0500.eps} \hfill
\includegraphics[height=0.29\textwidth]{slice1000.eps} \hfill
\includegraphics[height=0.29\textwidth]{slice1500.eps} \hfill
\includegraphics[height=0.29\textwidth]{legend.eps}
\caption{\label{fig:slice} Slices of fluid kinetic energy density $E/\Tc^4$ at $t=500\, T_\mathrm{c}^{-1}$, $t=1000 \,T_\mathrm{c}^{-1}$ and $t=1500\, T_\mathrm{c}^{-1}$ respectively, for the $\eta/T_c=0.15$, $N_\text{b}=988$ simulation.}
\end{centering}
\end{figure*}



In this section we present the main results from a campaign of numerical simulations, whose parameters are given in Table \ref{t:SimParsRuns}.  
As mentioned before, our simulations were carried out in a volume $(4800/\Tc)^3$.
A set of representative slices through the simulation are shown in Fig.~\ref{fig:slice}.
Our main results are derived from the set of simulations with latent heat to thermal energy ratio $\al_{\TN}\simeq 0.01$, which we characterise as a ``weak'' transition.  
Values of the friction parameter $\eta$ were chosen to give bubble growth proceeding by both detonation and deflagration, as well as one simulation tuned to the Jouguet case, where the bubble wall moves at the speed of sound. We have one ``intermediate'' strength transition, with $\strengthPar{\TN} \simeq 0.1$ where $\eta$ is chosen to give the same wall speed as the weak transition with $\eta/\Tc = 0.2$.  The ``weak scaled'' transition (employing the scaling of the previous section) is discussed later.

When plotting graphs, we focus on three representative cases, where the field-fluid coupling is $\eta/T_c=0.1$, $\eta/T_c=0.15$, and $\eta/T_c=0.2$,   
and the bubble wall speed is supersonic ($\vw = 0.83$), just subsonic ($\vw = 0.54$), and subsonic ($\vw = 0.44$). 
A complete set of graphs can be found in the supplementary material~\cite{supplementary}.

Our understanding of the transition developed in Section \ref{s:TheGWGen} shows that the important quantities for the overall gravitational wave energy density are the RMS fluid velocity $\fluidV$ and the fluid velocity scale $L_\text{f}$, and 
that the gravitational wave power spectrum is only indirectly dependent on the strength of the transition and the parameters in the potential.  Indeed, the gravitational wave power spectrum should be the same for parameters which give the same $\al_{\TN}\simeq 0.01$ (keeping the wall velocity constant). We can use the ``weak scaled'' run of Table \ref{t:SimParsRuns} to test this statement, where $\al_{\TN}$ is constant but the scalar bubble wall thickness is halved.  We test the effect of the strength of the transition with the ``intermediate'' run of Table \ref{t:SimParsRuns}.

We track the progress of the transition through  
the time evolution of the two quantities $\fieldV$ and $\fluidV$ defined in Eqs.\ (\ref{e:fluidVdef}),~(\ref{e:fieldVdef}).
We recall that the squares of these quantities give an estimate of the size of the shear stresses of the field and the fluid relative to the background fluid enthalpy density, and that the $\fluidV$ tends to the RMS\ fluid velocity for $\fluidV \ll 1$. We also note that the fraction of the fluid velocity power coming from rotational modes, $\fluidVperp$, is very small leading us to conclude that rotational fluid modes are not important in this system; we discuss this in more detail in Appendix~\ref{a:TraVelNeg}.

\begin{figure}
\begin{centering}
\includegraphics[width=0.4\textwidth,clip=true]{dimless-1000.eps}
\includegraphics[width=0.4\textwidth,clip=true]{dimless-37.eps}

\caption{Root mean square fluid velocity $\fluidV$ and root mean square scalar gradients $\fieldV$ for $N_\text{b} = 988$  (top row) and $N_\text{b} = 37$   (bottom row). }
\label{fig:timevolV} 
\end{centering}
\end{figure}



We see from Fig.~\ref{fig:timevolV} that $\fieldV$ grows and decays with the total surface area of the bubbles of the new phase, while the mean fluid velocity grows with the volume of the bubbles, and then stays constant once the bubbles have merged\footnote{We have no explicit viscosity, and the slight decreasing trend in some measurements of $\fluidV$ arises from the well-known numerical viscosity of donor-cell advection, $\nu_\text{num} \simeq \fluidV \delta x$.}.
This allows us to identify distinct phases of the transition: the collision phase, where $\fieldV$ grows and decays; and the subsequent acoustic phase where $\fluidV$ is approximately constant, and $\fieldV$ vanishes. 



\subsection{Length scales}

The analysis of Section \ref{s:TheGWGen} shows that the length scale of the velocity flow is an important determinant of the gravitational wave power spectrum. In Fig.\ \ref{fig:integralscale} we show the integral scales for the velocity and the gravitational radiation for runs at $N_\text{b}=37$ and $N_\text{b}=988$. During the collision phase, the bubbles expand and overlap, and hence the scales of the velocity field and the resulting gravitational radiation grow linearly in time. The gravitational radiation length scale is 2-3 times that of the velocity field. The scale of the velocity field stops growing as the bubbles collide and the scalar field decays to the vacuum, and stays constant during the acoustic phase. The scale imprinted on the gravitational radiation during the acoustic phase is close to that of the velocity field.


\begin{figure}[t]
\begin{centering}
\includegraphics[width=0.4\textwidth,clip=true]{integralscale-1000.eps}
\includegraphics[width=0.4\textwidth,clip=true]{integralscale-37.eps}
\caption{\label{fig:integralscale} Plot of integral scales $\xi_\text{f}$ and $\xi_\text{gw}$ associated with the fluid and gravitational wave power spectrum for simulations of interest for $N_\text{b} = 988$ (top) and $N_\text{b} = 37$ (bottom). Note the different $y$-axis plotting scales.}
\end{centering}
\end{figure}


\subsection{Velocity profile}

Given the discussions on initial bubble sizes in Section~\ref{sec:initconds}, it is important to bear in mind that the bubbles in our simulations expand in size by a factor of only around 10-100,
which is many orders of magnitude less than in a real phase transition.  One practical effect is that the profile of the velocity field around the bubbles does not reach its asymptotic scaling form, which can be expressed in terms of the previously introduced ratio $\xi = r/t$. In Fig.\ \ref{fig:profilecomparison} we showed the velocity profiles for the weak transition at $\eta/T_c=0.1$, $\eta/T_c=0.15$, and $\eta/T_c=0.2$, after times $t = 500/\Tc$ and $t = 1000/\Tc$.  These are approximately when most bubble collisions are happening, in the $N_\text{b} = 1000$ and $N_\text{b} = 37$ runs respectively.

We see that, at collision, the velocity profiles are qualitatively similar to their asymptotic forms in amplitude and shape, but differ in detail. In particular, the peak velocities are lower. 
This is particularly noticeable at the earlier time. We would therefore expect the RMS velocities $\fluidV$ measured in the simulations to be underestimates.  As the gravitational wave power spectrum depends on the fourth power of $\fluidV$, this is a significant source of uncertainty in deriving accurate predictions for the gravitational wave power spectrum.


\begin{table}[h!]
\begin{tabular}{ccccc}
     $\eta/T_c$ & $\vw$ &  $\fluidV$ &  $\sqrt{\frac{3}{4}\kappa_\text{1d}\al}$  &  $\sqrt{\frac{3}{4}\kappa_\text{Esp}\al}$   \\
\hline
                0.06  &  0.83 &        0.0052     &      0.0056      &      0.0063      \\
                0.1   &  0.68 &         0.0084     &       0.0085     &       0.0121      \\
                0.121 &  0.59 &         0.0116     &      0.0146     &      0.0192     \\
                0.15  &  0.54 &         0.0102     &       0.0103     &   0.0100      \\
                  0.2    & 0.44      &        0.0073     &   0.0066     &      0.0065        \\
                0.4   &  0.24 &         0.0036     &       0.0033     &      0.0036     \\
\hline
\end{tabular}
\caption{\label{tab:tablethree} Simulation parameters $\eta$ (field-fluid coupling), with the resulting bubble wall speed $\vw$, fluid RMS velocity $\fluidV$, for weak transitions with $N_\text{b}=988$, and the equivalent quantity $\sqrt{4\kappa\al/3}$ appearing in the envelope approximation (see Eq.\ \ref{e:KappaEquiv}). The efficiency parameter $\kappa$ is estimated in two ways: 
$\kappa_\text{1d}\al$ is estimated from the numerical spherically symmetric 1D fluid profiles at $t=1000/\Tc$, while 
$\kappa_\text{Esp}$ comes from the function 
$\kappa(\vw,\alpha)$ given in the Appendix of Ref.~\cite{Espinosa:2010hh}, using 
$\vw$ extracted from spherical 1D simulations at $t=1000/\Tc$.
   }
\end{table}


These considerations are tested in Table \ref{tab:tablethree}, where we compare our mean square fluid velocity parameter $\fluidV$ with $\sqrt{\frac{3}{4}\kappa\al}$, which should be equal according to the discussion in Section \ref{ss:ComEnvApp}. In the table, the efficiency parameter $\kappa$ is estimated in two ways: $\kappa_\text{1d}\al$ is estimated from integrating the numerical 1D fluid profiles at $t=1000/\Tc$, while $\kappa_\text{Esp}$ comes from the function $\kappa(\vw,\alpha)$ given in the Appendix of Ref.\ \cite{Espinosa:2010hh}, using $\vw$ extracted from 1D simulations at $t=1000/\Tc$. As can be seen, $\fluidV$ from the 3D simulations compares reasonably well to its estimate extracted from the 1D numerical profiles around the time of bubble collision, while the theoretical values are somewhat higher. It is remarkable that such a simple model for the mean square velocity, which omits all details of the bubble collisions, does so well.


\subsection{Power spectra}

In Figs.\ \ref{f:VelPSMulti} and \ref{f:GWPSMulti} we show velocity and gravitational wave power spectra at various times through the simulations, for weak transitions with $\eta/T_c=0.1$, $\eta/T_c=0.15$, and $\eta/T_c=0.2$, where the bubble wall speed is supersonic ($\vw = 0.83$), just subsonic ($\vw = 0.54$), and subsonic ($\vw = 0.44$). The same potential and fluid-field parameters are run with $N_\text{b}=988$ and $N_\text{b}=37$ bubbles, to show the effect of allowing a greater time for the fluid velocity around the expanding bubbles to approach their scaling profiles.  The power spectra develop in characteristic ways in the different phases of the transition, and one can see that if the simulation is stopped too early, a misleading impression of the power spectrum will be obtained.





\begin{figure*}[tb]
\begin{centering}
\includegraphics[width=0.325\textwidth,clip=true]{velps-0.1-1000.eps}
\hfill
\includegraphics[width=0.325\textwidth,clip=true]{velps-0.15-1000.eps}
\hfill
\includegraphics[width=0.325\textwidth,clip=true]{velps-0.2-1000.eps}\\
\vspace{2ex}
\includegraphics[width=0.325\textwidth,clip=true]{velps-0.1-37.eps}
\hfill
\includegraphics[width=0.325\textwidth,clip=true]{velps-0.15-37.eps}
\hfill
\includegraphics[width=0.325\textwidth,clip=true]{velps-0.2-37.eps}\\

\caption{Velocity power spectra, for weak transitions, at $\eta/T_c = 0.1$, $0.15$ and $0.2$ ($\vw = 0.83$, $0.54$ and $0.44$) for $N_\text{b} = 988$ (top row) and $N_\text{b}=37$ (bottom row). The large oscillations are due to all the bubbles being nucleated at exactly the same time. As in Fig.~\ref{f:GWPSMulti}, we note that the scales are standardised for all the plots, but that the phase transition has not necessarily finished by $2500/Tc$, the time of the latest curve.}

\label{f:VelPSMulti}

\end{centering}
\end{figure*}


\begin{figure*}[tb]
\begin{centering}
\includegraphics[width=0.325\textwidth,clip=true]{gwps-0.1-1000.eps}
\hfill
\includegraphics[width=0.325\textwidth,clip=true]{gwps-0.15-1000.eps}
\hfill
\includegraphics[width=0.325\textwidth,clip=true]{gwps-0.2-1000.eps}\\
\vspace{2ex}
\includegraphics[width=0.325\textwidth,clip=true]{gwps-0.1-37.eps}
\hfill
\includegraphics[width=0.325\textwidth,clip=true]{gwps-0.15-37.eps}
\hfill
\includegraphics[width=0.325\textwidth,clip=true]{gwps-0.2-37.eps}\\

\caption{Gravitational wave power spectra, for weak transitions, at $\eta/T_c = 0.1$, $0.15$ and $0.2$ ($\vw = 0.83$, $0.54$ and $0.44$) for $N_\text{b} = 988$ (top row) and $N_\text{b}=37$ (bottom row). Note that the axes and time intervals are the same for all plots, which means that in some cases the latest ($2500/\Tc$) curve is from before the completion of the phase transition. }

\label{f:GWPSMulti}

\end{centering}
\end{figure*}


\subsubsection{Collision phase}

Looking first at the velocity power spectra, the most striking feature is their periodic modulation. This is not a physical feature, and is due to the bubbles being nucleated all at the same time. We have checked that spreading the nucleation times reduces this modulation, and it is not expected to be a feature of the velocity power spectrum of a realistic bubble nucleation distribution in the infinite volume limit. In Appendix \ref{a:BubNucTim} we show the effect of allowing nucleation over a time of about $200/\Tc$.

Once the fluid shells of the nearest pair of bubbles begin to overlap, gravitational waves are generated, at a scale controlled by the size of the bubbles. The overlap of the fluid shells is quickly followed by the collisions of the bubble walls, and gravitational radiation is generated by the scalar field as well. The bubbles continue to grow and to collide, and as a result the length scales of the velocity field and the gravitational radiation get larger (see Fig.~\ref{fig:integralscale}). This effect can be seen in the power spectra, where the curves show a peak moving up and to the left with time.

In our simulations there is generally more energy in the scalar field than in the fluid to begin with, and so the gravitational radiation from the scalar field dominates the early phases (see Fig.~\ref{f:FluVsAll}).  However, when scaled to a real deflagration or detonation in the early universe, most of the latent heat of the transition goes to the fluid, and the radiation from the scalar field can be neglected. It is only in the case of a runaway bubble wall that the scalar field takes most of the latent heat.  We discuss the scaling to real transitions in Section \ref{ss:Extrapolate}, and we plan to study runaway transitions elsewhere.

\begin{figure}
\begin{centering}
\includegraphics[width=0.4\textwidth,clip=true]{fieldfluidcomparison.eps}
\caption{Power spectra for $\eta/T_c=0.2$, comparing fluid-only (dashed) and total (solid) GW power at intervals of $500/\Tc$. The power laws visible in the `total GW power' case are dominated by the gradient energy of the scalar field. This source, however, is short-lived. We conjecture that it can be calculated by means of the envelope approximation.}
\label{f:FluVsAll}
\end{centering}
\end{figure}

However, it is interesting to study the difference between a fluid-only gravitational wave power spectrum, and one sourced by both fluid and field (Fig.~\ref{f:FluVsAll}).
There one sees evidence of a $k^{-1}$ power spectrum in arising from the scalar field during the collision phase (solid lines), which is later dominated by the gravitational waves from the fluid. As the scalar field energy density is confined to a thin shell, it is reasonable to suppose that its contribution can be adequately computed in the envelope approximation in the collision phase. We will investigate this conjecture elsewhere.


\subsubsection{Acoustic phase}

Eventually, the low-temperature phase spreads throughout the volume, the scalar field domain walls disappear, and fluid velocity perturbations are left behind.  We call this the acoustic phase of the transition, as the fluid perturbations are primarily compressive (longitudinal) modes (see Appendix \ref{a:TraVelNeg}). During the acoustic phase, the length scale of the fluid perturbations and the gravitational waves remains constant.

For the simulations with fewer bubbles ($N_\text{b}=37$), we see from the lower row of Fig.~\ref{f:VelPSMulti} that the envelope of the velocity power spectrum has an approximate power-law envelope beyond the peak. This power-law envelope is also visible at $N_\text{b} = 988$ at $\eta/\Tc=0.2$, where the bubble wall speed is lower. In both cases at $\eta/\Tc=0.2$, the bubbles expand longer before collision, and we expect the velocity field to be closer to the asymptotic form. We note that the power law is approximately $k^{-1}$ at $\eta/\Tc=0.2$, and appears steeper for lower couplings.  However, we are not confident that we have reached the asymptotic form for bubble wall speeds above $\vw = 0.44$.

At low wavenumbers, the velocity power spectrum behaves as a power of $k$, and arguments based on the analyticity properties of the Fourier transform of a longitudinal vector field in Ref.~\cite{Caprini:2007xq} show that it should go as $k^5$. This is just visible in the first few bins of the simulations with $N_\text{b}=988$. Larger simulations are required to properly check the long-wavelength behaviour.

We see that the gravitational wave power spectrum grows linearly with time in the acoustic phase, maintaining its shape, except at the lowest wavenumbers. A power-law behaviour can be seen emerging beyond the peak, especially in the simulations with $N_\text{b} = 37$. 
The power-law is approximately $k^{-3}$ for the weak deflagrations at $\eta/\Tc=0.2$ and $\eta/\Tc=0.4$ (the power spectra for the latter can be found in the supplementary material for this work) for both $N_\text{b} = 988$ and $N_\text{b} = 37$, which gives us confidence that we are close to the true power law. 
However, a power-law can be seen only for $N_\text{b} = 37$ for $\eta \le 0.15$.  Without further simulations at larger $\Rbc$ we cannot properly determine the long-wavelength behaviour of the growing acoustic phase power spectrum in these cases.


\begin{figure}[t]
\begin{centering}
\includegraphics[width=0.4\textwidth,clip=true]{samerate.eps}\\
\caption{\label{fig:timevolGW} Time series of $\rho_\text{GW} L_{\mathrm{f}}^{-1} [(\bar \epsilon+ \bar p)^{-2}\fluidV^{-4}]$, showing the evolution of the gravitational wave energy density relative to an estimate of the square of the final fluid shear stresses. We take the fluid length scale $L_\mathrm{f}$ to be the integral scale $\xi_\mathrm{f}$. Some oscillation about the constant curve is caused by long-wavelength sloshing of the fluid or the infrared behaviour of the gravitational wave power, discussed later, but the striking feature is the scalable linearity of the signal across a factor of three for $R_*$. Only fluid contributions to the gravitational wave power are included here. The early-times steep growth is best explained by the violent behaviour when the two shocks overlap. This phase is not well explained by our random velocity field model.
}
\end{centering}
\end{figure}

We argued earlier that the gravitational wave density parameter $\Omega_\text{GW} = \rho_\text{GW}/\bar\epsilon$ is proportional to $L_\text{f}$, the fluid velocity length scale,  and the square of the volume-averaged fluid energy density $(\bar\epsilon + \bar p)^2 {\overline U}_\mathrm{f}^4$. We plot the scaled gravitational wave energy density in Figure~\ref{fig:timevolGW}. This plot shows nicely parallel, linear growth of gravitational wave power when rescaled by these quantities at late times.  The coincidence of the slopes is greatly improved over the equivalent figure in Ref.~\cite{Hindmarsh:2013xza}, thanks to the larger simulation volumes, longer run times, and above all the replacement of the average bubble separation at collision by the fluid integral scale.

The improved coincidence of the slopes is one of the major results of the paper.  It establishes the existence of an $\mathrm{O}(1)$ parameter $8\pi\tilde\Omega_\text{GW}$ for a wide range of relevant transitions, and shows that the gravitational wave energy density from a phase transition can be understood in terms of simple features of the velocity field created by the dynamics of the bubble collision. 


\subsection{Extrapolating to a real phase transition}
\label{ss:Extrapolate}

Our simulations are necessarily limited in volume, duration, and resolution. We now discuss how they can be extrapolated to the real universe.  In particular we would like to extrapolate the gravitational radiation power spectrum, expressed as a fraction of the critical density. 

There are three physical length scales in the system: the average bubble separation $\Rbc$, the size of the initial bubble of the broken phase $\Rc$, 
and the bubble wall width $\ell$.  They are all set by the dimensional scale of the effective potential, which one can chose to be the critical temperature $\Tc$, and various combinations of the dimensionless couplings $\gamma$, $A$ and $\lambda$.  In a real transition, the average bubble separation is much larger than the wall width because of the exponential factor in the tunnelling rate, whose argument is set by the ratio of the energy of the critical bubble (see Ref.~\cite{Kirzhnits:1976ts}) to the critical temperature. This is generally a large number.

There are also two physical time scales to consider: the lifetime of the fluid flow $\tau_\text{v}$, and the duration of the phase transition, which is of order $\be^{-1}$, the inverse of the tunnelling rate parameter~(\ref{e:TunRat}).  The duration of the phase transition is also of order $\Rbc/\vw$, the time it takes for bubbles of average separation to collide.

Finally, there are also scales set by the background cosmology: the Hubble rate at the phase transition $\Hc$ (and the Hubble length), and the gravitational constant $G$. The Hubble rate, the gravitational constant and the critical temperature $\Tc$ are related via the Friedmann equation. Our simulations are performed in a Minkowski background, as the duration of the transition is assumed to be comparable to the Hubble time. Therefore $\Tc$ and $G$ can be chosen independently. The role of $G$ is purely to set the scale of of the gravitational perturbations. As mentioned earlier, we use units $\Tc = 1$ and $G=1$.

The observable of interest is the gravitational wave power spectrum, expressed as a fraction of the total density~(\ref{e:GWPowSpe}). It is clear from that formula that the relevant scales are the fluid flow length scale (set by the average bubble separation), the fluid flow lifetime, and the Hubble rate. The power spectrum is determined by the ratios of the fluid flow lifetime to the Hubble time, and the fluid flow scale to the Hubble scale. The role of the bubble wall width is to provide a short-distance cut-off on the power spectrum. A physical transition has $\ell < \Rc \ll \Rbc$, with $\Hc\Rbc$ of order 
$10^{-2}$ at the electroweak scale.
Our simulations assume that the bubble separation is much less than the Hubble length, and that the transition rate is much larger than the Hubble rate, so that expansion can be neglected.

The fluid flow lifetime affects only the amplitude of the acoustically generated gravitational waves. In our simulations one sees that after the transition has completed, the power spectrum grows linearly with time while maintaining its shape.
Hence, apart from a trivial scaling, the relevant parameter for the gravitational wave power spectrum is the fluid flow scale, provided that the wall width and the critical bubble size are much less that the bubble separation. 
The effect of too large a ratio $\ell/\Rbc$ is that there is insufficient dynamic range to observe the power law behaviour of the power spectrum; the effect of too large a ratio $\Rc/\Rbc$ is that there is insufficient time for the fluid flow to approach its asymptotic self-similar profile, which results in too low a value for $\fluidV$. It also tends to obscure the power law behaviour. 
We have seen in our simulations that the ratio $\ell/\Rbc$ needs to be of order $10^{-3}$ in order to reliably distinguish the power law.  
Given our computing resources, this means we are not able to determine the shape of the power spectrum at wave numbers much less than the peak.

In order to test the approach to physical ratios we should explore a scaling of the parameters which shrinks the ratios $\ell/\Rbc$ and $\Rc/\Rbc$
to zero. Such a scaling was given in Eq.~(\ref{scaling}). Its only effect is to alter the width and surface tension of the bubble wall, and hence shrink the size of the critical bubble and the bubble wall width independently of the bubble separation. 
We carried out a simulation scaled with $r=2$ (so that the bubble wall was half the width) and parameters given in Table \ref{t:SimParsRuns} corresponding to the deflagration, 
and compare the resulting gravitational wave power spectra in Figure~\ref{fig:scaledcomp}.
The power spectra are substantially similar, but the $k^{-3}$ power law is clearer in the scaled run where $\ell/\Rbc$ is smaller.  This is consistent with our discussion above, and lends further confidence to our identification of the index of the power law in this case.


\begin{figure}
\begin{centering}
\includegraphics[width=0.4\textwidth,clip=true]{scaledcomp.eps}
\caption{Power spectra with $N_\text{b}=988$ comparing the weak phase transition parameters and friction $\eta/T_c=0.2$ with the results from an equivalent run with the scaled parameters. For clarity, only the power spectra at the end of the phase transition ($t=2500/\Tc$) are shown. Note the $y$-axis scale is different to that used in Fig.~\ref{f:GWPSMulti}, in order to highlight the differences between the two power spectra.}
\label{fig:scaledcomp}
\end{centering}
\end{figure}

Note that the scaling of Eq.~(\ref{scaling}) also reduces the surface tension $\sigma$ as it reduces the bubble wall width, and hence the relative contribution of the scalar field to the total gravitational wave source tensor (\ref{e:TauDef}), as the following argument makes clear. 
 The scalar field's source tensor $\tau^\phi$  is proportional to the product of $\sigma$ with the area per unit volume of the phase boundary, and the area per unit volume is at most of order $1/\Rbc$, which is unaffected by the scaling.  Hence 
 $\tau^\phi \to r^{-1}\tau^\phi$. 
 At the same time, the scale of the fluid source tensor $\tau^\text{f}$ is set by the latent heat of the transition ${\cal L}$, which is independent of $r$.  
Hence the the relative importance of the scalar field to the fluid goes as $\ell/\Rbc$ as it decreases towards physical values.

This is consistent with the argument given in the introduction that the scalar field contributes negligibly to the gravitational waves, as the ratio of the energy in the scalar field to the energy in the fluid goes as the ratio of the volume in the phase boundary to the total volume, which is  $\ell/\Rbc$.

This argument assumes that the bubble walls travel at constant speed, so that the effective surface energy is constant.  Hence if the bubble walls are weakly coupled to the plasma, they can continue to accelerate until they collide~\cite{Bodeker:2009qy}. In this ``run-away'' scenario, scalar fields can contribute importantly to the gravitational radiation.

\section{Discussion and conclusions}

In this paper we have reported on new numerical simulations of the production of gravitational radiation at a first order phase transition in the early universe.
Following standard methods, we model the contents of the universe as a scalar order parameter coupled to a relativistic fluid, with a thermal effective potential~(\ref{e:ScaPot}) and dissipative coupling~(\ref{e:CouTer}). This model captures the essential physics of the transition, which proceeds by the nucleation and growth of bubbles of the low temperature phase. 

The most important parameters of the transition are the latent heat density relative to the total energy density $\strengthPar{T}$, which characterises the strength of the transition, the bubble wall velocity $\vw$, which is determined by $\strengthPar{T}$ and the field-fluid coupling $\eta$, and the bubble nucleation rate parameter $\be$, which determines the average bubble separation.

Most of our simulations are carried out at $\strengthPar{T} \simeq 0.01$, the order of magnitude expected at an electroweak transition, from which we can extrapolate to other values.  We check our extrapolation with a smaller number of simulations at $\strengthPar{T} \simeq 0.1$, and with a scaling argument which changes parameters in the potential without affecting $\strengthPar{T}$.  We simulate	 for a range of phase boundary speeds $\vw$, covering deflagrations and detonations.  Instead of fixing the bubble nucleation rate parameter $\be$, we directly fix the average bubble separation $\Rbc$ by nucleating $N_\text{b} = V/\Rbc^3$ bubbles simultaneously.

We concentrate on the gravitational waves generated by the fluid motion, as the vast majority of the latent heat of the transition is transformed into thermal and kinetic energy of the fluid. 
We show that the gravitational wave density parameter~(\ref{e:OmgwEqn}) is proportional to the fourth power of the mean square fluid velocity, the ratio of lifetime of the source to the Hubble time, and the ratio of length scale of the source to the Hubble length. 

Our results confirm those of a more limited set of simulations reported in Ref.~\cite{Hindmarsh:2013xza}. The fluid kinetic energy is mostly in the form of sound waves generated by the compression or rarefaction of the fluid around the advancing phase boundary. Some rotational flow is generated by the collisions, but at a subdominant level. 
The sound waves remain for as long as we simulate, long after the phase transition completes.  

It was shown in Ref.~\cite{Hindmarsh:2013xza} that when viscosity is included 
the viscous damping time is much longer than the Hubble time  for most phase transitions of interest.
It was argued that the lifetime of the source, the shear stress generated by the sound waves, is approximately the Hubble time. In this paper we detail the calculation which shows that the lifetime parameter $\tau_\text{v}$, controlled by the decay and decorrelation of the shear stresses,  is in fact exactly the Hubble time.

The length scale of the source, approximately the average bubble separation in Ref.~\cite{Hindmarsh:2013xza}, is here measured directly from the fluid flow.
With this refinement, we show that the proportionality constant $\tilde\Omega_\text{GW}$ in the gravitational wave density parameter equation~(\ref{e:OmgwEqn}) varies little between phase transitions with different strength and bubble wall speeds. Indeed, our measurements show that $8\pi\tilde\Omega_\text{GW} = 0.8 \pm 0.1$, where the uncertainly is the root mean square fluctuation between simulations.

Our new simulations are carried out on larger lattices, and give a wider dynamic range between the physical scales set by the average bubble separation and the bubble wall width. We further widen the dynamic range by nucleating all bubbles at the same time, at the slight cost of introducing ``ringing'' in the velocity power spectrum.  With the increased dynamic range we are able to establish clear power laws for both velocity and gravitational wave power spectra between the physical scales. 
For the transitions with $\vw = 0.44$ or below, they are $k^{-1}$ and $k^{-3}$ respectively, where $k$ is the wave number, and steeper for the transitions with higher bubble wall speeds. In order to discern these power laws, we show it is important that the fluid velocity profile around the advancing bubble wall has sufficient time to approach its asymptotic self-similar form.

The $k^{-3}$ (or steeper) power law for gravitational waves contrasts with the prediction of $k^{-1}$ from the standard envelope approximation, which assumes that all the energy in the system is concentrated in a thin shell at the bubble wall, and that the radiation is produced only when the shells interact. We see signs of a $k^{-1}$ power spectrum generated by the scalar field in the initial phase of bubble collision, but this component is subdominant in our simulations, and would be completely negligible when extrapolated to the scale separation in a thermal phase transition.

The envelope approximation generically predicts far less gravitational radiation than is actually produced. This under-prediction stems from the incorrect modelling of the source as being the colliding bubble walls.  Instead, the main source is the overlapping sound waves which are left behind after the transition has completed.  We argued in \cite{Hindmarsh:2013xza} that this means that the gravitational wave energy density is boosted by the ratio of the lifetime parameter of the shear stress to the duration of the collision, which goes parametrically as $(\vwL_\text{f}\Hc)^{-1}$.  In this paper we studied the numerical factor in this ratio by a careful comparison of the quantities in the envelope approximation formula (\ref{e:EnvAppFor}) to the acoustic generation formula (\ref{e:OmgwEqn}).  We show that the numerical factor is of order unity, and hence we can confirm that the gravitational wave signal is boosted by the ratio of the Hubble time to the phase transition duration, which is two orders of magnitude or more for a typical first order electroweak transition \cite{Hogan:1984hx,Enqvist:1991xw}.

Our simulations shed new light on gravitational waves from phase transitions in the early universe.  They show that the envelope approximation needs to be replaced, both as a model and as a formula.  Instead, we should model the gravitational wave generation in terms of overlapping sound waves.

With this new ``acoustic'' model of gravitational wave generation, we have developed a quantitative understanding of the gravitational wave density parameter (\ref{e:OmgwEqn}), as a function of the mean square fluid velocity and the mean bubble separation.  We can estimate the mean fluid velocity from hydrodynamic considerations \cite{Espinosa:2010hh}, and the mean bubble separation from the nucleation rate parameter $\be$ \cite{Enqvist:1991xw}.
We have numerically determined that the gravitational wave power spectrum is a power law on the high wavenumber side of the peak, and shown that it is steeper than the $k^{-1}$  indicative of a vacuum transition.  Hence potential future observations of such a gravitational wave spectrum will allow us to distinguish between a thermal and a vacuum transition. 

Much remains to be done. We noted that we need larger simulations to trace out the shape of the power spectrum at wavenumbers lower than the peak value, and to determine the index of the power spectrum for the transitions with faster bubble walls. They may also help in the search for bubble wall instabilities identified in~\cite{Link:1992dm,Huet:1992ex,Megevand:2013yua}.
We also need to develop a theoretical understanding of the shape or the power spectrum, and most importantly making accurate quantitative predictions for future gravitational wave observatories. 

