Prepared for submission to JCAP

Polarization of gravitational waves from helical MHD turbulent sources

Alberto Roper Pol^{a,b,1} Sayan Mandal^{c,b} Axel Brandenburg^{d,e,f,b} Tina Kahniashvili^{f,b,g,h}

aUniversité de Paris, CNRS, Astrophysique et Cosmologie, Paris, F-75013, France bFaculty of Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili Ave, Tbilisi, GE-0194, Georgia cPhysics and Astronomy Department, Stony Brook University, Stony Brook, New York 11794, USA dNordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden eDepartment of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden fMcWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, PA 15213, USA gAbastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia hDepartment of Physics, Laurentian University, Ramsey Lake Road, Sudbury, ON P3E 2C, Canada

E-mail: roperpol@apc.in2p3.fr; Revision: 1.59

Abstract. We use direct numerical simulations of decaying primordial hydromagnetic turbulence with helicity to compute the resulting production of gravitational waves (GWs) and their degree of polarization. The turbulence sourcing field, which can be magnetic or kinetic, is either given initially, or it is produced by driving it with an electromotive force applied during a short initial time interval. In both types of simulations, we find a clear dependence of the degree of circular polarization of the resulting GWs on the fractional magnetic/kinetic helicity. Furthermore, the spectral peak is higher in the helical case, which facilitates its observational detection with the planned Laser Interferometer Space Antenna (LISA). In both cases, the low frequency tail is shallower than what was expected based on earlier analytical work. For frequencies above the peak, the two types of simulations show opposite trends in that the spectral energy is higher (lower) in the helical case than in the nonhelical one if an initial spectrum is given (driven). The degree of polarization approaches zero at frequencies below the peak, but stays finite at higher frequencies if the initial magnetic field is given to be fully helical all wavenumbers, and it approaches zero if the initial field is driven by a fully helical monochromatic forcing function.

Contents

1 Introduction 1

2 The model 2

3 Results 4

4 Possibility of detection with LISA 7

5 Conclusions 7

A Analytical model for polarization 9

B Beltrami field 9

1 Introduction

Primordial magnetic fields produced during phase transitions can drive stochastic background gravitational waves (GWs) in a broad range of frequencies, from nHz to Hz, depending on the energy scale of the phase transition (T ∼ 100 GeV for the electroweak, and T ∼ 100 MeV for the QCD phase transition), and the characteristic scale of the magnetic field k_{∗}, which is taken to be a fraction of the Hubble scale at the time of magnetic field generation. These signals could in future be observable with gravitational wave detectors, e.g., pulsar timing arrays (PTA), or the planned Laser Interferometer Space Antenna (LISA). These magnetic fields could have survived until the present time [14] and are thus a leading candidate for explaining the lower limits on the strength of intergalactic magnetic fields inferred from the non-observation of cascade photos from TeV photons of blazars interacting with the cosmic background light [1]. Such fields are expected to be turbulent and partially helical.

Helical magnetic fields produce circularly polarized GWs [16]. However, the degree of polarization of GWs has been a matter of uncertainty owing to the approximations made in the analytical calculations available to date. The work of Ref. [11] suggested that the circular polarization of GWs would be about 80% for a maximally helical magnetic field, by assuming Kolmogorov-type turbulence, with spectral index of −11/3 for the magnetic power spectrum and −14/3 for the helicity power spectrum, following Ref. [12]. The peak is observed at twice the magnetic spectral peak. For smaller values of the fractional magnetic helicity, the maximum degree of polarization also diminishes. Ref. [16] showed that the maximum circular polarization depends on the relation between the magnetic field and magnetic helicity power spectra, and it can actually reach nearly 100% in the case when the spectral indices for both are equal, which we call a Moiseev-type spectrum, following Ref. [17]. On the other hand, this type of spectrum generates a degree of polarization that does not decay to zero for large wave numbers, but it stays constant and is approximately equal to the fractional helicity. More recently, the two types of spectra mentioned above have been studied in Ref. [13] in the context of detectability with LISA. We show that the circular degree of polarization computed from direct numerical simulations does not follow any of the previous analytical models, even though it retains some similarities, depending on the how the magnetic field is generated/driven at the initial time.

In two recent papers [3, 4], the authors have described the implementation of a GW solver into the Pencil Code [5] and have presented numerical computations for primordial hydrodynamic and hydromagnetic turbulence from phase transitions, and their subsequent generation of a stochastic gravitational wave background, also computed numerically. Some of those solutions were fully helical, but helicity was not a focus of their studies. In Ref. [6] the degree of circular polarization for kinetically and magnetically forced turbulence was presented. We extend the study in the present work by including the effect of fractional helicity to the magnetic field power spectrum, and the variation of the polarization for different assumptions of magnetic field production.

We present here two types of simulations, similar to the two types of hydromagnetic simulations presented in Ref. [4]: one where a primordial magnetic field was assumed given as initial condition, and one where a magnetic field was generated by an electromotive force E(x, t) that depends on time t and position x. This driving was then applied during a short time interval, and then switched off, i.e, set to zero, E = 0, such that turbulence decays for later times. By using suitably scaled variables and conformal time, the governing equations describing the evolution of GWs and turbulent magnetic fields in an expanding universe in the radiation era can be brought into a form that is best suited for numerical simulations [3].

We begin by summarizing our approach and the equations solved. We then present the magnetic and GW energy spectra, the degree of polarization in relation to the magnetic helicity, and compare with the sensitivity limits for LISA. We conclude by highlighting unique properties of the GW spectra and their polarization properties.

Throughout this work, electromagnetic quantities are expressed in Lorentz–Heaviside units where the vacuum permeability is unity. Einstein index notation is used so summation is assumed over repeated indices. Latin indices i and j refer to spatial coordinates 1 to 3. The linear polarization modes are indicated by λ = +, ×.

2 The model

We use the linear polarization modes +, × to describe the two gauge independent components of the tensor-mode perturbations, h̃_{ij}(k) = h̃_{+}(k)e_{ij}^{+}(k) + h̃_{×}(k)e_{ij}^{×}(k), where the tilde indicates that this decomposition is performed in Fourier space. The linear polarization basis tensors are

e_{ij}^{+}(k) = e_{i}^{1}e_{j}^{1} − e_{i}^{2}e_{j}^{2}, e_{ij}^{×}(k) = e_{i}^{1}e_{j}^{2} + e_{i}^{2}e_{j}^{1},  (2.1)

where e^{1} and e^{2} form a basis with the units vectors k/|k| [3]. We solve the non-dimensional GW equation in the radiation era for the scaled strains h̃_{λ}, where λ = +, ×, using conformal time, normalized to unity at the initial time of magnetic energy generation t_{∗}, and comoving wave vector, normalized by the corresponding scale 1/(ct_{∗}) [3, 4]

(∂_{t}^{2} + k^{2}) h̃_{λ}(k, t) = \frac{6}{t}T̃_{λ}^{TT}(k, t),  (2.2)

where T̃_{λ}^{TT} is the comoving stress tensor, projected into the traceless and transverse (TT) gauge, described by the linear polarization modes, and normalized by the energy density at t_{∗}. The scaled strains are tensor-mode perturbations over the Friedmann-Lemaître-Robertson-Walker metric tensor. Hence, the metric tensor is such that the line element is ds^{2} = a^{2}(−dt^{2} + δ_{ij} + h_{ij}/a) [3, 4]. During radiation-dominated epoch, the equation of state is p = \frac{4}{3}ρ, which leads to a linear evolution of the scale factor a with t.

The stress is composed of magnetic and kinetic contributions. In physical space it is

T_{ij}(x) = \frac{4}{3} \frac{ρu_{i}u_{j}}{1 − u^{2}} − B_{i}B_{j} + (\frac{ρ}{3} + \frac{B^{2}}{2} ) δ_{ij},  (2.3)

where u is the plasma velocity, ρ is the energy density, and B is the magnetic field. The total enthalpy is w = p + ρ = \frac{4}{3}ρ. Since T_{ij} refers to comoving and normalized stress tensor, the MHD fields ρ, u, B are accordingly normalized and comoving.

The non-dimensional and comoving MHD equations for an ultrarelativistic gas in a flat expanding universe in the radiation-dominated era after the electroweak phase transition are given by [3, 4, 7, 8, 14]

\frac{∂ ln ρ}{∂t} = −\frac{4}{3} (∇ · u + u · ∇ ln ρ) + \frac{1}{ρ} [u · (J × B) + ηJ^{2}] ,
\frac{∂u}{∂t} = −u · ∇u + \frac{u}{3} (∇ · u + u · ∇ ln ρ) + \frac{2}{ρ}∇ · (ρνS)

−\frac{1}{4}∇ ln ρ − \frac{u}{ρ} [u · (J × B) + ηJ^{2}] + \frac{3}{4ρ}J × B,

\frac{∂B}{∂t} = ∇ × (u × B − ηJ + E), J = ∇ × B,

where S_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) − \frac{1}{3}δ_{ij}∇ · u are the components of the rate-of-strain tensor with commas denoting partial derivatives, J is the current density, ν is the kinematic viscosity, and η is the magnetic diffusivity.

The electromotive force, E, is used to model the magnetic field generation with fractional magnetic helicity. In some of our runs we will omit this term, and we will start with a fully developed turbulent spectrum of the magnetic field as our initial condition. On the other hand, we also consider magnetic fields generated dynamically using E. Specifically, in the latter case, we write

E(x, t) = Re{N f̃(k(t)) exp[ik · x + iφ]},  (2.4)

where the wave vector k(t) and the phase φ(t) change randomly from one time step to the next. This forcing function is therefore white noise in time and consists of plane waves with average wavenumber k_{∗} such that |k| lies in an interval k_{∗} − δk/2 ≤ |k| < k_{∗} + δk/2 of width δk. Here, N is a normalization factor. The Fourier amplitudes of the forcing are

f̃_{i} = (δ_{ij} − iσϵ_{ijl}k_{l}/k) f̃_{j}^{(0)}/\sqrt{1 + σ^{2}},  (2.5)

where f̃^{(0)}(k) = (k × e)/[k^{2} − (k · e)^{2}]^{1/2} is a nonhelical forcing function. Here, e is an arbitrary unit vector that is not aligned with k. Note that |f|^{2} = 1. The values σ = 0 and σ = 1 correspond to nonhelical and maximally helical cases, but we also consider intermediate values, which correspond to fractional helicities. The forcing is only enabled during an arbitrarily short time interval 1 ≤ t ≤ t_{max}. We chose t_{max} = 1.1. The resulting fractional magnetic helicity k_{∗}⟨A·B⟩/⟨B^{2}⟩ is given approximately by 2σ/(1+σ^{2}), being A the magnetic vector potential, and angle brackets denoting ensemble average over stochastic realizations.

When the magnetic field is given as initial condition, we generate a random threedimensional vector field in Fourier space,

B_{i}(k) = B_{0} (P_{ij}(k) − iσϵ_{ijl}k_{l}/k) g_{j}(k) g_{0}(k),  (2.6)

where g_{j}(k) is the Fourier transform of a δ-correlated vector field in three dimensions with Gaussian fluctuations, i.e., g_{i}(x)g_{j}(x′) = δ_{ij}δ^{3}(x − x′), k_{∗} is now identified with the initial wavenumber of the energy-carrying eddies, and g_{0}(k) determines the spectral shape [14].

g_{0}(k) = \frac{k_{∗}^{−3/2}(k/k_{∗})^{α/2−1}}{[1 + (k/k_{∗})^{2(α+5/3)}]^{1/4}} .  (2.7)

The power spectrum of magnetic fields is computed from \frac{1}{2}⟨B^{2}⟩. Hence, it is proportional to k^{2}g_{0}^{2}(k) in Fourier space, which correspond to a spectral index α in the low wave number limit (subinertial range) and Kolmogorov-type spectral slope −5/3 in the high wave number range. Here, the k_{∗}^{−3/2} prefactor ensures that the resulting magnetic energy is independent of the value of k_{∗}.

The equal time correlation function of the magnetic field, assuming statistical homogeneity and isotropy, and a Gaussian-distribution in space, is

⟨B_{i}^{∗}(k, t)B_{j}(k′, t)⟩ = (2π)^{6}δ^{3}(k − k′)[P_{ij} \frac{E_{M}(k, t)}{4πk^{2}} + iϵ_{ijl}k_{l} \frac{H_{M}(k, t)}{8πk^{2}}],  (2.8)

where P_{ij} = δ_{ij} −k_{i}k_{j}/k^{2}, and E_{M}(k, t), H_{M}(k, t) are the magnetic and helical power spectra, respectively. We work here with magnetic energy spectra per linear wave number interval.

The GW energy density (in physical units) is

E_{GW}(t) = \frac{c^{2}}{32πG}⟨ḣ_{ij}^{phys}ḣ_{ij}^{phys}⟩,  (2.9)

where h_{ij}^{phys} = h_{ij}/a are the physical strains, and a dot represents derivative with respect to physical time. In terms of the normalized and comoving units used in Eq. (2.2), the ratio of comoving GW energy density to critical energy density Ω_{GW} = E_{GW}/E_{crit}^{0} is

a^{4}Ω_{GW}(t) = \frac{1}{12}(\frac{H_{∗}}{H_{0}} )^{2}⟨\frac{∂h_{ij}}{∂t} \frac{∂h_{ij}}{∂t} + \frac{1}{t^{2}} h_{ij}h_{ij} − \frac{2}{t}h_{ij} \frac{∂h_{ij}}{∂t} ⟩,  (2.10)

where E_{crit}^{0} = 3H_{0}^{2}c^{2}/(8πG), with H_{0} = 100 h_{0} km s^{−1} Mpc^{−1} ≈ 3.241 × 10^{−18} h_{0} s^{−1} being the Hubble rate at the present time, and h_{0} takes into account the uncertainties in its exact value. To simplify the notation, we will use prime to denote derivative with respect to t.

The equal time correlation function for the strains and their derivatives (again assuming the statistical fields to be homogeneous, isotropic and Gaussian) can be expressed as

⟨h′_{ij}(k, t)h′_{ij}(k′, t)⟩ = \frac{1}{4}[M_{ijlm}S_{h′}(k, t) + iA_{ijlm}A_{h′}(k, t)],
⟨h_{ij}(k, t)h_{ij}(k′, t)⟩ = \frac{1}{4}[M_{ijlm}S_{h}(k, t) + iA_{ijlm}A_{h}(k, t)],
⟨h_{ij}(k, t)h′_{ij}(k′, t)⟩ = ⟨h′_{ij}(k, t)h_{ij}(k′, t)⟩ = \frac{1}{4}[M_{ijlm}S_{mix}(k, t) + iA_{ijlm}A_{mix}(k, t)], (2.11)

3 Results

We have computed solutions for a range of values of σ, using both given and driven initial fields. In Fig. 1 we compare energy and helicity spectra of the magnetic field and the resulting GWs for σ = ±1 and ±0.01. The spectra are averaged over the time interval 1.1 ≤ t ≤ 1.2,

Figure 1. Magnetic and GW energy and helicity spectra.

which is just after the maximum magnetic energy has been reached and the GW energy begins to fluctuate around an approximately statistically steady state.

We see that in all cases, the causal k^{4} magnetic energy spectrum [9] is established. This implies a flat spectrum for the GW energy [4, 15]. For k > k_{∗}, the slope of the magnetic energy spectrum is a little steeper This is because of the finite time driving during the rather short time interval, 1 ≤ t ≤ 1.1. The k^{−5/3} Kolmogorov-type spectrum would emerge if the driving was continued over a longer time interval. than the expected k^{−5/3} spectrum.

In Fig. 1, the magnetic helicity spectra are multiplied by k/2, which allows us to see that the spectral realizability condition [10],

|kH_{M}(k, t)/2| ≤ E_{M}(k, t),  (3.1)

is nearly saturated at k = k_{∗}, so the magnetic energy is here nearly fully helical for σ = ±1; see the upper two panels of Fig. 1. For σ = ±0.01 (lower two panels) there is a clear separation between E_{M}(k, t) and kH_{M}(k, t)/2 and the sign of the magnetic helicity tends to fluctuate noticeably. We also see a systematic sign flip at higher k both in H_{M} and in H_{GW}. Such sign flips are typical of decaying helical turbulence are a consequence of the fact that the fractional magnetic helicity is there already extremely small.

In Fig. 2 we plot the GW polarizations

P(k) = \frac{Ξ_{GW}(k)}{Ω_{GW}(k)},  (3.2)

Figure 2. (a) given initial field with different fractional helicity. (b) driven initial field with different fractional helicity. + analytical formula.

where

Ξ_{GW}(k) = \frac{k}{3H_{∗}^{2}} ∫_{4π} Im (˙h̃_{+}^{phys} ˙h̃_{×}^{phys,∗}) k^{2} dΩ_{k} ,  (3.3)

with ∫ Ξ_{GW}(k) d ln k = Ξ_{GW} being the total normalized helical energy density. We also define the fractional magnetic helicity P_{M}(k) = kH_{M}(k)/2E_{M}(k).

We see that P_{GW}(k) reaches ±1 at k ≈ 2k_{∗} ≈ 1200 when σ = ±1. This is more than what was found in Ref. [11]. This discrepancy can probably be explained by the use of the approximations made in the analytic calculations. Toward smaller k, there are systematic fluctuations around P_{GW}(k) = 0, especially for k = 200.

We also plot the total polarization integrated over all wave numbers, i.e.,

P̃_{GW} = \frac{Ξ_{GW}}{Ω_{GW}} , P̃_{M} = \frac{H_{M}}{2ξ_{M}E_{M}} .  (3.4)

The result is shown in Fig. 3. We see that the numerical data can well be fitted with the curve

P̃_{GW} ≈ P̃_{0} tan P̃_{M},  (3.5)

Figure 3. GW polarization versus magnetic polarization.

where P̃_{0} ≈ 0.6 for the driven case and P̃_{0} ≈ 0.76 for the case with an initial magnetic field. This relation between P̃_{GW} and P̃_{M} is unexpected and the agreement with the actual data points is remarkable. On the other hand, for a Beltrami field, which is monochromatic and defined only for a single value of k, we have an exact analytical solution for both P̃_{GW} and P̃_{M}; see Appendix A. The resulting dependence of P̃_{GW} on P̃_{M} has a rather different shape; compare with the dotted line for the Beltrami field. This is surprising, but it may be related to the smaller fractional polarization away from the peak near k_{∗}; see Fig. 2. This feature becomes more pronounced for smaller values of σ. On the other hand, close to the maximally helical fields, P̃_{GW} shows a much more dramatic increase with P̃_{M}.

4 Possibility of detection with LISA

To assess the observational prospects of detecting GWs for different values of σ, we now plot Ω_{GW} for σ = 1, 0.7, 0.5, 0.3, and 0.1; see Fig. 4. The result for σ = 1 was already shown in the top left panel of Fig. 1. The spectra are similar, but we now also see that for smaller values of σ, the jump in Ω_{GW}(f) near the peak is less pronounced, so for larger frequencies, i.e., to the right of the peak, Ω_{GW}(f) increases (decreases) for smaller (larger) values of σ.

For smaller frequencies, to the left of the peak, which corresponds to the subinertial range, we have the aforementioned shallow spectrum with Ω_{GW}(f) ∝ f, which is approximately independent of the value of σ. Again, positive (negative) values are indicated by a red (blue) symbols.

Finally, we plot in Fig. 5 the energy spectra together with the helicity spectra Σ_{GW}. In this logarithmic representation, the fractional GW polarization appears still rather noticeable for σ = ±0.1. We see that for fractional magnetic helicity, the GW power decreases and the values of Ξ_{GW}(f) are now close to the LISA sensitivity limit, Ξ_{GW}^{LISA}(f). On the other hand, the polarization signal is still rather clean in the sense that the amount of statistical fluctuations in the sign of Ξ_{GW}^{LISA}(f) is weak.

5 Conclusions

Our numerical simulations have confirmed that there is a direct correspondence between the sign of magnetic helicity and the sign of circular polarization of the GWs produced

Figure 4. Diagnostic diagram of h_{0}^{2}Ω_{GW}(f) for σ = 1 (red), 0.7 (orange), 0.5 (green), 0.3 (blue), 0.1 (black) for the driven case (left) and the case with initial condition (right). The LISA sensitivity curve is shown as a green dash-dotted line and slopes proportional to f^{−4} and f^{−5}, as well as f^{−8/3} are shown for orientation.

Figure 5. Diagnostic diagram of h_{0}^{2}Ω_{GW}(f) and h_{0}^{2}Ξ_{GW}(f) for given initial magnetic fields with (a) σ = ±1 and (b) ±0.1, and initially driven magnetic fields with (c) σ = ±1 and (d) ±0.1.

from the resulting magnetic stress. Our work has also revealed a new functional form of the dependence of circular polarization of GWs on the fractional magnetic helicity. The dependence is shallower for weakly helical magnetic fields and steeper for nearly fully helical fields. The precise form of this curve can be important if one wants to infer the magnetic helicity from circular polarization measurements of GWs. Cosmological magnetic fields may well be close to fully helical, because the magnetic helicity is a conserved quantity while the magnetic energy decays, so the ratio always increases until it reaches nearly 100%.

The subinertial range of the GW spectrum is hardly affected by the presence or absence of magnetic helicity. However, the peak of the GW spectrum is higher in the case with magnetic helicity. This is another reason why magnetic helicity is beneficial for enhancing the detection prospects of GWs, provided LISA is indeed able to measure the GW energy at the peak frequency. On the other hand, at larger wave numbers, corresponding to the inertial range of the underlying turbulence, the trend is the other way around, i.e., the spectral GW energy is lower in the presence of helicity and the slope steeper.

A Analytical model for polarization

B Beltrami field

Simple magnetic field with partial helicity

B(x, t) = B_{0}Θ(t) ( \begin{matrix}0 \\ σ sin k_{0}x \\ cos k_{0}x \end{matrix} ) ,  (B.1)

where σ ∈ [0, 1] is a parameter to modify the helicity of the magnetic field, k_{0} ∈ R is a characteristic wavenumber, and Ω_{M} = B_{0}^{2}/2 is the amplitude. The Heaviside function Θ(t) is included to indicate that this field is zero for t ≤ 0.

In Fourier space this is

B̃(k, t) = \frac{1}{2}B_{0}Θ(t) ( \begin{matrix}0 \\ iσ[δ^{3}(k + k_{0}) − δ^{3}(k − k_{0})] \\ δ^{3}(k + k_{0}) + δ^{3}(k − k_{0}) \end{matrix} ) ,  (B.2)

where k_{0} = (k_{0}, 0, 0), and δ^{3}(k) is the Dirac’s delta function in three dimensions. We have used the convention

B̃(k, t) = ∫ B(x, t)e^{−ik·x}dx, B(x, t) = \frac{1}{(2π)^{3}} ∫ B̃(k, t)e^{−ik·x}dk.  (B.3)

The magnetic spectrum is

E_{M}(k, t) = \frac{1}{2} ∫_{Ω_{1}} B̃(k, t) · B̃^{∗}(k, t) dΩ_{1} = \frac{Ω_{M}}{2} Θ(t)(1 + σ^{2})δ(k − |k_{0}|),  (B.4)

where Ω_{1} = 2 is the one-dimensional solid angle, and δ(k) is the one-dimensional Dirac’s delta function.

The total magnetic energy density is

E_{M}(t) = ∫_{0}^{∞} E_{M}(k, t) dk = \frac{Ω_{M}}{2} Θ(t)(1 + σ^{2}).  (B.5)

The magnetic helicity spectrum is

H_{M}(k, t) = ∫_{Ω_{1}} B̃(k, t) · (ik × B̃(k, t))^{∗} dΩ_{1} = 2Ω_{M}Θ(t)|k_{0}|σδ(k − |k_{0}|).  (B.6)

The total magnetic helicity density is

H_{M}(t) = ∫_{0}^{∞} H_{M}(k, t) dk = 2Ω_{M}Θ(t)|k_{0}|σ.  (B.7)

The magnetic spectral peak is

k_{M} = \frac{E_{M}(t)}{∫_{0}^{∞} k^{−1}E_{M}(k, t) dk} = |k_{0}|.  (B.8)

The fractional magnetic helicity is

ε_{h} = \frac{⟨B(x, t) · (∇ × B(x, t)⟩}{k_{M}⟨B^{2}(x, t⟩} = \frac{H_{M}(t)}{2k_{M}E_{M}(t)} = \frac{2σ}{1 + σ^{2}} ,  (B.9)

which reduces to +1, 0, and −1, for σ = +1, 0, and −1, respectively.

The stress tensor for magnetic fields is

T_{ij}(x, t) = −B_{i}B_{j} + \frac{1}{2}δ_{ij}B^{2},  (B.10)

where

B^{2}(x, t) = Ω_{M}Θ(t) [(1 + σ^{2}) + cos 2k_{0}x(1 − σ^{2})] ,  (B.11)

and

− B_{i}B_{j}(x, t) = −Ω_{M}Θ(t) ( \begin{matrix}0 & 0 & 0 \\ 0 & σ^{2}(1 − cos 2k_{0}x) & σ sin 2k_{0}x \\ 0 & σ sin 2k_{0}x & 1 + cos 2k_{0}x \end{matrix} ) .  (B.12)

Combining both we obtain the stress tensor

T_{ij}(x, t) = \frac{1}{2}Ω_{M}Θ(t) ( \begin{matrix}T_{11}(x) & 0 & 0 \\ 0 & T_{22}(x) & T_{12}(x) \\ 0 & T_{12}(x) & −T_{22}(x) \end{matrix} ) ,  (B.13)

where

T_{11}(x) = (1 + σ^{2}) + (1 − σ^{2}) cos 2k_{0}x,
T_{22}(x) = (1 − σ^{2}) + (1 + σ^{2}) cos 2k_{0}x,
T_{12}(x) = 2σ sin 2k_{0}x.  (B.14)

Since GWs are produced by linear perturbations over the metric tensor, and the stress tensor components are also perturbations over background fields, constant modes (k = 0) do not produce GWs. For this reason, we can neglect the terms that are homogeneous in space. This can be done by computing T̃_{ij}^{k}(k, t) = T̃_{ij}(k, t)(1 − δ^{3}(k)). The stress obtained is the same as in Eq. (B.13), replacing T_{11}(x) and T_{22}(x) by

T_{11}^{k}(x) = (1 − σ^{2}) cos 2k_{0}x, T_{22}^{k}(x) = (1 + σ^{2}) cos 2k_{0}x.  (B.15)

In Fourier space, the resulting stress tensor is

T̃_{ij}^{k}(k, t) = \frac{1}{2}Ω_{M}Θ(t) ( \begin{matrix}T̃_{11}^{k}(k) & 0 & 0 \\ 0 & T̃_{22}^{k}(k) & T̃_{12}(k) \\ 0 & T̃_{12}(k) & −T̃_{22}^{k}(k) \end{matrix} ) ,  (B.16)

where

T̃_{11}^{k}(k) = \frac{1}{2}(1 − σ^{2})[δ^{3}(k + 2k_{0}) + δ^{3}(k − 2k_{0})],

T̃_{22}^{k}(k) = \frac{1}{2}(1 + σ^{2})[δ^{3}(k + 2k_{0}) + δ^{3}(k − 2k_{0})],

T̃_{12}(k) = iσ[δ^{3}(k + 2k_{0}) − δ^{3}(k − 2k_{0})].  (B.17)

The transverse and traceless projection of the stress tensor is

T̃_{ij}^{TT}(k, t) = (P_{il}P_{jm} − \frac{1}{2}P_{ij}P_{lm})T̃_{lm}^{k}(k, t),  (B.18)

where P_{ij} = δ_{ij} − k̂_{i}k̂_{j} is the projection operator being k̂ = k/\sqrt{k^{2}}.

The only wave vectors that leads to non-zero stress in Fourier space are k = 2k_{0} = (±2k_{0}, 0, 0). The projection operator for these wave vectors is P_{ij}(±2k_{0}) = δ_{ij} − δ_{i1}δ_{j1}.

References

[1] A. Neronov and I. Vovk, Science 328, 73 (2010).

[2] A. S. Monin and A. M. Yaglom. Statistical Fluid Mechanics Mechanics of Turbulence, vol. 2. MIT press, Cambridge, 1975.

[3] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky, and S. Mandal, Geophys. Astrophys. Fluid Dynam. 114, 130 (2020). The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence

[4] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili, and A. Kosowsky, , Phys. Rev. D 102, 083512 (2020).

[5] The Pencil Code. https://github.com/pencil-code. DOI:10.5281/zenodo.2315093

[6] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal, and A. Roper Pol, , Phys. Rev. Lett. submitted, arXiv:2011.05556 (2020).

[7] A. Brandenburg, K. Enqvist, and P. Olesen, Phys. Rev. D 54, 1291 (1996).

[8] R. Durrer and A. Neronov, Astron. Astrophys. Rev. 21, 62 (2013).

[9] R. Durrer and C. Caprini, J. Cosmol. Astropart. Phys. 0311, 010 (2003).

[10] H. K. Moffatt, Magnetic Field Generation in Electrically Conducting Fluids. Cambridge: Cambridge Univ. Press (1978).

[11] L. Kisslinger and T. Kahniashvili, Phys. Rev. D 92, 043006 (2015). Polarized gravitational waves from cosmological phase transitions

[12] R. H. Kraichnan, Phys. Fluids 8, Inertial-range spectrum of hydromagnetic turbulence (1965).

[13] J. Ellis, M. Fairbairn, M. Lewicki, V. Vaskonen, and A. Wickens, J. Cosmol. Astropart. Phys. 10 (2020) 032.

[14] A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze, and T. Vachaspati, Phys. Rev. D 96, 123528 (2017).

[15] A. Brandenburg and S. Boldyrev, Astrophys. J. 892, 80 (2020). The turbulent stress spectrum in the inertial and subinertial ranges

[16] T. Kahniashvili, G. Gogoberidze, and B. Ratra, Phys. Rev. Lett. 95, 151301 (2005).

[17] S. S. Moiseev and O. G. Chetiani, Sov. Phys. JETP 83, 192 (1996).