Numerical Simulations of Gravitational Waves from Early-Universe Turbulence

TL;DR

This paper uses direct numerical simulations to study gravitational waves generated by early-universe turbulence, revealing a universal spectral form and implications for detectability by LISA at electroweak scales.

Contribution

It provides the first non-analytic numerical analysis of gravitational wave production from early-universe turbulence, showing a universal spectral shape and updated detectability predictions.

Findings

Universal spectral form below the peak frequency.

Larger low-frequency gravitational wave signals than previously estimated.

Detectability of electroweak-scale turbulence by LISA depends on energy in plasma motions.

Abstract

We perform direct numerical simulations of magnetohydrodynamic turbulence in the early universe and numerically compute the resulting stochastic background of gravitational waves and relic magnetic fields. These simulations do not make the simplifying assumptions of earlier analytic work. If the turbulence is assumed to have an energy-carrying scale that is about a hundredth of the Hubble radius at the time of generation, as expected in a first-order phase transition, the peak of gravitational wave power will be in the mHz frequency range for a signal produced at the electroweak scale. The efficiency of gravitational wave (GW) production varies significantly with how the turbulence is driven. Detectability of turbulence at the electroweak scale by the planned Laser Interferometer Space Antenna (LISA) requires anywhere from 0.1% to 10% of the thermal plasma energy density to be in plasma motions or magnetic fields, depending on the model of the driving process. Our results predict a new universal form below the spectral peak frequency that is shallower than previously thought. This implies larger values of the GW energy spectra in the low-frequency range. This extends the range where turbulence is detectable with LISA to lower frequencies, corresponding to higher energy scales than the assumed energy-carrying scale.

Equations

$$\displaystyle\biggl{(}\partial^{2}_{t_{\rm phys}}+3H(t)\partial_{t_{\rm phys}}-c^{2}\mathchoice{\mbox{\boldmath$\displaystyle\nabla$}}{\mbox{\boldmath$\textstyle\nabla$}}{\mbox{\boldmath$\scriptstyle\nabla$}}{\mbox{\boldmath$\scriptscriptstyle\nabla$}}{}^{2}_{\text{phys}}\biggr{)}h_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t)$$

$$\displaystyle=\frac{16\pi G}{c^{2}}T^{\text{TT}}_{ij,\text{phys}}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t),$$

$$\left(\partial_{t}^{2}-c^{2}\mathchoice{\mbox{\boldmath$\displaystyle\nabla$}}{\mbox{\boldmath$\textstyle\nabla$}}{\mbox{\boldmath$\scriptstyle\nabla$}}{\mbox{\boldmath$\scriptscriptstyle\nabla$}}{}^{2}\right)h_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t)=\frac{16\pi G}{ac^{2}}T_{ij}^{\text{TT}}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t),$$

$${\mathcal{E}}_{\mathrm{crit}}^{\ast }=\frac{3H_{\ast }^2{c}^2}{8\pi G}\approx {\mathcal{E}}_{\mathrm{rad}}^{\ast }=\frac{{\pi }^2{g}_{\ast }(k_{\mathrm{B}}{T}_{\ast }{)}^4}{30(\hslash c{)}^3},$$

$$\left(\partial^{2}_{\bar{t}}-\bar{\mathchoice{\mbox{\boldmath$\displaystyle\nabla$}}{\mbox{\boldmath$\textstyle\nabla$}}{\mbox{\boldmath$\scriptstyle\nabla$}}{\mbox{\boldmath$\scriptscriptstyle\nabla$}}}{}^{2}\right)h_{ij}(\bar{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{}},\bar{t})\,=\,{6\over\bar{t}}\bar{T}_{ij}^{\rm TT}(\bar{\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{}},\bar{t}).$$

$$h_{c}^{2}(t)=\frac{1}{2}\bigl{\langle}\bigl{(}h_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t)\bigr{)}^{2}\bigr{\rangle}=\frac{1}{2t^{2}}\bigl{\langle}h^{2}_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t)\bigr{\rangle},$$

$${\cal E}_{\rm GW}(t)=\frac{c^{2}}{32\pi G}\bigl{\langle}\bigl{(}\partial_{t_{\rm phys}}h_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t)\bigr{)}^{2}\bigr{\rangle},$$

$${\Omega}_{\rm GW}(t)=\frac{1}{12t^{4}}\bigl{\langle}\bigl{(}\partial_{t}h_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t)-h_{ij}(\mathchoice{\mbox{\boldmath$\displaystyle x$}}{\mbox{\boldmath$\textstyle x$}}{\mbox{\boldmath$\scriptstyle x$}}{\mbox{\boldmath$\scriptscriptstyle x$}}{},t)/t\bigr{)}^{2}\bigr{\rangle}.$$

$$f=\frac{c{k}_{\mathrm{phys}}}{2\pi }=\frac{H_{\ast }{a}_0^{-1}}{2\pi }\bar{k}.$$

$$T_{ij}=\frac{4}{3}\frac{\rho u_{i}u_{j}}{1-\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}{}^{2}}-B_{i}B_{j}+(\rho/3+\mathchoice{\mbox{\boldmath$\displaystyle B$}}{\mbox{\boldmath$\textstyle B$}}{\mbox{\boldmath$\scriptstyle B$}}{\mbox{\boldmath$\scriptscriptstyle B$}}{}^{2}/2)\delta_{ij},$$

$$\frac{\partial \ln \rho }{\partial t}$$

$$\displaystyle{\partial\mathchoice{\mbox{\boldmath$\displaystyle u$}}{\mbox{\boldmath$\textstyle u$}}{\mbox{\boldmath$\scriptstyle u$}}{\mbox{\boldmath$\scriptscriptstyle u$}}{}\over\partial t}$$

$$\displaystyle{\partial\mathchoice{\mbox{\boldmath$\displaystyle B$}}{\mbox{\boldmath$\textstyle B$}}{\mbox{\boldmath$\scriptstyle B$}}{\mbox{\boldmath$\scriptscriptstyle B$}}{}\over\partial t}$$

$$H_{\ast }\approx 2.066\times 10^{10}{\mathrm{s}}^{-1}{\left(\frac{k_{\mathrm{B}}{T}_{\ast }}{100\mathrm{GeV}}\right)}^2{\left(\frac{g_{\ast }(T_{\ast })}{100}\right)}^{1/2},$$

$$a_0\approx 1.254\times 10^{15}\left(\frac{k_{\mathrm{B}}{T}_{\ast }}{100\mathrm{GeV}}\right){\left(\frac{g_{\mathrm{S}}(T_{\ast })}{100}\right)}^{1/3},$$

Full text

††preprint: NORDITA-2019-024

I Introduction

A period of turbulence in the early universe can produce a stochastic background of gravitational waves (GWs). The turbulence that produces GW radiation might arise from the dynamics of a first-order phase transition Witten:1984rs ; Hogan:1986qda ; Kamionkowski:1993fg , from the dynamics of primordial magnetic fields Deryagin , or from the dynamical coupling of primordial magnetic fields and the highly conducting primordial plasma Brandenburg:1996fc ; Christensson:2000sp ; Kahniashvili:2010gp ; BKMRPTV17 . Analytic estimates suggest that turbulence generated by an electroweak phase transition can produce GWs within the detectable amplitude and frequency range of the Laser Interferometer Space Antenna (LISA) if the turbulent energy density is roughly 1% of the total energy density of the Universe at that time Kosowsky:2001xp ; KGR05 ; Caprini:2006jb ; Gogoberidze:2007an .

However, the aforementioned analytic estimates make a number of simplifying assumptions. Turbulence is assumed to be hydrodynamic with a typical Kolmogorov power spectrum and a duration set by a small fraction of the Hubble time, omitting the effect of the expansion of the Universe during the turbulent period. The inclusion of magnetic fields can extend the frequency range of the resulting GWs due to the transfer of power to larger scales PFL76 ; Brandenburg:1996fc . These turbulence models depend on the temporal correlation function of the turbulent velocity field, which was assumed in earlier works and not computed from magnetohydrodynamic (MHD) simulations. An accurate treatment of these effects is essential for reliably establishing the spectral shape of the resulting GW background and its detectability with upcoming detectors Romano:2016dpx . A proper understanding of turbulent sourcing of GWs is especially relevant for using LISA to constrain the parameter space of a first-order phase transition HHRW17 .

If primordial magnetic fields were present during the early universe, they could dynamically enhance turbulent plasma motions and serve as an additional source of GWs Deryagin ; Kahniashvili:2008er ; Caprini:2009yp . Such magnetic fields can persist until the present epoch. Lower bounds on the strength of magnetic fields obtained from observations of TeV blazar spectra NV10 are suggestive of the existence of these primordial fields.

Numerical simulations are required to make progress beyond previous analytic estimates, as pointed out in a recent report of the LISA cosmology working group LISA19 . We present here the results of direct numerical simulations of MHD turbulence and the resulting stochastic GW spectra. Given that the turbulent energy densities are below 10% of the total energy density of the Universe, the bulk motions in our simulations are subrelativistic, but the equation of state is still a relativistic one. We use the Pencil Code PC , a sixth-order finite-difference code using third-order time stepping for the MHD equations and a novel approach for numerically solving the GW equation, which is discussed in a separate paper Pol:2018pao .

The present paper is arranged as follows. Section II presents the equations that describe the production of GWs and the dynamics of the magnetic and velocity fields during the radiation-dominated epoch of the early universe. The initial conditions and the setup of the simulations are presented in Sec. III. The results of the numerical simulations are presented and compared with previous analytic estimates, and the prospects of detectability with LISA are discussed in Sec. IV.

Electromagnetic quantities are expressed in Lorentz-Heaviside units where $\mu_{0}=1$. Einstein index notation is used, so summation is assumed over repeated indices. Latin indices $i$ and $j$ refer to spatial coordinates 1 to 3.

II Equations

We assume the evolution of the background universe to be described by the spatially flat, homogeneous, and isotropic Friedmann-Lemaître-Robertson-Walker metric $g_{ij}=a^{2}\delta_{ij}$, with $a$ being the scale factor. The expansion of the Universe described by the temporal evolution of $a$ leads to a dilution of radiation energy density and magnetic fields, and to the damping of the GW amplitude. It is convenient to scale out the effects of expansion by using conformal time $t$, defined as $dt=dt_{\rm phys}/a$, and comoving coordinates \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{}=\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{}_{\rm phys}/a and MHD fields. The physical coordinates and time are expressed as \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{}_{\rm phys} and $t_{\rm phys}$. The ultrarelativistic equation of state $p=\rho c^{2}/3$ is valid during radiation domination. This leads to a linear evolution of the scale factor with conformal time $t$ as the solution to the Friedmann equations Friedmann1922 .

II.1 Gravitational wave equation

We consider small tensor-mode perturbations $a^{2}h_{ij}^{\rm phys}$ over the background metric $g_{ij}$. The GW equation is then Gri74 ; Mukhanov

[TABLE]

where $c$ is the speed of light, $G$ is Newton’s gravitational constant, and $H=a^{-2}\partial_{t}a$ is the Hubble rate. The transverse and traceless stress tensor $T_{ij,{\rm phys}}^{\rm TT}$ sources the gauge-free metric perturbations MTW73 .

The introduction of comoving coordinates, conformal time, comoving stress tensor $T_{ij}^{\rm TT}=a^{4}T_{ij,{\rm phys}}^{\rm TT}$, and scaled strains $h_{ij}=ah_{ij}^{\rm phys}$ simplifies Eq. (1) to Pol:2018pao

[TABLE]

where the omitted damping term $h_{ij}\,a^{-1}\partial_{t}^{2}a$ is zero during radiation-domination.

The stress tensor $\bar{T}_{ij}^{\rm TT}$ is normalized by the energy density at the initial time of turbulence generation $t_{*}$, which, during the radiation-dominated era, is

[TABLE]

where $T_{*}$, $g_{*}$, and $H_{\ast}$ are the temperature, the number of relativistic degrees of freedom, and the Hubble rate, respectively, $k_{\rm B}$ is the Boltzmann constant, and $\hbar$ is the reduced Planck constant.

We also normalize conformal time by $t_{*}$ and define $\bar{t}=t/t_{*}$. Because of the linear evolution of the scale factor $t_{\ast}=H_{\ast}^{-1}$, where $a_{\ast}=1$ has been used. The scale factor is just $a=\bar{t}$. Note that in this convention, $a_{0}\neq 1$ at the present time. Similarly, the space coordinates are normalized as \bar{\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}}{}=\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{}H_{\ast}/c. Equation (2) reduces then to the normalized GW equation,

[TABLE]

This is the wave equation that we solve within the Pencil Code. From now on, we omit overbars and always refer to normalized variables, unless stated otherwise.

The use of normalized variables has the advantage that our simulations can readily be applied to different energy scales where our equations are applicable, i.e., after the electroweak phase transition and within the radiation-dominated epoch. Through a simple set of steps, the final diagnostic frequency diagram can then easily be assembled from the normalized results.

II.2 Gravitational wave characteristics

The characteristic amplitude of GWs is defined as Mag00

[TABLE]

where angle brackets denote averaging over the physical volume, and the second equality is true during the radiation-dominated epoch with the normalization described above. The spectrum of the characteristic amplitude is defined following Ref. Mag00 , such that $\int h_{\rm c}^{2}(k,t)\,d\ln k=h_{\rm c}^{2}(t)$; see details in Ref. Pol:2018pao . The integration to compute $h_{\rm c}$ is performed over wave numbers $k$ from 0 to $\infty$. Since the spectral function $h_{\rm c}(k)$ is defined to be integrated in $\ln k$, the limits of integration become $-\infty$ to $\infty$. Note that $k$ refers to the normalized wave number $\bar{k}$, consistently given by $\bar{k}=ck/H_{\ast}=ack_{\rm phys}/H_{\ast}$, where again we omit the overbar from now on.

In the absence of turbulent sources, and neglecting details of the GW transfer function due to evolving relativistic degrees of freedom and transitions between radiation, matter, and dark energy dominations (see Ref. SS18 ), the characteristic amplitude $h_{\rm c}(k,t)$ dilutes due to the expansion of the Universe as $a^{-1}$. Hence, the relic observable amplitude at the present time $h_{\rm c}(k)$ is the amplitude at the end of the simulation $t_{\rm end}$ diluted by a factor $t_{\rm end}/a_{0}$. Note that $t_{\rm end}$ is assumed to be within the radiation-dominated epoch, such that the computed numerical results and the described normalization are valid. The value of the scale factor $a_{0}$ is obtained assuming adiabatic expansion of the Universe, i.e., such that $g_{\rm S}\,T^{3}a^{3}$ stays constant, where $g_{\rm S}$ is the number of adiabatic degrees of freedom.

The physical energy density carried by the GWs ${\cal E}_{\rm GW}(t)$ is defined as Mag00

[TABLE]

which we normalize by the radiation energy density, ${\Omega}_{\rm GW}(t)={\cal E}_{\rm GW}(t)/{\cal E}_{\rm rad}^{\ast}$. In terms of conformal time $t$ and scaled strains h_{ij}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{},t) during the radiation-dominated epoch, our normalization leads to

[TABLE]

The GW energy density spectrum ${\Omega}_{\rm GW}(k,t)$ is defined as in Ref. Mag00 , such that $\int{\Omega}_{\rm GW}(k,t)\,d\ln k={\Omega}_{\rm GW}(t)$; see details in Ref. Pol:2018pao . This is the standard normalization that we use within the Pencil Code. However, when we are interested in the observable relic signal, it is useful to normalize by the critical energy density at the present time, ${\cal E}_{\rm crit}^{0}=(3H_{0}^{2}c^{2})/(8\pi G)$, where $H_{0}=100\,h_{0}\,{\rm km}\,\,{\rm s}^{-1}\,{\rm Mpc}^{-1}\approx 3.241\times 10^{-18}\,h_{0}\,{\rm s}^{-1}$ is the Hubble rate at the present time. We use $h_{0}^{2}{\Omega}_{\rm GW}$ to quote our results independently of the uncertainties in the value of $h_{0}$ Mag00 . The GW energy density dilutes due to the expansion of the Universe as $a^{-4}$. Hence, the relic observable at the present time ${\Omega}_{\rm GW}(k)$ is the energy density in Eq. (7) at the end of the simulation $t_{\rm end}$, reduced by a factor $(H_{\ast}/H_{0})^{2}(t_{\rm end}/a_{0})^{4}\sim g_{\ast}^{-1/3}$. During the radiation-dominated epoch, $a^{4}{\Omega}_{\rm GW}$ is constant. However, since $a$ does not evolve as $t$ all the way to the present time, this factor is not unity. Finally, we express the GW amplitude $h_{\rm c}(f)$ and the energy spectra ${\Omega}_{\rm GW}(f)$, which are observables at the present time—as a function of the physical frequency, shifted to the present time. The frequency is related to $\bar{k}$ through

[TABLE]

II.3 MHD equations

The GW equation is sourced by the stress tensor $T_{ij}$. In particular, we consider GWs sourced by MHD turbulence. Starting with initial conditions for the plasma velocity and magnetic fields at the starting time of the turbulence period, we numerically solve for the dynamics of early-universe MHD turbulence using the Pencil Code. At each time step, we compute the spatial Fourier components of the stress tensor of a relativistic perfect fluid,

[TABLE]

where $\textstyle u$ is the plasma velocity and $\textstyle B$ is the magnetic field. The MHD fields $\rho=a^{4}\rho_{\rm phys}$ and \mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}=a^{2}\mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}_{\rm phys} are expressed as comoving variables. We emphasize that in MHD, the Faraday displacement current is omitted. This implies that electric fields do not contribute to the stress tensor Brandenburg:1996fc ; Subramanian15 .

The MHD equations for an ultrarelativistic gas in a flat expanding universe in the radiation-dominated era after the electroweak phase transition are given by Brandenburg:1996fc ; BKMRPTV17 ; DN13

[TABLE]

where {\sf S}_{ij}={\textstyle{1\over 2}}(u_{i,j}+u_{j,i})-{\textstyle{1\over 3}}\delta_{ij}\mathchoice{\mbox{\boldmath\displaystyle\nabla}}{\mbox{\boldmath\textstyle\nabla}}{\mbox{\boldmath\scriptstyle\nabla}}{\mbox{\boldmath\scriptscriptstyle\nabla}}{}\cdot\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}{} are the components of the rate-of-strain tensor with commas denoting partial derivatives, $\nu$ is the kinematic viscosity, and $\eta$ is the magnetic diffusivity. Energy can be injected into velocity and magnetic fields through ponderomotive and electromagnetic forces ${\cal F}$ and ${\cal E}$, respectively.

All variables are normalized with the appropriate powers of the radiation energy density and the Hubble rate, both at the time of generation: $\bar{\rho}=\rho c^{2}/{\cal E}_{\rm rad}^{\ast}$, \bar{\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}{}}=\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}{}/c, \bar{\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}{}}=(c/H_{*})\,\mathchoice{\mbox{\boldmath\displaystyle J}}{\mbox{\boldmath\textstyle J}}{\mbox{\boldmath\scriptstyle J}}{\mbox{\boldmath\scriptscriptstyle J}}{}/\sqrt{{\cal E}_{\rm rad}^{\ast}}, \bar{\mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}}=\mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}/\sqrt{{\cal E}_{\rm rad}^{\ast}}, $\bar{\eta}=H_{*}\eta/c^{2}$, $\bar{\nu}=H_{*}\nu/c^{2}$, \bar{\mbox{\boldmath{\cal F}}{}}=\mbox{\boldmath{\cal F}}{}/(H_{\ast}c), and \bar{\mbox{\boldmath{\cal E}}{}}=\mbox{\boldmath{\cal E}}{}/\sqrt{c^{2}{\cal E}_{\rm rad}^{\ast}}, where the overbars have been dropped on Eqs. (9)–(12). In addition, similar to $\rho$ and $\textstyle B$, the current density $\textstyle J$ is comoving, i.e., scaled with $a^{3}$. The physical value of the magnetic diffusivity $\eta$ at the electroweak phase transition is given in Eq. (9) of Ref. BSRKBFRK17 : $\eta\approx 4\times 10^{-9}\,(k_{\rm B}T_{\ast}/100\,{\rm GeV})^{-1}\,{\rm cm}^{2}/\,{\rm s}$. This corresponds to $9.2\times 10^{-20}$ in our normalized units.

The energy densities of the magnetic and velocity fields are computed as {\Omega}_{\rm M}(t)=\langle\mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}^{2}\rangle/2 and {\Omega}_{\rm K}(t)=\langle\rho\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}{}^{2}\rangle/2. We define the magnetic and kinetic energy spectrum such that $\int\Omega_{\rm M,K}(k,t)\,d\ln k=\Omega_{\rm M,K}(t)$. Here, $\Omega_{\rm M,K}(k,t)$ are the spectra in terms of logarithmic wave number intervals. They are defined analogously to ${\Omega}_{\rm GW}(k,t)$; see Refs. Mag00 ; DC03 ; Pol:2018pao .

III Numerical simulations

To compute the resulting GW production, we evolve the strains in Fourier space using Eq. (4), assuming a constant source term during the length of one time step of the MHD evolution. This assumption is accurate for time steps small enough to guarantee numerical stability of the MHD equations, and it allows much longer time steps than what is required by a direct numerical simulation; see Sec. 2.6 of Ref. Pol:2018pao for a discussion of this new method, which is described there as approach II.

It turns out that the GW energy production ceases some time after the kinetic and/or magnetic energies have started to decay. The GW spectrum is then statistically steady. We continue our simulations to gather sufficient statistics to compute accurate averages for the GW spectra in comoving variables (more specifically, until the GW spectra start to fluctuate around a steady state and we have computed at least one period). In all our simulations, this occurs well within the radiation-dominated era. The last time of the numerical simulations is what we have previously called $t_{\rm end}$.

To study the sensitivity to initial conditions, we have performed several sets of simulations with different physical models for driving plasma motions. The motivation for the different types of primordial magnetic fields obtained below is given in Ref. BKMRPTV17 , where their subsequent evolution and observational constraints are discussed. The physical magnetic diffusivity $\eta$ of the early universe is much smaller than what we can accurately simulate. For similar reasons, thermal and radiative diffusion were not included in the equations above. We fix the viscosity $\nu=\eta$ and choose it to be as small as possible, but still large enough such that the inertial range of the computed spectra is appropriately resolved BKMRPTV17 . If the much smaller physical values were used instead, much larger numerical resolution would be required and, the inertial range of the turbulence would extend to higher frequencies. Those higher frequencies are of little observational interest since the GW amplitude at those frequencies would be very low, as seen from our spectra shown below.

Our full set of runs is summarized in Table LABEL:Tsummary. For all of the calculations, we assume \mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}{}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{})=\bm{0} initially. In set I (runs ini1–3), \mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{}) is initialized as a fully helical (indicated by “y” under “hel”) Gaussian random field with magnetic energy spectrum $\Omega_{\rm M}(k)\propto k^{5}$ for $k<k_{\ast}$, corresponding to a solenoidal causally generated field, and $\Omega_{\rm M}(k)\propto k^{-2/3}$ (Kolmogorov spectrum) for $k>k_{\ast}$ where $k_{\ast}$ is the wave number at which the magnetic energy is injected. In set II (runs hel1–4 and noh1–2), \mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{})=\bm{0} initially, but it is then numerically driven by applying an electromotive force ${\cal E}$ during $1\leq t\leq t_{\max}$ in the induction Eq. (12) consisting of random and nearly monochromatic waves around wave number $k_{\ast}$.

The driving force field is taken as either fully helical (runs hel1–4) or nonhelical (noh1–2), using a forcing term quantified by ${\cal E}_{0}$ described in Refs. Bran01 ; BDS02 ; see Table LABEL:Tsummary for values of ${\cal E}_{0}$ and $t_{\max}$. The initial number of eddies per horizon length at the driving scale is $N\equiv k_{\ast}/2\pi$, usually taken to be between 1 and 100 for the first-order electroweak phase transition TWW92 .

We compute the decaying MHD turbulent motions, with no forcing term, for times $t>t_{\max}$, where $t_{\max}$ has been defined for set II of runs, and it can be assumed to be $t_{\max}=1.0$ for the runs corresponding to set I. We arrange the simulations such that the maximum total magnetic energy density ${\Omega}_{\rm M}^{\rm max}$ integrated over all wave numbers is a specified fraction of the radiation energy density. We take values in the range $10^{-3}$ to $10^{-1}$. The lower limit is required for the turbulence to be the dominant source of GWs during a first-order phase transition according to analytic estimates Nicolis:2003tg , and the higher limit is imposed on magnetic fields due to their effect on big bang nucleosynthesis Kahniashvili:2009qi ; Yamazaki:2012jd . Recently, values up to 10% have been obtained for magnetogenesis lattice simulations ZVF19 .

We also consider set III (runs ac1–3) with initial \mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{})=\mathchoice{\mbox{\boldmath\displaystyle u}}{\mbox{\boldmath\textstyle u}}{\mbox{\boldmath\scriptstyle u}}{\mbox{\boldmath\scriptscriptstyle u}}{}(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{})={\bf 0}, using irrotational or “acoustic” hydrodynamic turbulence. In this case, the forcing ${\cal F}$ appears as an additional term in the momentum Eq. (11), which acts during $1\leq t\leq t_{\max}$. The forcing term, with amplitude ${\cal F}_{0}$, is computed as the gradient of a combination of Gaussian random potentials \phi\propto\exp[-(\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{}-\mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{}_{i})^{2}/R^{2}] centered at random positions \mathchoice{\mbox{\boldmath\displaystyle x}}{\mbox{\boldmath\textstyle x}}{\mbox{\boldmath\scriptstyle x}}{\mbox{\boldmath\scriptscriptstyle x}}{}_{i} of the domain. This results in a number of eddies $N=(\pi R)^{-1}$ MB06 .

We choose solenoidal and irrotational forcing fields in sets II and III, respectively, for comparison purposes. In all of our runs, the size of the cubic domain $L$ is taken to be $L=2\pi/N$, such that the lower wave number in the computed spectra corresponds to $N$. We evolve the dynamical equations on a mesh of $1152^{3}$ grid points.

IV Results

IV.1 The turbulent GW energy spectrum

We begin by discussing the spectrum for a run of set I. Figure 1 shows the resulting magnetic field and GW energy spectra for such a case with a Kolmogorov-type spectrum. In this case, the magnetic field has a Batchelor ($k^{4}$) spectrum in the subinertial range and a Kolmogorov ($k^{-5/3}$) spectrum in the inertial range. For the resulting stress, this corresponds to a white noise ($k^{2}$) spectrum in the subinertial range and to the same Kolmogorov power law in the inertial range; see the Appendix and Ref. BB19 . The GW energy density shows a spectrum proportional to $k^{-2}$ with respect to the spectrum of the stress. This can be understood by comparing the order of $k$ of the different terms in Eq. (1). The third and fourth terms of the equation in Fourier space are c^{2}k_{\rm phys}^{2}\tilde{h}^{\rm phys}_{ij}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t), and (16\pi G/c^{2})\tilde{T}_{ij,{\rm phys}}^{\rm TT}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t)/t. Therefore, if one assumes these terms to be of the same order, then k^{4}\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t)\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t)\sim\tilde{T}_{ij,{\rm phys}}^{\rm TT}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t)\tilde{T}_{ij,{\rm phys}}^{\rm TT}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t)/t^{2}.

We shell integrate both sides to obtain a term proportional to the spectrum, $k^{4}h_{\rm c}^{2}(k,t)\sim{\Omega}_{\rm T}(k,t)/t^{2}$. We define the stress spectrum ${\Omega}_{\rm T}(k)/k=E_{\rm T}(k)$ analogously to ${\Omega}_{\rm M}$; see the Appendix. On the other hand, the first term is \partial_{t_{\rm phys}}^{2}\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t) in Fourier space, and \omega^{2}\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},\omega)=c^{2}k^{2}_{\rm phys}\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},\omega) if we Fourier transform this term also in time. This leads again to a similar relation, although now in frequency space: $k^{4}h_{\rm c}^{2}(k,\omega)\sim\int{\Omega}_{\rm T}(k,t)\,e^{-i\omega t}/t^{2}\,dt$. The ${\Omega}_{\rm GW}(k,t)$ spectrum is computed by shell integration of \partial_{t_{\rm phys}}\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t)\,\partial_{t_{\rm phys}}\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},t), which is c^{2}k_{\rm phys}^{2}\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},\omega)\tilde{h}_{ij}^{\rm phys}(\mathchoice{\mbox{\boldmath\displaystyle k}}{\mbox{\boldmath\textstyle k}}{\mbox{\boldmath\scriptstyle k}}{\mbox{\boldmath\scriptscriptstyle k}}{},\omega) in the frequency domain. Hence, we have the asymptotic relation ${\Omega}_{\rm GW}(k,\omega)\sim k^{2}h_{\rm c}^{2}(k,\omega)$. This leads to the observed behavior ${\Omega}_{\rm GW}(k)\sim{\Omega}_{\rm T}(k)/k^{2}$ for any fixed time or frequency.

Figure 1 shows a GW spectrum ${\Omega}_{\rm GW}(k)$ that asymptotically falls off faster by a $k^{2}$ factor than the magnetic spectrum ${\Omega}_{\rm M}(k)$ in the inertial range. This is explained by noting that ${\Omega}_{\rm M}(k)$ and ${\Omega}_{\rm T}(k)$ follow the same power law in the inertial range. Hence, ${\Omega}_{\rm GW}(k)\sim{\Omega}_{\rm T}(k)/k^{2}\sim{\Omega}_{\rm M}(k)/k^{2}$. For wave numbers below the spectral peak $k_{\rm GW}\approx 2k_{\ast}$, the GW spectrum, ${\Omega}_{\rm GW}(k)/k$, becomes essentially flat. This small-$k$ behavior in ${\Omega}_{\rm GW}(k)$ can be traced back to the white noise ($k^{2}$) spectrum of the magnetic stress $E_{\rm T}(k)$, which seems to emerge even when the magnetic field itself has a spectrum $E_{\rm M}(k)={\Omega}_{\rm M}(k)/k$ steeper than $k^{2}$ in the subinertial range; see the Appendix. This argument shows that the scaling of ${\Omega}_{\rm GW}(k)$ with $k^{3}$ obtained in previous analytical estimates as in, e.g., Ref. Gogoberidze:2007an , is not expected for the turbulent developed spectrum. The subinertial power law ${\Omega}_{\rm GW}(k)\sim k$ is a novel result from our simulations that was not obtained in previous analytical estimates.

The characteristic amplitude has the following asymptotic behavior: $h_{\rm c}(k)\sim{\Omega}_{\rm T}^{1/2}(k)/k^{2}\sim{\Omega}_{\rm GW}^{1/2}(k)/k$ for a fixed instant of time. Looking at Fig. 1 of Ref. Gogoberidze:2007an , we see that their subinertial range slope in $h_{\rm c}(k)$ is 1/2. This slope in $h_{\rm c}(k)$ corresponds to the +3 slope in ${\Omega}_{\rm GW}(k)$ mentioned above. At high frequencies, our spectrum $h_{\rm c}(k)$ has a slope of $-7/3$ corresponding to a magnetic spectrum of Kolmogorov type. This agrees with what has been obtained in recent analytic work NSS18 ; SSS19 , although earlier work Gogoberidze:2007an reported a slope of $-10/3$, which we would obtain if we used a small magnetic Reynolds number, which results in a $k^{-8/3}$ Golitsyn-type spectrum for the magnetic field ${\Omega}_{\rm M}(k)$; see Table LABEL:Tslopes.

The emergence of a flat GW spectrum in Fig. 1 makes one wonder how this can be reconciled with the principle of causality. We recall that the reason for the magnetic energy spectrum to have a $k^{4}$ subinertial range is indeed causality combined with the fact that the magnetic field is solenoidal. A flat GW energy spectrum would imply that there must have been GW energy immediately at the largest possible length scales, which would be unphysical. Therefore, we have performed numerical simulations with very small time steps to study how the novel low-$k$ spectrum develops at initial times. In Fig. 2, we show that initially, the GW spectrum is indeed proportional to $k^{2}$ and that similar spectra are also being reproduced in 10 and 50 times larger domains. It is only during the first few time steps of the numerical simulation that the GW spectrum is still proportional to $k^{2}$, until the $k^{0}$ scaling extends over the range between the stirring scale and the lower wave number in our simulations. We observe the development of this flat spectrum for the different sizes of the numerical domain. This rules out the possibility that this scaling is due to numerical artifacts, indicating that the flat spectrum is physical and emerges only later.

The time it takes for the change of slope to occur below the horizon scale is much shorter than the time it takes for the GW spectrum to become stationary. Therefore, we conclude that the $+2$ slope in ${\Omega}_{\rm GW}(k)/k$ is not relevant for the characterization of the signal. To demonstrate this further, we show in Fig. 3 that at small $k$, $\Omega_{\rm GW}(k,t)/k$ grows with $t$ proportional to $k^{2}(t-t_{*})^{2}$, where $t_{*}=1$ in normalized units, and reaches a constant level that is independent of $k$ and is given by the white noise spectrum of the source at large scales. This is demonstrated for wave numbers as small as a few times the Hubble horizon wave number, $k=1$, $2$, $4$, and $8$.

IV.2 Spectra from the electroweak phase transition

The GW energy density $h_{0}^{2}{\Omega}_{\rm GW}(f)$ and the characteristic strain amplitude $h_{\rm c}(f)$ are shown in Fig. 4 for runs ini1–3 as a function of the frequency $f$, all shifted to the present time as defined in Sec. II.1. These are obtained by scaling the computed normalized GW spectra to the physical spectra produced at the electroweak scale. If we take $T_{\ast}=100\,{\rm GeV}$ for the temperature, and $g_{\ast}\approx g_{\rm S}=100$ for the number of relativistic and adiabatic degrees of freedom at the electroweak phase transition, the Hubble rate is [see Eq. (3)]

[TABLE]

where the proportionality factors $T$ and $g_{\ast}$ are kept as a parameter due to the uncertainty of the exact values. Our convention of setting $a_{\ast}=1$ leads to the following value of $a_{0}$ computed assuming adiabatic expansion of the Universe [see text above Eq. (6)],

[TABLE]

where we have used the values $g_{\rm S}=3.91$ and $T_{0}=2.73\,{\rm K}$ at the present time.

For different scenarios, the results scale in the following way: The frequency shifts proportional to $Tg_{\ast}^{1/2}g_{\rm S}^{-1/3}$, the strain amplitude varies with $T^{-1}g_{\rm S}^{-1/3}$, and the GW energy density with $g_{\ast}g_{\rm S}^{-4/3}$. Likewise, the GW strain amplitude is proportional to the stirring scale $N^{-3/2}$ and the frequency is proportional to $N$ Gogoberidze:2007an .

The slopes of ${\Omega}_{\rm GW}(f)$ in Fig. 4 are consistent with those in Fig. 1, where ${\Omega}_{\rm GW}(k)$ was shown, due to the dispersion relation; see Eq. (8). The spectrum of $h_{\rm c}(f)$ shows scaling with $f^{-1/2}$ for low frequencies and with $f^{-7/3}$ in the inertial range, as expected. As discussed above, we expect the subinertial slope of $-1/2$ to eventually turn over a slope of $+1/2$ as $f$ and time decrease, due to the lack of causality on scales larger than the horizon. However, the simulation domains are smaller than the horizon scale, so this turnover is not observed.

The analytic approximation in Ref. Gogoberidze:2007an gives a peak value $h_{\rm c}\approx 4\times 10^{-20}$ at $1\,{\rm mHz}$ for their largest Mach number of unity (see Fig. 1 of Ref. Gogoberidze:2007an ). By comparison, for our run ini1, the spectrum shows an intermediate peak at $3\,{\rm mHz}$ with $h_{\rm c}\approx 0.7\times 10^{-21}$; see Fig. 4.

In runs ini1–3, GWs are produced by the sudden emergence of a magnetic field. In reality, this will be a gradual process, as modeled by sets II and III of runs. The time evolution of $\Omega_{i}$ (for $i=\,$GW, K, or M) integrated over all wave numbers is shown in Fig. 5 for ini1–3, hel1–2, and ac1. In all these cases, the GW energy density saturates at a value ${\Omega}_{\rm GW}^{\rm sat}$ shortly after the sourcing energy density has reached its maximum value $\Omega_{\rm M,K}^{\rm max}$; see Table LABEL:Tsummary.

In Fig. 6, we observe that runs hel1–3, noh1, and ac1 present steeper GW spectra at high frequencies than in runs ini1–3. The monochromatic forcing produces a spike in ${\Omega}_{\rm M}(f)$ at $f_{\ast}$ and a sharp drop in ${\Omega}_{\rm GW}(f)$ and $h_{\rm c}(f)$ beyond $f_{\rm GW}\approx 2f_{\ast}$, for the magnetic runs. For the acoustic runs, we observe a smooth bump on the spectra of ${\Omega}_{\rm GW}(f)$ and $h_{\rm c}(f)$ near the spectral peak $f_{\ast}$. Again, the spectra have the same low-frequency tail, which underlines its universal nature. Also, for a given input energy, $\Omega_{\rm M,K}^{\rm max}$, we obtain larger values of ${\Omega}_{\rm GW}^{\rm sat}$ for acoustic than for vortical turbulence. This case was already studied in Refs. HHRW15 ; HHRW17 ; NSS18 . These features could not be captured by previous analytical estimates, and the power laws in the inertial range also seem to be affected by how the turbulent fields are driven at initial times.

For a given type of initial condition and stirring scale, the final energy density in GWs has the expected quadratic dependence on the source energy density to a very good approximation as shown in Fig. 7. The efficiency of GW production varies significantly with the type of initial conditions; for the same total source energy, the cases with forced acoustic compression lead to a factor of around $200$ more GW energy than those with a sudden magnetic field (ini1–3), while the cases with forced nonhelical magnetic fields are $10$ times more efficient than the latter. We also observe that nonhelical forcing fields are about a factor of $1.6$ more efficient than helical magnetic fields. The detailed reasons behind these significant variations in efficiency are unclear, but they imply that accurate predictions of GW production from cosmological phase transitions will require a detailed model of how latent heat is converted to plasma and magnetic field energies. The comparison of efficiency in GW generation between acoustic and rotational turbulence is a subject of further investigation.

IV.3 LISA Detectability

The projected sensitivity curve for the LISA space mission was plotted in Figs. 4 and 6 along with GW spectra from our runs. The plotted sensitivity assumes a mission of four years RCC18 ; SC19 , and corresponds to the power law sensitivity obtained assuming a threshold signal-to-noise ratio (SNR) of 10 Caprini:2019pxz . The cases ini2–3, each with a turbulent energy input of around 1% of the total radiation energy density, produce signals with a SNR below 10, while ini1, with a turbulent energy input of around 10%, could be detectable. An energy input of about 8% is required to obtain a GW spectrum with a SNR above 10 for runs with an initial helical magnetic field when $N=100$, and 3% when $N=10$. On the other hand, the runs with forced magnetic fields would be detectable for an energy input of approximately 4%, according to our results. Acoustic forced turbulence has been shown to be the more efficient case considered, even though it leads to a GW spectral peak closer to the forcing peak, which slightly reduces the prospects of detection for $T_{\ast}=100\,{\rm GeV}$ and $N=100$. An energy input of around 0.5% would be enough in this case for the GW signal to have a SNR above 10.

V Conclusions

In the early universe, hydrodynamic and MHD turbulence can be an efficient source of GWs. Our direct numerical simulations have shown that the GW energy produced by the turbulence depends quadratically on the energy of the turbulence at the time turbulence is strongest. In the inertial range of the turbulence, the slope of the GW spectrum is steeper than the slope of the magnetic energy spectrum by a factor of $k^{2}$. For a magnetic energy spectrum of Kolmogorov type of the form ${\Omega}_{\rm M}\sim k^{-2/3}$, the GW energy spectrum is of the form ${\Omega}_{\rm GW}\sim k^{-8/3}$. In the subinertial range, however, where the magnetic energy spectrum is expected to be proportional to $k^{5}$, the GW energy spectrum is not proportional to $k^{3}$, as naively expected, but proportional to $k$.

The shallow subinertial range spectrum for the GW energy also implies a novel $f^{-1/2}$ low-frequency spectrum for $h_{\rm c}(f)$. This would enhance the detectability of such a signal compared to the $f^{1/2}$ spectrum obtained from earlier analytic models typically assumed in recent analyses as, e.g., in Ref. croon18 .

Comparing vortical MHD turbulence with irrotational turbulence driven by spherical expansion waves, we find that at similar turbulent energies, irrotational turbulence appears to drive GW energy more efficiently than vortical MHD turbulence. This may be connected with the temporal correlations of the turbulence. Depending on the specific dynamical evolution during the symmetry breaking process, the GWs produced by primordial turbulence may be detectable with LISA when the fraction of radiation energy converted into turbulent energy exceeds a value between 0.5% and 10%. In addition to the total GW energy density, the spectral shape is also affected by the dynamical evolution of the magnetic and/or velocity fields during the phase transition. The specific features around the spectral peak and the power law in the inertial range require numerical simulations to be accurately described.

Now scheduled for launch in the mid 2030s, LISA may provide crucial insight into fundamental physics during the first picoseconds of cosmic evolution.

Data availability—The source code used for the simulations of this study, the Pencil Code, is freely available from Ref. PC . The simulation setup and the corresponding data are freely available from Ref. DATA .

Acknowledgements.

Support through the NSF Astrophysics and Astronomy Grant Program (Grants No. 1615940 and No. 1615100) and the Shota Rustaveli NSF (Georgia) (Grant No. FR/18-1462) are gratefully acknowledged. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm.

Appendix: Spectrum of the source

We compare in Fig. 8 the shell-integrated spectrum of the magnetic field, $E_{\rm M}(k)={\Omega}_{\rm M}(k)/k$ defined such that $\int_{0}^{\infty}E_{\rm M}(k)\,dk=\int_{-\infty}^{\infty}{\Omega}_{\rm M}(k)\,d\ln k$, with the spectrum of the stress tensor $T_{ij}$. The latter spectrum $E_{\rm T}(k)={\Omega}_{\rm T}(k)/k$ is defined such that $\int_{0}^{\infty}E_{\rm T}(k)\,dk=\langle T_{ij}T_{ij}\rangle/2$. Note that, in the absence of fluid motions, this corresponds to the spectrum of the squared magnetic field, whose integral over all wave numbers gives \langle(\mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}^{2})^{2}\rangle/2 instead of just \langle\mathchoice{\mbox{\boldmath\displaystyle B}}{\mbox{\boldmath\textstyle B}}{\mbox{\boldmath\scriptstyle B}}{\mbox{\boldmath\scriptscriptstyle B}}{}^{2}\rangle/2.

In the $k$ range where $E_{\rm M}(k)$ has a Batchelor ($k^{4}$) spectrum, the stress spectrum $E_{\rm T}(k)$ is white noise ($k^{2}$). The $k^{4}$ spectrum is caused by a white noise ($k^{2}$) spectrum of the vector potential. Recent work BB19 shows that the spectrum of the stress $E_{\rm T}(k)$ can never be steeper than that of white noise, and that the peak of the stress spectrum shifts to $2k_{\ast}$, being $k_{\ast}$ the position of the spectral peak of the magnetic field, as it is observed in Fig. 8. In the inertial range, we observe both spectra $E_{\rm M}(k)$ and $E_{\rm T}(k)$ to possess the same Kolmogorov scaling $k^{-5/3}$. These results have also been confirmed analytically by calculating the spectrum of the stress as the autocorrelation function of the two turbulence spectra under the assumption that the underlying fields are Gaussian distributed BB19 . As inferred from Eq. (4), the relevant spectrum related to the GW energy is the stress spectrum, instead of the magnetic spectrum.