# The scalar, vector, and tensor modes in gravitational wave turbulence
simulations

Axel Brandenburg1,2,3,4, Grigol Gogoberidze4, Tina Kahniashvili3,4,5, Sayan
Mandal6,4111Corresponding author; all authors listed alphabetically., Alberto
Roper Pol7,4, Nakul Shenoy8 1 Nordita, KTH Royal Institute of Technology and
Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden  2
Department of Astronomy, AlbaNova University Center, Stockholm University,
10691 Stockholm, Sweden  3 Department of Physics and McWilliams Center for
Cosmology, Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, PA 15213,
USA  4 Faculty of Natural Sciences and Medicine, Ilia State University, 3-5
Cholokashvili Street, 0194 Tbilisi, Georgia  5 Department of Physics,
Laurentian University, Sudbury, ON P3E 2C, Canada  6 Physics and Astronomy
Department, Stony Brook University, 100 Nicolls Road, Stony Brook, New York
11794, USA  7 Université de Paris, CNRS, Astrophysique et Cosmologie, Paris,
F-75013, France  8 Department of Physics, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>()

###### Abstract

We study the gravitational wave (GW) signal sourced by primordial turbulence
that is assumed to be present at cosmological phase transitions like the
electroweak and quantum chromodynamics phase transitions. We consider various
models of primordial turbulence, such as those with and without helicity,
purely hydrodynamical turbulence induced by fluid motions, and
magnetohydrodynamic turbulence whose energy can be dominated either by kinetic
or magnetic energy, depending on the nature of the turbulence. We also study
circularly polarized GWs generated by parity violating sources such as helical
turbulence. Our ultimate goal is to determine the efficiency of GW production
through different classes of turbulence. We find that the GW energy and strain
tend to be large for acoustic or irrotational turbulence, even though its
tensor mode amplitude is relatively small at most wave numbers. Only at very
small wave numbers is the spectral tensor mode significant, which might
explain the efficient GW production in that case.

Keywords: gravitational waves, magnetohydrodynamics, turbulence, cosmological
phase transitions

##  1 Introduction

Gravitational waves (GWs) can provide a window into the earliest epochs of the
universe when it was opaque to electromagnetic radiation (see Ref. [1],
Chapters: 21–22). Being massless and decoupled from other physical
interactions (weak, strong, and electromagnetic) GWs propagate almost222We
discard the damping due to neutrino free streaming [2, 3, 4, 5, 6, 7, 8] or
from anisotropic stresses [9] freely after their generation preserving their
spectral characteristics and diluting only with the expansion of the universe.
In particular, cosmological GW backgrounds produced by magnetohydrodynamic
(MHD) turbulence conserve the characteristics of the source, i.e., the
spectral peak and the energy density of the produced GW signal can be related
to the characteristic scale and total energy density of the sourcing magnetic
or velocity field at the time of generation (with a precision of about one
order of magnitude for the energy density depending on the duration of the
turbulence sourcing) [10].

If GWs are generated from the early universe at energy scales around the
electroweak (EW) and quantum chromodynamics (QCD) phase transitions are
detected, they will carry information about the fundamental physical processes
at play when the Universe was only $10^{-11}$ and $10^{-5}$ seconds old,
respectively. The Laser Interferometer Space Antenna (LISA), slated to launch
in 2034, will be sensitive to the GWs generated around the
$100\,\mathrm{GeV}-1\,\mathrm{TeV}$ scale (see Refs. [11, 12, 13] for
pioneering works and Ref. [14] for a review and reference therein). It is also
conceivable that the GWs generated at the EW epoch can be detected in atomic
interferometry experiments [15], and those originated at the QCD epoch (or
around) at typical frequencies of $10^{-7}$–$10^{-9}$ Hz in pulsar timing
array (PTA) experiments [16, 17, 18]. Indeed PTAs provide detection of GW
signals in the nHz regime making the GW signal arising at QCD energy scale
accessible. Recently, the NANOGrav collaboration has reported the detection of
a signal, which could be compatible with a stochastic GW background [19]. If
this detection is confirmed in the future, it could be a manifestation of the
stochastic GW background originated from the epoch around the QCD energy scale
[20, 21, 22, 24, 23, 25, 26]. Interestingly, the source of these GWs might be
the presence of primordial magnetic fields at the QCD scale [21]333The
commonly discussed sources of GWs at the nHz regime are mergers of
supermassive black holes [19, 27, 28]..

Stochastic relic GW backgrounds can be generated through different processes
in the early universe, such as amplification of quantum-mechanical
fluctuations during inflation, topological defects, cosmic strings, and
cosmological phase transitions, see Ref. [29] for a review and references
therein. In the case of cosmological phase transitions GWs are generated
either from the collision of bubbles during first order phase transitions,
like the EW or QCD phase transitions (see Refs. [12, 13] for pioneering
works), or from the hydrodynamic and magnetohydrodynamic (MHD) turbulence
produced in association with the generation of primordial magnetic fields
(PMFs) and turbulence at those early times; see Refs. [9, 30, 31] for
pioneering works, and for a recent review, see Sec. 8 of Ref [29] and
references therein.

The stochastic magnetic fields and the turbulent velocity fields of the
primordial plasma contribute to the anisotropic stress tensor that sources the
primordial GWs. In addition there might be parity violating (or CP violating)
processes in the early universe that lead to the generation of helical
magnetic fields, or helical turbulent motions, and thus produce circularly
polarized GWs [32]. In particular, the detection of circularly polarized GWs
will shed light on phenomena of fundamental symmetry breaking in the early
universe, such as parity violation, and potentially can serve as an
explanation of the lepto- and baryogenesis problem (see Refs. [33, 34, 35, 36]
for pioneering works, and Ref. [37] for a recent review).

In this paper, we study the scalar, vector, and tensor mode decomposition of
different turbulence sources with the objective to help understand their
efficiency in producing GWs. More precisely we have considered different
classes of turbulence that may be realized in the early universe. In
particular we have studied kinetically and magnetically dominant helical and
non-helical cases as well as two types of acoustic turbulence. The latter type
consists of interacting sound waves. At large enough Mach numbers, the flow is
dominated by shocks. Notably it was shown that the efficiency of sound waves
as a source generating relic GWs is not only comparable but exceeds that from
bubble collisions; see Ref. [38] for a pioneering work, and Ref. [39] for a
recent study and references therein. However the generation mechanism related
to sound waves requires a first order PT. The ultimate goal of the work
presented here is to determine the physical conditions in the early universe
that reach maximal efficiency of the GW signal production without being
limited by the first order PT consideration. We also study helical vs non-
helical sources behavior, and we confirm that, for distinguishing between
parity conserving and parity violating sources, the measurement of
polarization degrees in the low frequency tail of the spectrum could be used
to identify the inverse cascade (which applies to the helical turbulence case)
versus the slow inverse transfer.

This paper is arranged as follows. In Sec. 2, we briefly review the
decomposition of the MHD stress tensor sourcing the GWs into the scalar,
vector, and tensor components, and provide a background on the modeling of
turbulent sources. This is followed by detailed numerical modeling for several
types of turbulent hydrodynamic and/or MHD flows in Sec. 3. In Sec. 4, we
describe in brief the degree of circular polarization of both the source and
the generated GWs, and present numerical simulations of the degree of
polarization of the source. Finally in Sec. 5, we present our conclusions.
Throughout this work, we set $\hbar=c=k_{B}=1$, and use the metric signature
$(-,+,+,+)$. We also set the permeability of free space to unity, i.e.,
$\mu_{0}=1$, expressing the electromagnetic quantities in Lorentz-Heaviside
units. Repeated indices denote a summation, unless otherwise specified. Greek
indices run over spacetime coordinates, while Latin indices run over spatial
coordinates.

##  2 Formalism

The stress tensor that sources GWs in the early universe receives
contributions from PMFs [9], as well as the turbulent velocity field of the
primordial plasma [40, 41], and the dynamics of phase transitions [11, 12,
13]. In the following subsection, we discuss the decomposition of the stress
tensor into its scalar, vector, and tensor (i.e., the transverse-traceless)
components444This is not to be conflated with the transverse-traceless gauge
in which the gravitational waves are computed far from the source., following
the prescription by Lifshitz [42], and in the next subsection, we will present
a brief overview of the different types of turbulent sources.

###  2.1 The scalar-vector-tensor decomposition of the stress tensor

According to the Lifshitz prescription, any symmetric rank-2 tensor
$\lambda_{ij}(\mathbf{x})$ in $d$-dimensions can be decomposed as [44]

$\lambda_{ij}=L\delta_{ij}+\nabla_{\langle
i}\nabla_{j\rangle}\lambda+\nabla_{(i}\bar{\lambda}_{j)}+\bar{\lambda}_{ij},$
(1)

where $\bar{\lambda}_{i}$ is divergenceless, $\nabla_{i}\bar{\lambda}_{i}=0$,
$\bar{\lambda}_{ij}$ is divergenceless and traceless,

$\nabla_{i}\bar{\lambda}_{ij}=\nabla_{j}\bar{\lambda}_{ij}=0,\quad\bar{\lambda}_{ii}=0,$
(2)

and

$\displaystyle\nabla_{\langle i}\nabla_{j\rangle}\lambda$
$\displaystyle\equiv\left(\nabla_{i}\nabla_{j}-\frac{1}{d}\delta_{ij}\nabla^{2}\right)\lambda,$
(3) $\displaystyle\nabla_{(i}\bar{\lambda}_{j)}$
$\displaystyle\equiv\frac{1}{2}\big{(}\nabla_{i}\bar{\lambda}_{j}+\nabla_{j}\bar{\lambda}_{i}\big{)},$
(4)

correspond to the removal of trace, and the symmetrization operators,
respectively. This decomposition is commonly referred to as a scalar-vector-
tensor (SVT) decomposition. Extracting the components in real space is
algebraically involved, and is discussed in detail in A.

The process, however, is less cumbersome in momentum space where the analog of
Eq. (1) can be written as,

$\tilde{\lambda}_{ij}=\tilde{L}\delta_{ij}-k_{\langle
i}k_{j\rangle}\tilde{\lambda}^{S}+ik_{(i}\tilde{\lambda}^{V}_{j)}+\tilde{\lambda}^{T}_{ij},$
(5)

where tildes denote quantities in Fourier space555We also adopt the notational
convention that a function $f(\mathbf{r})$ of real space coordinates
$\mathbf{r}$ is expressed in terms of its Fourier transform
$\tilde{f}(\mathbf{k})$ as $f(\mathbf{r})=(2\pi)^{-3}\int
d^{3}\mathbf{k}\,e^{-i\mathbf{k}\cdot\mathbf{r}}\tilde{f}(\mathbf{k})$, and
the inverse transform is defined accordingly., $k_{i}$ is the $i$-th component
of $\mathbf{k}$, and the superscripts $S$, $V$, and $T$ denote the scalar,
vector, and tensor contributions respectively. Analogous to the real space
expressions, we have
$k_{i}\tilde{\lambda}^{V}_{i}=0,k_{i}\tilde{\lambda}^{T}_{ij}=k_{j}\tilde{\lambda}^{T}_{ij}=0$,
and $\tilde{\lambda}^{T}_{ii}=0$, $k_{\langle i}k_{j\rangle}\equiv
k_{i}k_{j}-k^{2}\delta_{ij}/d$ with $k=|\mathbf{k}|$, and the parentheses in
the indices denote a symmetrization as in Eq. (3). $\tilde{L}$ is proportional
to the trace, i.e., $\tilde{L}=\tilde{\lambda}_{ii}/d$. The other scalar part
is extracted out as

$\tilde{\lambda}^{S}=-\frac{d}{d-1}\frac{\hat{k}_{i}\hat{k}_{j}}{k^{2}}\tilde{\lambda}_{\langle
ij\rangle}=-\frac{d}{d-1}\frac{\hat{k}_{i}\hat{k}_{j}}{k^{2}}\left(\tilde{\lambda}_{ij}-\frac{1}{d}\delta_{ij}\tilde{\lambda}_{ii}\right),$
(6)

where $\hat{k}_{i}\equiv k_{i}/k$, and the vector contribution turns out to be

$\displaystyle\tilde{\lambda}^{V}_{i}$
$\displaystyle=-\frac{2i}{k^{2}}\left(k_{j}\tilde{\lambda}_{\langle
ij\rangle}+\frac{d-1}{d}k_{i}k^{2}\tilde{\lambda}^{S}\right)$ (7)
$\displaystyle=-\frac{2i}{k^{2}}\left(k_{j}\tilde{\lambda}_{\langle
ij\rangle}-k_{i}\hat{k}_{p}\hat{k}_{q}\tilde{\lambda}_{\langle
pq\rangle}\right).$

Finally, the tensor part can thus be extracted from Eq. (5) as

$\displaystyle\tilde{\lambda}^{T}_{ij}$
$\displaystyle=\tilde{\lambda}_{\langle
ij\rangle}-\frac{d}{d-1}\hat{k}_{\langle
i}\hat{k}_{j\rangle}\hat{k}_{p}\hat{k}_{q}\tilde{\lambda}_{\langle pq\rangle}$
(8)
$\displaystyle-\left[\left\\{\hat{k}_{i}\hat{k}_{p}\tilde{\lambda}_{\langle
jp\rangle}+(i\longleftrightarrow
j)\right\\}-2\hat{k}_{i}\hat{k}_{j}\hat{k}_{r}\hat{k}_{s}\tilde{\lambda}_{\langle
rs\rangle}\right].$

The more conventional way to extract the transverse-traceless part of the
stress tensor is through the four-point projection operator
$P_{ijlm}(\hat{k})$, defined as
$P_{ijlm}(\hat{k})\equiv\left[P_{il}P_{jm}-\frac{1}{2}P_{ij}P_{lm}\right](\hat{k})$
with $P_{ij}(\hat{k})=\delta_{ij}-\hat{k}_{i}\hat{k}_{j}$, as is done, for
example, in Ref. [45] (see their Eq. (22)). It can be shown in a
straightforward way that the result of this operation gives the same
expression as in Eq. (8). We now apply this formalism to the energy-momentum
tensor for the magnetic and the plasma velocity fields.

###  2.2 SVT Decomposition of Magnetic Stress

For a magnetic field $\mathbf{B}$, the stress tensor can be written as

$\tau_{ij}(\mathbf{x})=B_{i}(\mathbf{x})B_{j}(\mathbf{x})-\frac{1}{2}\delta_{ij}B^{2}(\mathbf{x}),$
(9)

where $B_{i}$ are the spatial components of the three-dimensional vector
$\mathbf{B}$. The magnetic field $\mathbf{B}$ is solenoidal, i.e.,
$\nabla\cdot\mathbf{B}=0$. This tensor can be decomposed into its SVT
contributions using Eq. (1),

$\displaystyle\tau_{ij}$ $\displaystyle=T\delta_{ij}+\nabla_{\langle
i}\nabla_{j\rangle}\tau+\nabla_{(i}\bar{\tau}_{j)}+\bar{\tau}_{ij},$ (10)

where $\nabla_{i}\bar{\tau}_{i}=0$,
$\nabla_{i}\bar{\tau}_{ij}=\nabla_{j}\bar{\tau}_{ij}=0$, and
$\bar{\tau}_{ii}=0$. An example is given in B for a one-dimensional Beltrami
field. The term proportional to the trace is $T=-B_{i}B_{i}/6\equiv-B^{2}/6$,
and the other scalar component is (see A for the derivation of these
expressions for general $\lambda_{ij}$)

$\tau=-\frac{1}{2}\int dV^{\prime}\int
dV^{\prime\prime}\frac{\displaystyle{\nabla^{\prime\prime
2}\big{[}B^{2}(\mathbf{x}^{\prime\prime})\big{]}}-3\nabla^{\prime\prime}_{j}B_{i}(\mathbf{x}^{\prime\prime})\nabla^{\prime\prime}_{i}B_{j}(\mathbf{x}^{\prime\prime})}{16\pi^{2}|\mathbf{x}-\mathbf{x}^{\prime}||\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}|},$
(11)

and the vector contribution is

$\bar{\tau}_{j}=\int
dV^{\prime}\,\frac{\displaystyle{\nabla^{\prime}_{j}\big{\\{}B^{2}(\mathbf{x}^{\prime})}+2\nabla^{\prime
2}\tau(\mathbf{x}^{\prime})\big{\\}}-3B_{i}(\mathbf{x}^{\prime})\nabla^{\prime}_{i}B_{j}(\mathbf{x}^{\prime})}{6\pi|\mathbf{x}-\mathbf{x}^{\prime}|}.$
(12)

Note that the solenoidality of the magnetic field has been used to obtain Eqs.
(11–12). The tensor mode is then easily extracted by subtracting the scalar
and vector contributions,

$\bar{\tau}_{ij}=B_{\langle i}B_{j\rangle}-\nabla_{\langle
i}\nabla_{j\rangle}\tau-\nabla_{(i}\bar{\tau}_{j)},$ (13)

where, as before, $B_{\langle
i}B_{j\rangle}\equiv\left(B_{i}B_{j}-\delta_{ij}B^{2}/3\right)$. In Fourier
space, the stress tensor in Eq. (9) becomes a convolution,

$\displaystyle\tilde{\tau}_{ij}(\mathbf{k})$
$\displaystyle=\frac{1}{(2\pi)^{3}}\left(\tilde{B}_{i}\star\tilde{B}_{j}\right)(\mathbf{k})\equiv\frac{1}{(2\pi)^{3}}\int
d^{3}\mathbf{p}\,\tilde{B}_{i}(\mathbf{p})\tilde{B}_{j}(\mathbf{k}-\mathbf{p}),$
(14)

where $\tilde{B}_{i}(\mathbf{k})$ is the Fourier transform of
$B_{i}(\mathbf{x})$, and we have omitted the isotropic term in Eq. (9). For
the stress-tensor in Eq. (10), the expression analogous to Eq. (5) is,

$\tilde{\tau}_{ij}=\tilde{T}\delta_{ij}-k_{\langle
i}k_{j\rangle}\tilde{\tau}^{S}+ik_{(i}\tilde{\tau}^{V}_{j)}+\tilde{\tau}^{T}_{ij},$
(15)

with
$k_{i}\tilde{\tau}^{V}_{i}=0,k_{i}\tilde{\tau}^{T}_{ij}=k_{j}\tilde{\tau}^{T}_{ij}=0$,
and $\tilde{\tau}^{T}_{ii}=0$. The term proportional to the trace is
$\tilde{T}=(2\pi)^{-3}(\tilde{B}_{s}\star\tilde{B}_{s})(\mathbf{k})/3$, while

$\tilde{\tau}^{S}(\mathbf{k})=-\frac{1}{(2\pi)^{3}}\frac{3}{2}\frac{\hat{k}_{i}\hat{k}_{j}}{k^{2}}\tilde{B}_{\langle
i}\star\tilde{B}_{j\rangle}(\mathbf{k}),$ (16)

where, continuing the spirit of using angle bracket notations,
$\tilde{B}_{\langle
i}\star\tilde{B}_{j\rangle}(\mathbf{k})\equiv(\tilde{B}_{i}\star\tilde{B}_{j})(\mathbf{k})-\delta_{ij}(\tilde{B}_{s}\star\tilde{B}_{s})(\mathbf{k})/3$.
For clarity, we suppress the functional dependence on $\mathbf{k}$ henceforth.
The vector contribution turns out to be

$\tilde{\tau}^{V}_{i}=-\frac{2ik_{j}}{k^{2}}\left[\frac{1}{(2\pi)^{3}}\tilde{B}_{\langle
i}\star\tilde{B}_{j\rangle}+k_{\langle i}k_{j\rangle}\tilde{\tau}^{S}\right],$
(17)

and finally, the tensor part can thus be extracted as

$\tilde{\tau}^{T}_{ij}=\frac{1}{(2\pi)^{3}}\tilde{B}_{\langle
i}\star\tilde{B}_{j\rangle}+k_{\langle
i}k_{j\rangle}\tilde{\tau}^{S}-ik_{(i}\tilde{\tau}^{V}_{j)}.$ (18)

###  2.3 SVT Decomposition of Kinetic Stress

For a fluid with velocity $\mathbf{u}$, the stress tensor analogous to Eq. (9)
can be written as

$\sigma_{ij}(\mathbf{x})=w\gamma^{2}u_{i}u_{j}-p\delta_{ij},$ (19)

where $w=p+\rho$ is the enthalpy666In a magnetized plasma, an important
quantity associated with the magnetic field is the Alfvén velocity
$\mathbf{v}_{\rm A}$ defined as $\mathbf{v}_{\rm A}=\mathbf{B}/\sqrt{w}$., $p$
and $\rho$ are the pressure and density of the fluid respectively,
$\gamma=(1-u^{2})^{-1/2}$ is the Lorentz factor, and the spacetime dependence
of all the fields are suppressed777If we have a solenoidal flow, i.e.,
$\nabla\cdot\mathbf{u}=0$, then the anisotropic stress tensor for the fluid
velocity $\gamma u_{i}$ is analogous to that of magnetic fields.. The SVT
decomposition of this tensor, analogous to Eq. (10), can be written as

$\sigma_{ij}=S\delta_{ij}+\nabla_{\langle
i}\nabla_{j\rangle}\sigma+\nabla_{(i}\bar{\sigma}_{j)}+\bar{\sigma}_{ij},$
(20)

and as before, $\nabla_{i}\bar{\sigma}_{i}=0$,
$\nabla_{i}\bar{\sigma}_{ij}=\nabla_{j}\bar{\sigma}_{ij}=0$, and
$\bar{\sigma}_{ii}=0$. The mathematical details are more involved compared to
the case of magnetic fields, since the velocity fields are allowed to have a
longitudinal component, as opposed to the case with magnetic fields. In
addition, we consider non-relativistic fluid velocities ($\gamma\sim 1)$ for
simplicity of the expressions in the analytical results. However, this term is
recovered in the numerical treatment. We omit the general expressions for the
S, V, and T components of the stress tensor in real space, which are provided
in C.

The stress tensor can be written in $k$-space in a manner analogous to Eq.
(14) as

$\displaystyle\tilde{\sigma}_{ij}(\mathbf{k})$
$\displaystyle=\frac{1}{(2\pi)^{6}}\left((\tilde{p}+\tilde{\rho})\star\tilde{u}_{i}\star\tilde{u}_{j}\right)(\mathbf{k})-\tilde{p}(\mathbf{k})\delta_{ij}$
(21) $\displaystyle=\frac{1}{(2\pi)^{6}}\int d^{3}\mathbf{p}\int
d^{3}\mathbf{q}\,\bigg{[}(\tilde{p}+\tilde{\rho})(\mathbf{p})\tilde{u}_{i}(\mathbf{q})\tilde{u}_{j}(\mathbf{k}-\mathbf{p}-\mathbf{q})\bigg{]}$
$\displaystyle-\tilde{p}(\mathbf{k})\delta_{ij}.$

This is again decomposed as,

$\tilde{\sigma}_{ij}=\tilde{S}\delta_{ij}-k_{\langle
i}k_{j\rangle}\tilde{\sigma}^{S}+ik_{(i}\tilde{\sigma}^{V}_{j)}+\tilde{\sigma}^{T}_{ij},$
(22)

where $k_{i}\tilde{\sigma}^{V}_{i}=0,k_{i}\tilde{\sigma}^{T}_{ij}=0$, and
$\tilde{\sigma}^{T}_{ii}=0$. $\tilde{S}$ is again proportional to the trace,
$\tilde{S}=-\tilde{p}(\mathbf{k})+(2\pi)^{-6}((\tilde{p}+\tilde{\rho})\star\tilde{u}_{i}\star\tilde{u}_{i})(\mathbf{k})/3$,
and the other scalar part is

$\tilde{\sigma}^{S}(\mathbf{k})=-\frac{1}{(2\pi)^{6}}\frac{3}{2}\frac{\hat{k}_{i}\hat{k}_{j}}{k^{2}}(\tilde{p}+\tilde{\rho})\star\tilde{u}_{\langle
i}\star\tilde{u}_{j\rangle}(\mathbf{k}),$ (23)

and as in Eqs. (17) and (18), we have

$\tilde{\sigma}^{V}_{i}=-\frac{2ik_{j}}{k^{2}}\left[\frac{1}{(2\pi)^{6}}(\tilde{p}+\tilde{\rho})\star\tilde{u}_{\langle
i}\star\tilde{u}_{j\rangle}+k_{\langle
i}k_{j\rangle}\tilde{\sigma}^{S}\right],$ (24)

and,

$\tilde{\sigma}^{T}_{ij}=\frac{1}{(2\pi)^{3}}(\tilde{p}+\tilde{\rho})\star\tilde{u}_{\langle
i}\star\tilde{u}_{j\rangle}+k_{\langle
i}k_{j\rangle}\tilde{\sigma}^{S}-ik_{(i}\tilde{S}^{V}_{j)}.$ (25)

###  2.4 Modeling of Turbulent Sources

As mentioned above, the primordial plasma that generates the GWs is highly
turbulent, and contains magnetic fields. It is therefore essential to study
the MHD evolution of the plasma and understand the various kinds of turbulent
initial conditions. The purpose of this subsection is to provide an elementary
background of MHD turbulence, in preparation for the numerical results
discussed in Sec. 3.

Turbulent flows are characterized by the magnetic and fluid Reynolds numbers,
$\mathrm{Re}_{M}$ and $\mathrm{Re}$ respectively, defined as

$\mathrm{Re}_{M}\equiv\frac{vL}{\eta},\quad\quad\mathrm{Re}\equiv\frac{vL}{\nu},$
(26)

for a typical fluid velocity $v$ on the length scale $L$, and $\eta$ and $\nu$
being the magnetic diffusivity and kinematic viscosity respectively.
$\mathrm{Re}_{M}$ characterizes the relative importance of the diffusion term
over the induction term, while $\mathrm{Re}$ tells us about the importance of
the nonlinear velocity term relative to the viscous dissipation term.

In MHD simulations, one generally studies the evolution of the energy
spectra888These energy spectra are defined such that the mean energy density
$\mathcal{E}_{i}$ is given by $\mathcal{E}_{i}=\int dk\,E_{i}(k)$. $E_{i}(k)$
of the velocity and magnetic fields (for $i=K,M$ respectively). In a turbulent
medium, the energy spectra have a range of scales where the spectral shape is
given by the Kolmogorov spectrum, i.e., $E_{i}(k)\propto k^{-5/3}$ – this is
called the inertial subrange; see [46] for a review. Scales smaller than the
inertial range are dominated by viscous eddies, while larger scales have
energy carrying eddies. Turbulence involves nonlinear dynamics, which leads to
an energy cascade whereby energy is transferred from larger to smaller scales
(direct cascade) or vice versa (inverse cascade). The latter, for example, is
interesting in the context of evolution of primordial magnetic fields, since
it can explain the presence of magnetic fields correlated at sufficiently
large scales in the universe after they evolve through MHD evolution of the
plasma.

The physics governing the MHD evolution of the magnetic and velocity fields of
the plasma affects the instantaneous decay exponents of the mean energies. For
example, if we have helical magnetic fields in a highly conducting plasma,
magnetic helicity is conserved throughout the evolution; in that case,
magnetic energy decays with time $t$ as $\mathcal{E}_{M}(t)\propto t^{-2/3}$.
On the other hand, for non-helical fields, one has $\mathcal{E}_{M}(t)\propto
t^{-1}$. The correlation length of the turbulence, $\xi_{\rm M}$, which is
defined as a weighted average over the spectrum as $\xi_{\rm
M}=\mathcal{E}_{i}^{-1}\int dk\,k^{-1}E_{i}(k)$, tends to increase in a
corresponding fashion. The increase is fastest for helical turbulence, where
$\xi_{\rm M}\propto t^{2/3}$, and slower for nonhelical turbulence, where
$\xi_{\rm M}\propto t^{1/2}$ [47]. This can have a noticeable effect on the
resulting GW production [10].

##  3 Numerical scalar–vector–tensor decomposition of the stress tensor

###  3.1 Numerical modeling setup

In the following, we present the results of a numerical SVT decomposition of
the stress tensor. We consider four distinct types of isotropic hydrodynamic
and MHD turbulence similar to those considered in recent numerical studies of
GW production [48, 10]. These flows are distinguished by different types of
forcing; a summary is given in Table LABEL:Tsummary. The forcing energy is
injected at wave numbers around $k_{\ast}$, which we arrange to be six times
larger than the smallest wavenumber $k_{1}=2\pi/L$ in our triply-periodic
computational domain of side length $L$. For most of the cases, we choose
$k_{1}=100\,H_{*}$), and in two cases, we use $k_{1}=0.3\,H_{*}$). Energy is
dissipated through a Laplacian diffusion operator with coefficients $\nu$ and
$\eta$ for the viscosity and the magnetic diffusivity, which act on the
smallest resolved length scales. We use conformal time normalized to the
Hubble time at the time of source initiation $H_{\star}^{-1}$, which
corresponds to $t=1$ in our units. Following the recent work of Ref. [10], we
apply a forcing term (in the momentum equation or in the induction equation
for kinetic or magnetic turbulence, respectively) for one Hubble time until
$t=2$. The forcing has an initial stage $1<t<2$ when the forcing term is
present, an intermediate stage $2<t<3$ when it decreases linearly, and a final
stage $t>3$ when it is absent. The forcing is random and applied at each time
step. This allows to sustain stationary turbulence for a finite time during
the intermediate stage, which leads to freely decaying turbulence in the final
stage. For Run A, we apply forcing in the form of spherical expansion waves
through the gradients of randomly placed Gaussians with amplitude $f_{0}$
[49], which are irrotational, i.e., curl-free. This case was also studied in
Ref. [48], and is similar to their Run ac1, except that in their case the
forcing was only applied until $t=1.1$, and the intermediate stage is not
considered. In Run B, we also apply irrotational forcing, but in the form of
plane waves, again with amplitude $f_{0}$. The remaining four runs also have
plane wave forcings, but are all vortical, calculated as a curl of a randomly
oriented vector with amplitude $f_{0}$ and can have helicity, which is
quantified by the parameter $\sigma\in[0,1]$. Run C has nonhelical forcing,
while Run D has a helical forcing. Runs E and F have magnetic fields and the
forcing is applied in the induction equation, again first without helicity
(Run E) and with helicity (Run F). In Runs G and H, forcing is applied at
small $k$ ($k_{*}=2$) during a time interval $1\leq t\leq 2$, using either
vortical (Run G) or irrotational (Run H) forcing. In all cases, we present the
results of the SVT decomposition in the form of power spectra per linear wave
number interval for the projected part of the stress tensor. The efficiency of
GW production will be related to the contribution by the tensor mode with
respect to the scalar and vector modes of the turbulent source. The energy
injection wave numbers are $k=k_{*}\equiv 600$ (Runs A–F) and 2 (Runs G and
H). The resulting spectra of the stress peak at $k=2k_{*}\approx 1200$ and 4,
respectively. Note that we have used normalized wave numbers $k$ such that
$k=1$ corresponds to the Hubble scale $H_{*}^{-1}$ at the time of generation.

###  3.2 Comparison for Acoustic and Vortical Turbulence

In Table LABEL:Tsummary, we list the maximum kinetic or magnetic energies,
${\cal E}_{i}^{\max}$ for $i={\rm K}$ and ${\rm M}$, respectively. These are
normalized by the radiation energy, just like the resulting GW energy, ${\cal
E}_{\rm GW}^{\rm sat}$, which refers to the GW energy averaged over
oscillations around the saturated stationary value. In earlier work of Ref.
[48], the saturation GW energy ${\cal E}_{\rm GW}^{\rm sat}$ was found to be
proportional to the square of the ratio ${\cal E}_{i}^{\rm max}/k_{*}$. To
characterize the efficiency in GW energy conversion, we quote in Table
LABEL:Tsummary the factor $q=k_{*}\,({\cal E}_{\rm GW}^{\rm sat})^{1/2}/{\cal
E}_{i}^{\max}$. In Ref. [48], the values of $q$ were found to be between one
and ten, but here we find even larger values.

We begin with the case of acoustic turbulence in Runs A and B. This type of
turbulence is unusual in that it processes very little vorticity. In Fig. 1,
we show the power spectra of the projected stress tensor for acoustic
turbulence. Those flows have always zero helicity. The spectrum has a
subinertial range proportional to $k^{3/2}$ and an inertial range which is
proportional to $k^{-2}$ for Run A and proportional to $k^{-1.8}$ for Run B,
similar to what is expected for acoustic turbulence [50]. The scalar and
vector components dominate over the tensor components at almost all wave
numbers, although the scalar component becomes smaller than the tensor
component at small wave numbers ($k<400$). In addition, at small wave numbers,
the vector components dominate over the scalar component.

Table 1:  Summary of the runs discussed in the paper.

Run $i$ $k_{\ast}$ ${\cal E}_{i}^{\max}$ $\sigma$ wave ${\cal E}_{\rm GW}^{\rm
sat}$ $\nu$ ($=\eta$) $f_{0}$ $q$ A K 600 $3.8\times 10^{-3}$ 0 expan
$1.5\times 10^{-8}$ $10^{-6}$ $1$ 19 B K 600 $2.4\times 10^{-3}$ 0 plane
$1.6\times 10^{-8}$ $2\times 10^{-7}$ $2\times 10^{-3}$ 32 C K 600 $8.0\times
10^{-3}$ 0 plane $4.5\times 10^{-9}$ $2\times 10^{-7}$ $4\times 10^{-1}$ 4.3 D
K 600 $1.1\times 10^{-2}$ 1 plane $6.3\times 10^{-9}$ $2\times 10^{-7}$
$4\times 10^{-1}$ 5.0 E M 600 $5.7\times 10^{-3}$ 0 plane $1.5\times 10^{-9}$
$2\times 10^{-7}$ $6\times 10^{-4}$ 4.1 F M 600 $1.7\times 10^{-2}$ 1 plane
$5.1\times 10^{-9}$ $2\times 10^{-7}$ $6\times 10^{-4}$ 2.5 G K 2 $4.2\times
10^{-2}$ 0 plane $8.3\times 10^{-4}$ $2\times 10^{-2}$ $7\times 10^{-1}$ 1.4 H
K 2 $4.7\times 10^{-2}$ 0 plane $1.1\times 10^{-3}$ $1\times 10^{-2}$ $4\times
10^{-1}$ 1.4

Figure 1:  SVT decomposition for Runs A (left) and B (right) for irrotational
turbulence produced by spherical expansion waves and plane wave forcing,
respectively.

The situation is quite different for vortical turbulence. We consider here
cases without kinetic helicity (Run C) and with kinetic helicity (Run D).
These two runs are similar in that they have a much shallower inertial range
spectrum proportional to $k^{-1}$. In addition, we also consider two cases
with magnetic fields, which are driven by an electromotive force. The
functional form is the same as for the forcing in the momentum equation, but
the amplitude factor $f_{0}$ tends to be much smaller while still leading to
comparable stresses. This is a consequence of this forcing occurring
underneath the curl of the induction equation. Run E has null magnetic
helicity and Run F has maximum magnetic helicity; see Fig. 3. In these two
cases, we find a somewhat steeper inertial range spectrum proportional to
$k^{-5/3}$, but for Run F there is a clear uprise of the spectrum at low
wavenumber, which is likely a consequence of an inverse cascade; see also [40]
for hydromagnetic decay simulations showing inverse cascading. Those cases are
rather similar in that now the tensor component is dominant; see Fig. 3.

Figure 2:  SVT decomposition for Runs C (left) and D (right) for vortical
hydrodynamic turbulence without and with helicity, respectively.

Figure 3:  SVT decomposition for Runs E (left) and F (right) with magnetic
forcing without and with helicity, respectively.

As in the case of acoustic turbulence, the vector mode dominates over the
scalar mode at the largest length scales ($k<400$), with Run F being the
exception (where the vector and scalar modes are comparable, and both much
smaller than the tensor mode). For vortical hydrodynamic turbulence, the
tensor and vector modes are comparable in the inertial range, while for
magnetic turbulence the tensor mode dominates over this regime. Hence, in the
magnetic turbulence case, the spectrum of the total unprojected stress tensor
agrees surprisingly well with that of the tensor mode (for all $k$ in the
helical case, and at large $k$ for the non-helical case), so the computational
effort for doing such a projection could potentially be ignored in a first
approximation, as was done in Refs. [51, 52]. This is, however, not the case
for acoustic turbulence, where the projection is essential.

It is worth noting that all spectra of Figs. 1–3 show a subinertial range
proportional to $k$ or $k^{3/2}$, which is shallower than the generic $k^{2}$
spectrum. In the present case, the subinertial range is not long enough to
make strong claims, but it is useful to recall that deviations from a $k^{2}$
spectrum of the stress have been found when there are departures from
Gaussianity of the underlying magnetic or velocity fields [53]. Although the
departures from Gaussianity may not yet be very strong in the present
simulations, it is worth noting that real astrophysical magnetic fields have
much larger magnetic Reynolds numbers and may well lead to much more
significant departures.

###  3.3 On the Relative GW Strength for Acoustic and Vortical Turbulence

In this connection, it is interesting to note that the apparent strength of
GWs for acoustic turbulence cannot straightforwardly be explained by the SVT
decomposition. If there was such a connection, one would have expected the GW
energy to be less for acoustic turbulence than for vortical turbulence,
because the tensor mode (i.e., the actual source for the GWs) is weaker there.
However, this is not the case [48]. We must therefore conclude that this
difference between acoustic and vortical turbulence is explained by other
factors, most likely the detailed aspects of the temporal variation of the
stress. In fact, to drive GWs, the temporal stress frequency must be close to
the frequency of light at each spatial Fourier mode and it seems plausible
that this is more easily satisfied for acoustic turbulence. It can be argued
that only the parallel component of the perturbations in magnetic fields
contributes to the energy spectrum of the stress tensor of $B^{2}$ sourcing
the GWs. Thus the acoustic modes of turbulence are more efficient in
transferring energy to the GWs than the vortical modes, and thus GWs sourced
by acoustic turbulence carry more energy compared to those by vortical
turbulence.

In Fig. 4 we show the positions of Runs A–F in a ${\cal E}_{\rm GW}^{\rm sat}$
versus ${\cal E}_{\rm M}^{\max}$ diagram. For orientation, we also show the
old data points from the Ref. [48]. We see that the new acoustic or
irrotational turbulence simulations are now even slightly more efficient in
driving GWs: Run B has not only the smallest kinetic energy of all the runs,
but it also has a relatively large GW energy.

Larger GW energies can be obtained by decreasing $k_{*}$. In Fig. 5, we show
the results for $k_{*}=2$ and compare vortical and irrotational hydrodynamic
forcings. Similarly to the cases with larger $k_{*}$, we see, again, that the
tensor mode of acoustic turbulence is weaker (second panel of Fig. 5), but
this is only true of the initial range, which contributes little to the total
GW energy production. Most of the GW energy comes from wave numbers $k$ below
the peak, and there the tensor mode is, similarly to the cases with larger
$k_{*}$, the largest (left panel of Fig. 5 around the peak, and right panel
for low $k$) or, at least comparable to the vector mode, which is the largest
(left panel of Fig. 5 for low $k$).

Figure 4:  Positions of Runs A–F in a diagram showing ${\cal E}_{\rm GW}^{\rm
sat}$ versus ${\cal E}_{\rm M}^{\max}$, For orientation the data points of the
Ref. [48] are shown as gray symbols for acoustic (ac), helical (hel),
initialized (ini) turbulence.

Figure 5:  SVT decomposition for Runs E (left) and F (right) with $k_{\ast}=2$
using irrotational and vortical forcings, respectively.

##  4 Relation Between Stress and Strain Spectra

###  4.1 GW Equation for the Radiation Dominated Epoch

In this section, we connect the observable quantities related to GWs to
properties of the source obtained in the previous section. GWs are considered
to be small tensor perturbations over the Friedmann-Lemaître-Robertson-Walker
(FLRW) metric, and the dynamics of these GWs is given by the linearized
Einstein equations. The perturbed metric is taken to be,

$ds^{2}=-dt_{\rm P}^{2}+a^{2}(t_{\rm P})[\delta_{ij}+h_{ij,{\rm
P}}(\mathbf{x},t)]\,dx_{i}dx_{j},$ (27)

where $\mathbf{x}$ are the comoving coordinates, $t_{\rm P}$ is the physical
time, $a(t_{\rm P})$ is the scale factor of the universe, and $h_{ij,{\rm P}}$
are the physical strains. The equation governing the evolution of these GWs is
given by [43, 44],

$\left[\frac{\partial^{2}}{\partial t_{\rm P}^{2}}+3H(t_{\rm
P})\frac{\partial}{\partial t_{\rm
P}}-\frac{\nabla^{2}}{a^{2}}\right]h_{ij,{\rm P}}(\mathbf{x},t_{\rm P})=16\pi
GT^{\rm TT}_{ij,{\rm P}}(\mathbf{x},t_{\rm P}),$ (28)

where $G$ is Newton’s constant, $H(t_{\rm P})=da/(a\,dt_{\rm P})$ is the
Hubble rate, $T^{\rm TT}_{ij,{\rm P}}$ is the physical stress tensor, the
superscript TT denotes the the transverse and traceless (TT) projection, and
$\nabla$ signifies derivatives with respect to the comoving coordinates
$\mathbf{x}$.

The energy density in the magnetic fields and plasma velocity dilutes with the
expansion of the universe, and to scale out the effects of this expansion, we
need to express the GW equation in terms of comoving stress tensor
$T_{ij}=a^{4}\,T_{ij,{\rm P}}$, the scaled strain $h_{ij}=a\,h_{ij,{\rm P}}$,
and the conformal time $t$ defined as $dt_{\rm P}=a\,dt$. Also in the
radiation dominated epoch, the scale factor is linear in conformal time, and
the damping term in the above equation vanishes. Eq. (28) can then be
expressed as

$h^{\prime\prime}_{ij}-\nabla^{2}h_{ij}=\frac{16\pi G}{a}T^{\rm
TT}_{ij}(\mathbf{x},t),$ (29)

where the primes denote derivatives with respect to the conformal time $t$.
The analytic solution of the GW equation is discussed in several references;
see, for example, Refs. [51, 32]. The full numerical solution for a turbulent
source is also discussed in Refs. [54, 48].

The four-dimensional power spectrum $H_{ijlm}(\mathbf{k},t,\tau)$ of the
energy density tensor is defined as

$\frac{1}{w^{2}}\left<\tilde{T}^{*}_{ij}(\mathbf{k},t)\tilde{T}_{lm}(\mathbf{k}^{\prime},t+\tau)\right>=(2\pi)^{3}\delta^{3}(\mathbf{k}-\mathbf{k}^{\prime})H_{ijlm}(\mathbf{k},t,\tau),$
(30)

where $\tilde{T}_{ij}(\mathbf{k},t)$ is the Fourier transform of
$T_{ij}(\mathbf{x},t)$, we have used the reality condition
$\tilde{T}_{ij}(-\mathbf{k},t)=\tilde{T}^{*}_{ij}(\mathbf{k},t)$. As mentioned
previously, $w\equiv\rho+p$ is the enthalpy density of the primordial plasma,
where $\rho$ and $p$ being the energy density and the pressure of the
primordial plasma, and all quantities are comoving.

In the standard framework there are two degrees of freedom for $h_{ij}$ (a
$3\times 3$ tensor has 9 components, the symmetry condition $h_{ij}=h_{ji}$
leaves 6 degrees of freedom, the transverse nature of GWs, i.e.,
$k_{i}h^{ij}=0$, reduces that to 3, while the traceless condition $h_{ii}=0$
further reduces the total number of the degrees of freedom to 2) that
correspond to two possible states of GWs polarization. In the case of a parity
even source, these two states are equal to each other, while a parity odd
sources generates circularly polarized GWs, where there is an excess of one
state over the other. The circular polarization of GWs can be understood
intuitively by working in the circular polarization basis, which can be
constructed from the unit vector $\hat{{\bf k}}$, and two more unit vectors
$\hat{{\bf e}}^{1}$ and $\hat{{\bf e}}^{2}$ which are perpendicular to
$\hat{{\bf k}}$ and to each other. The basis is defined as [55]

$\hat{e}^{R,L}_{ij}(\hat{\mathbf{k}})=-\frac{1}{2}\big{(}\hat{e}_{1}\pm
i\hat{e}_{2}\big{)}_{i}\big{(}\hat{e}_{1}\pm i\hat{e}_{2}\big{)}_{j}.$ (31)

It is easy to check that $(e_{ij}^{L})^{*}=e_{ij}^{R}$, and vice versa; and
also that

$\displaystyle\hat{e}^{L}_{ij}(\hat{\mathbf{k}})\hat{e}^{R}_{ij}(\hat{\mathbf{k}})$
$\displaystyle=1,$ (32)
$\displaystyle\hat{e}^{L}_{ij}(\hat{\mathbf{k}})\hat{e}^{L}_{ij}(\hat{\mathbf{k}})$
$\displaystyle=\hat{e}^{R}_{ij}(\hat{\mathbf{k}})\hat{e}^{R}_{ij}(\hat{\mathbf{k}})=\hat{\mathbf{k}}_{i}\hat{e}^{R,L}_{ij}=\delta_{ij}\hat{e}^{R,L}_{ij}=0.$
(33)

The Fourier space stress tensor can now be decomposed into the left and right
handed components,

$\tilde{T}_{ij}(\mathbf{k},t)=\tilde{T}_{L}(\mathbf{k},t)\,e^{L}_{ij}(\hat{{\bf
k}})+\tilde{T}_{R}(\mathbf{k},t)\,e^{R}_{ij}(\hat{{\bf k}}).$ (34)

We can also construct symmetric and antisymmetric projectors [55],

$\displaystyle\mathcal{S}_{ijlm}(\hat{\mathbf{k}})$
$\displaystyle\equiv\frac{1}{2}\left[e^{L}_{ij}e^{R}_{lm}+e^{R}_{ij}e^{L}_{lm}\right],$
(35) $\displaystyle\mathcal{A}_{ijlm}(\hat{\mathbf{k}})$
$\displaystyle\equiv\frac{1}{2i}\left[e^{R}_{ij}e^{L}_{lm}-e^{L}_{ij}e^{R}_{lm}\right].$
(36)

These projection operators satisfy the properties that
$\mathcal{S}_{ijlm}(\hat{\mathbf{k}})\mathcal{S}_{ijlm}(\hat{\mathbf{k}})=\mathcal{A}_{ijlm}(\hat{\mathbf{k}})\mathcal{A}_{ijlm}(\hat{\mathbf{k}})=2$,
and
$\mathcal{S}_{ijlm}(\hat{\mathbf{k}})\mathcal{A}_{ijlm}(\hat{\mathbf{k}})=0$,
and pick out respectively the symmetric and antisymmetric components of a
fourth-order tensor.

The degree of polarization of the source can then be expressed as

$\mathcal{P}_{T}(k,t)=\frac{T^{\mathcal{A}}(k,t)}{T^{\mathcal{S}}(k,t)},$ (37)

where

$T^{\alpha}(k,t)=\frac{1}{2w^{2}}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int
d^{3}\mathbf{k}^{\prime}\alpha_{ijlm}(\hat{\mathbf{k}})\left<\tilde{T}^{*}_{ij}(\mathbf{k},t)\tilde{T}_{lm}(\mathbf{k}^{\prime},t)\right>,$
(38)

with $\alpha$ denoting either $\mathcal{S}$ or $\mathcal{A}$.

In terms of the polarization components of Eq. (34), we can write

$\displaystyle T^{\mathcal{S}}(k)$
$\displaystyle=\frac{1}{2w^{2}}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int
d^{3}\mathbf{k}^{\prime}\left<\tilde{T}^{*}_{R}(\mathbf{k})\tilde{T}_{R}(\mathbf{k}^{\prime})+\tilde{T}^{*}_{L}(\mathbf{k})\tilde{T}_{L}(\mathbf{k}^{\prime})\right>,$
(39) $\displaystyle T^{\mathcal{A}}(k)$
$\displaystyle=\frac{1}{2w^{2}}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int
d^{3}\mathbf{k}^{\prime}\left<\tilde{T}^{*}_{R}(\mathbf{k})\tilde{T}_{R}(\mathbf{k}^{\prime})-\tilde{T}^{*}_{L}(\mathbf{k})\tilde{T}_{L}(\mathbf{k}^{\prime})\right>,$
(40)

where the explicit $t$-dependence of the quantities in Eq. (39) is dropped for
brevity. $\mathcal{P}_{T}(k,t)$ thus captures the relative inequality of power
in the left and right handed components of the stress tensor. One can
alternatively work in the linear polarization basis, where the physical
meaning of the polarization degree of the source is not so straightforward.
The definition of the basis, as well as the expressions analogous to those
above are presented in D. For more details about the linear and circular
polarization bases, and the decomposition of the four-dimensional power
spectrum into its symmetric and antisymmetric parts, see Ref. [55], where the
stress tensor is sourced by magnetic fields. Notably, the analytical estimates
for the polarization degree differ significantly from the numerical
simulations due to the simplified treatment of the time-decorrelation
functions and the neglect of different scaling laws for the energy density and
helicity decays scalings; see Ref. [10] for more discussion.

Since GWs are sourced by the TT projected stress tensor, i.e., the tensor
component $\bar{\lambda}_{ij}$ of the stress tensor [cf. Eq. (1)], we built
the four-dimensional power spectrum $H_{ijlm}(\mathbf{k},t,\tau)$ from this TT
component and decomposed it into parity odd and even parts. It is of interest
to see whether it is possible to ascribe a degree of polarization analogous to
that in Eq. (37) to the scalar and vector components of the stress tensor. We
defer a discussion of that to E

###  4.2 Shell-integrated Spectra

In a manner analogous to Eq. (30), we can form a four-dimensional power
spectrum $Q_{ijlm}(\mathbf{k},t,\tau)$ from the time derivative of the strains
in Fourier space

$\frac{1}{32\pi G}\left<\dot{\tilde{h}}^{*}_{ij,{\rm
P}}(\mathbf{k},t)\dot{\tilde{h}}_{lm,{\rm
P}}(\mathbf{k}^{\prime},t+\tau)\right>=(2\pi)^{3}\delta^{3}(\mathbf{k}-\mathbf{k}^{\prime})Q_{ijlm}(\mathbf{k},t,\tau).$
(41)

For more details on the fourth order correlation functions and their Fourier
transforms, see Ref. [56].

The polarization degree of GWs is written as

$\mathcal{P}_{GW}(k,t)=\frac{H_{GW}(k,t)}{E_{GW}(k,t)},$ (42)

where $E_{GW}(k,t)$ and $H_{GW}(k,t)$ are respectively the spectral GW energy
density and the spectral GW helical energy density, obtained respectively by
the application of the projection operators of Eq. (35) on the LHS of Eq.
(41), in a way similar to Eq. (38).

In the circular polarization basis, we get (dropping the explicit
$t$-dependence)

$\displaystyle E_{GW}(k)$ $\displaystyle=\frac{1}{(2\pi)^{6}}\frac{1}{32\pi
G}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int d^{3}\mathbf{k}^{\prime}$ (43)
$\displaystyle\times\left<\dot{\tilde{h}}^{*}_{R}(\mathbf{k})\dot{\tilde{h}}_{R}(\mathbf{k}^{\prime})+\dot{\tilde{h}}^{*}_{L}(\mathbf{k})\dot{\tilde{h}}_{L}(\mathbf{k}^{\prime})\right>,$
$\displaystyle H_{GW}(k)$ $\displaystyle=\frac{1}{(2\pi)^{6}}\frac{1}{32\pi
G}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int d^{3}\mathbf{k}^{\prime}$ (44)
$\displaystyle\times\left<\dot{\tilde{h}}^{*}_{R}(\mathbf{k})\dot{\tilde{h}}_{R}(\mathbf{k}^{\prime})-\dot{\tilde{h}}^{*}_{L}(\mathbf{k})\dot{\tilde{h}}_{L}(\mathbf{k}^{\prime})\right>,$

where $\dot{\tilde{h}}_{ij,{\rm P}}(\mathbf{k},t)$ is decomposed in the
circular polarization basis,

$\dot{\tilde{h}}_{ij}(\mathbf{k},t)=\dot{\tilde{h}}_{L}(\mathbf{k},t)\,e^{L}_{ij}(\hat{{\bf
k}})+\dot{\tilde{h}}_{R}(\mathbf{k},t)\,e^{R}_{ij}(\hat{{\bf k}}).$ (45)

As expected, $\mathcal{P}_{GW}(k,t)$ captures the relative inequality of power
in the left and right handed components of the GW perturbations.

Figure 6:  Spectra of stress (solid line at maximum) and strain (dotted lines
for different times after the maximum) for Runs B, C, E, and F.

###  4.3 Numerical Results for the Strain and Stress Spectra

We now solve for the resulting GWs numerically using the approach described in
Refs. [54, 48]. Here we focus on a comparison of the results with irrotational
plane wave (Run B) and vortical plane wave forcings without (Run C) and with
(Run D) magnetic field, and with helicity (Run F). In each of the cases, we
show the spectra of both the TT-projected stress of the magnetic (kinetic)
fields $\bar{\tau}_{ij}$ ($\bar{\sigma}_{ij}$), and the strain $h$ scaled with
a $k^{2}$ factor. Here we denote the latter simply as ${\rm{Sp}}(k^{2}h)$, but
mean by this all of the components of the strain tensor. Likewise, we denote
the stress spectrum by ${\rm{Sp}}(\tau)$, and mean by this the TT-projected
stress tensor with the $16\pi G$ factor included, so that it would be equal to
${\rm{Sp}}(k^{2}h)$, if the second conformal time derivative of the strain in
Eq. (29) was absent. The result is shown in Fig. 6 for Runs B, C, E, and F.

We see from Fig. 6(d) that in the magnetic case with helicity, the spectra of
${\rm{Sp}}(k^{2}h)$ and ${\rm{Sp}}(\tau)$ are in almost perfect agreement with
each other at all wave numbers, except near $2k_{\ast}$. The situation changes
when there is no helicity and ${\rm{Sp}}(k^{2}h)$ is now slightly above
${\rm{Sp}}(\tau)$; see Fig. 6(c). The difference increases further when we
compare with purely hydrodynamic cases (Run C); see Fig. 6(b). Particularly
striking is a difference for irrotational turbulence, where we see significant
differences between the two spectra at all wave numbers (Run B); see Fig.
6(a).

Figure 7:  Spectra of stress (solid line at maximum) and strain (dotted lines
for different times after the maximum) for Runs G and H.

Our results suggest that, although there is a correspondence between
${\rm{Sp}}(k^{2}h)$ and ${\rm{Sp}}(\tau)$, there can also be a significant
departure and even a shift between them. Thus, even though the TT-projection
of the stress has resulted in a significant decrease in the projected stress
(see Fig. 1), the stress in the irrotational hydrodynamic cases have now other
qualities that make it more efficient for driving GWs. As already discussed
above, this could be caused by larger power at larger temporal frequencies.
This is supported by our findings that in the cases with the smallest
efficiency (Run D), the two spectra are almost equal, i.e.,
${\rm{Sp}}(\tau)\approx{\rm{Sp}}(k^{2}h)$.

###  4.4 Polarization Spectra of Stress and Strain

It is common to define the degree of polarization in terms of the strains
themselves, as opposed to our definition in terms of the time derivatives of
the strains (see, e.g., Eq. (6) of Ref. [32]). It can be argued that these
definitions are equivalent for GWs sourced by stationary turbulence, where one
can apply the aeroacoustic approximation. For stationary turbulence, one can
take the time derivatives out of the time averaging operation on the four-
point correlation function of the stress tensor (see Eq. (12) of Ref. [51]).
The same reasoning applies to the four-point correlation function of the GW
strain, since the strains are expected to have the same statistical properties
as the turbulent source. This means

$\displaystyle\left<\dot{h}^{\rm TT}_{ij,{\rm P}}(\mathbf{x},t)\dot{h}^{\rm
TT}_{lm,{\rm P}}(\mathbf{x}+\bm{\xi},t+\tau)\right>_{t}$
$\displaystyle=-\partial_{t}^{2}\left<h^{\rm TT}_{ij,{\rm
P}}(\mathbf{x},t)h^{\rm TT}_{lm,{\rm
P}}(\mathbf{x}+\bm{\xi},t+\tau)\right>_{t},$ (46)

where the subscript $t$ is used to denote that the angle brackets now also
include time averaging, in addition to the ensemble averaging. In frequency
space, the time derivatives become factors of frequency, and cancel out when
calculating the degree of polarization (see Eq. (42)).

Figure 8:  Degree of polarization of $T^{\rm TT}$ for (a) Runs K0, K01, K03,
K05, and K1 of Ref. [10], and (b) their Runs M01, M03, M05, and M1, as well as
the degree of polarization of $h$ (c) for Runs K0, K01, K03, K05, and K1, and
(d) for Runs M01, M03, M05, and M1. The inserts show the same, but in a
double-logarithmic representation. For (b), the straight line has a slope of
$-0.8$, while for (c) and (d) the two straight lines have a slope of $-0.55$.

Let us now discuss the degree of polarization of the source, ${\cal
P}_{T}(k)$, and of the strain. ${\cal P}_{h}(k)$, from numerical simulations.
In Fig. 8(a) and (b) we show ${\cal P}_{T}(k)$ for the source for the runs of
Ref. [10] for kinetic and magnetic driving. The run names for the kinetically
driven cases are K0, K01, K03, K05, and K1 for $\sigma=0$, $0.1$, $0.3$,
$0.5$, and $1$, respectively. The magnetically driven runs have analogous
names M0, M01, M03, M05, and M1.

Remarkably, in the kinetically driven case, the polarization of the source
never reaches 100%, and yet, the polarization of $h$ does reach nearly 100%.
This was already found in Ref. [10], and their results are convenient
comparison also reproduced in Fig. 8(c) and (d). Furthermore, in the
magnetically driven case, the source reaches nearly 100%, but not at the same
wave number where the GWs become 100% polarized. Instead, maximum polarization
happens here at the smallest scale.

##  5 Conclusions

We have computed the SVT decomposition for six different flows: kinetically
and magnetically driven flows, with and without helicity, as well as two types
of acoustically driven flows. In addition, we have considered two cases with
significantly smaller forcing wave numbers. All our flows had a turbulent
inertial range spectrum, as well as a subinertial range spectrum. However,
significant differences can be seen in the relative strengths of scalar,
vector, and tensor modes of the turbulent source. The largest contribution to
the tensor mode comes from the flows in Runs E and F. Both are magnetically
driven, although only one of them had helicity. In the presence of magnetic
helicity (Run F), an inverse cascade emerges, which leads to an additional
contribution to the tensor mode at small wave numbers. The scalar and vector
modes of the turbulence stress are small and of comparable strengths in the
helically driven case. Without helicity, however, the vector mode exceeds the
scalar mode at small and intermediate length scales. For helically driven
flows, there is no significant inverse cascade, and the tensor mode is no
longer so dominant. Instead, tensor and vector modes are now of comparable
strength. At smaller wave numbers, the vector mode does even become dominant.

For acoustically or irrotationally driven flows, the situation is drastically
different in that the tensor mode is now comparatively weak, while scalar and
vector modes are dominant and of similar strengths. Furthermore, there is very
little difference between two fundamentally different realizations of acoustic
turbulence: plane wave and spherical expansion wave forcings. Both show a very
strong “bottleneck” effect, i.e., a shallower spectrum at high wave numbers,
which was already found in Ref. [49]. Bottleneck effects are known to exist
also for vortical turbulence [57], but they are not as strong [59]. The planar
wave forcing was thought to be more similar to the vortical plane wave
forcing, but this is not the case. Even though the acoustically driven flows
have comparatively weak tensor modes, they are most effective in driving GWs.
The reason for this higher efficiency of acoustic flows in generating GWs is
not very well understood; but it might be interesting to note that at least in
a gravitationally stratified fluid, maximum acoustic wave power output was
found for flows with significant longitudinal modes compared to transverse
ones [58]. The planar wave-driven acoustic flows are even slightly more
effective (Run B, $q=32$), than spherical expansion wave-driven flows (Run A,
$q=19$).

Given that the presence of helicity had no significant effect on the results
of the SVT decomposition, it was of interest to look for differences in the
degree of polarization of the source and the resulting GWs. The difference is
indeed immense, but it results almost entirely from the inverse cascade, which
is only possible in the magnetically driven case. In these cases, the degree
of polarization of the source is virtually independent of the fractional
helicity. This is probably because, even for weak magnetic helicity, inverse
cascading is possible. At the driving scale, however, the degree of
polarization is never close to 100% – neither in the magnetically driven
cases, nor in the kinetically driven ones. In the latter, however, fractional
helicity of the velocity does lead to a smaller fractional helicity of the
stress. By contrast, in the magnetically driven cases, all spectra of the
polarization degree are similar and they seem to be controlled by some kind of
forward cascade from the largest scale downward to all smaller scales. This
includes even the driving scale of the system, whose degree of polarization
seems to be controlled by the stress on the largest scale of the domain.

Clearly, spectral power on the scale of the domain must be interpreted with
care, because the intention of periodic box simulations is to represent an
infinite and unbounded domain. The results will therefore be different if one
considered a larger domain. It is conceivable, then, that inverse cascading
also works in an infinitely extend domain, but, again, only on the scale of
the domain will the polarization degree be close to maximum.

The work presented here will be important for identifying the nature of the
turbulent motions in the early universe once a GW signal from turbulent
sources will be detected. More specifically, it can help to disentangle the
turbulent conditions in the early universe regarding the dependence on
irrotational (longitudinal) and vortical (solenoidal) flows.

Data availability—The source code used for the simulations of this study, the
Pencil Code, is freely available [60]. The simulation setups and the
corresponding data are freely available from Refs. [61, 62].

## Acknowledgments

Support through the Swedish Research Council, grant 2019-04234, and Shota
Rustaveli GNSF (grants 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.

## ORCID iDs

Axel Brandenburg: https://orcid.org/0000-0002-7304-021X
Grigol Gogoberidze: https://orcid.org/0000-0002-0506-333X
Tina Kahniashvili: https://orcid.org/0000-0003-0217-9852
Sayan Mandal: https://orcid.org/0000-0003-3129-5799
Alberto Roper Pol: https://orcid.org/0000-0003-4979-4430

##  Appendix A Helmholtz Theorem and the real space SVT decomposition

The scalar-vector-tensor (SVT) decomposition is the most generalized way to
decompose linearized perturbations of a symmetric rank-2 tensor based on how
these components transform under rotations. In this formulation, the evolution
equations for the scalar, vector, and tensor components of the metric
perturbations are decoupled from each other at the linear level. It is a more
generalized version of the Helmholtz decomposition by which a vector is
decomposed into divergence-free and curl-free parts.

In Sec. 2.1, we have provided the expressions for the SVT decomposition of a
general symmetric rank-2 tensor in momentum space. In this Appendix, we show
the corresponding decomposition in real space.

The Helmholtz theorem states that any 3-vector $b_{i}$ can be decomposed into
a scalar part $b$ and a vector part $\bar{b}_{i}$,

$b_{i}=\nabla_{i}b+\bar{b}_{i}.$ (47)

The first term on the RHS is curl-free and the second term is divergence-free,
$\nabla_{i}\bar{b}_{i}=0$. To extract this latter part, we need to solve for
$b$; we get the following condition on $b$,

$\nabla^{2}b=\nabla_{i}b_{i}(\mathbf{x}).$ (48)

This corresponds to a Poisson’s equation, whose general solution in three
dimensions is

$b=-\int
dV^{\prime}\,\frac{\nabla^{\prime}_{i}b_{i}(\mathbf{x}^{\prime})}{4\pi|\mathbf{x}-\mathbf{x}^{\prime}|},$
(49)

where the prime on the $\nabla$ indicates that the derivative is with respect
to the components of $\mathbf{x}^{\prime}$; we thus get,

$\bar{b}_{i}=b_{i}+\nabla_{i}\int
dV^{\prime}\,\frac{\nabla^{\prime}_{j}b_{j}(\mathbf{x}^{\prime})}{4\pi|\mathbf{x}-\mathbf{x}^{\prime}|}.$
(50)

Analogous to the case for vectors, we have SVT decomposition for tensors – any
symmetric rank-2 tensor $\lambda_{ij}(\mathbf{x})$ can be decomposed as in Eq.
(1). Taking the trace of Eq. (1), we get,

$L=\frac{1}{d}\lambda_{ii},$ (51)

which allows us to rewrite Eq. (1), using the condensed notation
$\lambda_{\langle ij\rangle}\equiv\lambda_{ij}-\delta_{ij}\lambda_{ss}/d$,

$\lambda_{\langle ij\rangle}=\nabla_{\langle
i}\nabla_{j\rangle}\lambda+\nabla_{(i}\bar{\lambda}_{j)}+\bar{\lambda}_{ij}.$
(52)

Acting on this with $\nabla_{i}$ and then with $\nabla_{j}$, we get,

$\nabla_{i}\nabla_{j}\lambda_{\langle
ij\rangle}=\frac{d-1}{d}\nabla^{2}\big{(}\nabla^{2}\lambda\big{)},$ (53)

and thus, we can solve for $\lambda$, by solving the Poisson’s Eq. (53) first
for $\nabla^{2}\lambda$, and then for $\lambda$ in $d$ dimensions111For notes
on deriving the Green’s function of $\nabla^{2}$ in $d$ dimensions see, for
example, the discussion in https://tinyurl.com/knnm7bbc.,

$\lambda=\frac{d}{d-1}\frac{\Gamma^{2}\left({\frac{d}{2}}-1\right)}{16\pi^{d}}\int
dV^{\prime}\int
dV^{\prime\prime}\frac{\nabla^{\prime\prime}_{i}\nabla^{\prime\prime}_{j}\lambda_{\langle
ij\rangle}(\mathbf{x}^{\prime\prime})}{|\mathbf{x}-\mathbf{x}^{\prime}|^{d-2}|\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}|^{d-2}}.\vspace*{1ex}$
(54)

To extract $\bar{\lambda}_{i}$, we insert this value of $\lambda$ into Eq. (1)
and act on it with $\nabla_{j}$,

$\nabla^{2}\bar{\lambda}_{i}=2\nabla_{j}\lambda_{\langle
ij\rangle}-\frac{2(d-1)}{d}\nabla_{i}(\nabla^{2}\lambda\big{)},$ (55)

where $\lambda$ is understood to be known from Eq. (54). Solving this,

$\vspace*{1ex}\bar{\lambda}_{i}=-\frac{\Gamma\left({\frac{d}{2}}-1\right)}{2\pi^{\frac{d}{2}}}\int
dV^{\prime}\frac{\nabla^{\prime}_{j}\lambda_{\langle
ij\rangle}(\mathbf{x}^{\prime})-(d-1)\nabla^{\prime}_{i}\nabla^{\prime
2}\lambda(\mathbf{x}^{\prime})/d}{|\mathbf{x}-\mathbf{x}^{\prime}|^{d-2}}.$
(56)

Finally, the traceless, divergenceless part of the energy-momentum tensor is,

$\displaystyle\bar{\lambda}_{ij}$ $\displaystyle=\lambda_{\langle
ij\rangle}-\frac{d}{d-1}\frac{\Gamma^{2}\left({\frac{d}{2}}-1\right)}{16\pi^{d}}\nabla_{\langle
i}\nabla_{j\rangle}\int dV^{\prime}\int
dV^{\prime\prime}\,\frac{\nabla^{\prime\prime}_{m}\nabla^{\prime\prime}_{n}\lambda_{\langle
mn\rangle}(\mathbf{x}^{\prime\prime})}{|\mathbf{x}-\mathbf{x}^{\prime}|^{d-2}|\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}|^{d-2}}$
(57)
$\displaystyle+\frac{\Gamma\left({\frac{d}{2}}-1\right)}{4\pi^{\frac{d}{2}}}\int
dV^{\prime}$
$\displaystyle\cdot\Bigg{[}\nabla_{j}\frac{\nabla^{\prime}_{m}\lambda_{\langle
im\rangle}(\mathbf{x}^{\prime})+\frac{\Gamma^{2}\left({\frac{d}{2}}-1\right)}{16\pi^{d}}\nabla^{\prime}_{i}\nabla^{\prime
2}\displaystyle{\int dV^{\prime\prime}\int
dV^{\prime\prime\prime}\,\frac{\nabla^{\prime\prime\prime}_{m}\nabla^{\prime\prime\prime}_{n}\lambda_{\langle
mn\rangle}(\mathbf{x}^{\prime\prime\prime})}{|\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}|^{d-2}|\mathbf{x}^{\prime\prime}-\mathbf{x}^{\prime\prime\prime}|^{d-2}}}}{|\mathbf{x}-\mathbf{x}^{\prime}|^{d-2}}$
$\displaystyle+(i\longleftrightarrow j)\Bigg{]}.$

##  Appendix B An example

It is instructive to consider an example for the formalism of Sec. 2.2. For a
Beltrami field, $\mathbf{B}=(0,\sin kx,\cos kx)$, we have

$\tau_{ij}={\textstyle\frac{1}{2}}\pmatrix{-1&0&0\cr 0&-\cos 2kx&\sin 2kx\cr
0&\sin 2kx&\cos 2kx}.$ (58)

One can easily verify that

$\bar{\tau}_{ij}={\textstyle\frac{1}{2}}\pmatrix{0&0&0\cr 0&-\cos 2kx&\sin
2kx\cr 0&\sin 2kx&\cos 2kx}$ (59)

is traceless and transverse. Using $\tau=(y^{2}+z^{2})/4$, we have
$\nabla_{\langle i}\nabla_{j\rangle}\tau={\rm diag}(-2,1,1)/6$. Next, using
$T=-1/6$, we see that $T\delta_{ij}+\nabla_{\langle
i}\nabla_{j\rangle}\tau={\rm diag}(-1/2,0,0)$. Thus, following Eq. (10), we
can write $\tau_{ij}$ as $T\delta_{ij}+\nabla_{\langle
i}\nabla_{j\rangle}\tau+\nabla_{(i}\bar{\tau}_{j)}+\bar{\tau}_{ij}$ with
$\bar{\tau}_{j}=0$. The Beltrami field has therefore only a spatially periodic
tensor mode, a spatially (unbounded) nonperiodic scalar mode, and no vector
mode.

##  Appendix C The SVT decomposition of the fluid stress tensor

In this section, we provide the expressions for the scalar, vector, and tensor
components of the plasma velocity stress tensor of Eq. (19). As before, one
gets the trace-term to be

$S=-p+\frac{1}{3}w\gamma^{2}u^{2},$ (60)

and we can then solve for $\sigma$,

$\sigma=\frac{3}{32\pi^{2}}\int dV^{\prime}\int
dV^{\prime\prime}\,\frac{\nabla^{\prime\prime}_{i}\nabla^{\prime\prime}_{j}\left\\{\left(p+\rho\right)(\mathbf{x}^{\prime\prime})\gamma^{2}(\mathbf{x}^{\prime\prime})u_{\langle
i}u_{j\rangle}(\mathbf{x}^{\prime\prime})\right\\}}{|\mathbf{x}-\mathbf{x}^{\prime}||\mathbf{x}^{\prime}-\mathbf{x}^{\prime\prime}|}.$
(61)

The vector part is also extracted in a straightforward manner,

$\bar{\sigma}_{i}=-\frac{1}{2\pi}\int
dV^{\prime}\frac{\nabla^{\prime}_{j}\left\\{\left(p+\rho\right)(\mathbf{x}^{\prime})\gamma^{2}(\mathbf{x}^{\prime})u_{\langle
i}u_{j\rangle}(\mathbf{x}^{\prime})\right\\}-2\nabla^{\prime}_{i}\nabla^{\prime
2}\sigma(\mathbf{x}^{\prime})/3}{|\mathbf{x}-\mathbf{x}^{\prime}|}.$ (62)

The tensor part finally comes out by subtracting the scalar and vector parts
from the full tensor,

$\bar{\sigma}_{ij}=w\gamma^{2}u_{\langle i}u_{j\rangle}-\nabla_{\langle
i}\nabla_{j\rangle}\sigma-\nabla_{(i}\bar{\sigma}_{j)}.$ (63)

##  Appendix D The linear polarization basis

In a manner similar to the circular polarization basis, one can also define
the linear polarization basis [55],

$e^{+}_{ij}(\hat{{\bf k}})=e^{1}_{i}e^{1}_{j}-e^{2}_{i}e^{2}_{j},\quad
e^{\times}_{ij}(\hat{{\bf k}})=e^{1}_{i}e^{2}_{j}+e^{2}_{i}e^{1}_{j},\\\ $
(64)

which is orthogonal, but not normalized, and has the properties

$e^{+}_{ij}(\hat{{\bf k}})e^{+}_{ij}(\hat{{\bf k}})=e^{\times}_{ij}(\hat{{\bf
k}})e^{\times}_{ij}(\hat{{\bf k}})=2,\quad e^{+}_{ij}(\hat{{\bf
k}})e^{\times}_{ij}(\hat{{\bf k}})=0.$ (65)

In addition, the fact that $\hat{{\bf e}}^{1}$, $\hat{{\bf e}}^{2}$, and
$\hat{{\bf k}}$ form a right-handed orthonormal basis leads to the relations

$-\epsilon_{ipq}\hat{k}_{p}e^{\times}_{qj}(\hat{{\bf k}})=e^{+}_{ij}(\hat{{\bf
k}}),\quad\epsilon_{ipq}\hat{k}_{p}e^{+}_{qj}(\hat{{\bf
k}})=e^{\times}_{ij}(\hat{{\bf k}}).$ (66)

As in Eq. (34), we can decompose the stress tensor into $+$ and $\times$
polarization modes,

$\tilde{T}_{ij}(\mathbf{k},t)=\tilde{T}_{+}(\mathbf{k},t)\,e^{+}_{ij}(\hat{{\bf
k}})+\tilde{T}_{\times}(\mathbf{k},t)\,e^{\times}_{ij}(\hat{{\bf k}}).$ (67)

The projection operators of Eq. (35) can be written in this basis as

$\displaystyle\mathcal{S}_{ijlm}(\hat{\mathbf{k}})$
$\displaystyle\equiv\frac{1}{2}\left[e^{+}_{ij}e^{+}_{lm}(\hat{\mathbf{k}})+e^{\times}_{ij}e^{\times}_{lm}(\hat{\mathbf{k}})\right],$
(68) $\displaystyle\mathcal{A}_{ijlm}(\hat{\mathbf{k}})$
$\displaystyle\equiv\frac{1}{2}\left[e^{\times}_{ij}e^{+}_{lm}(\hat{\mathbf{k}})-e^{+}_{ij}e^{\times}_{lm}(\hat{\mathbf{k}})\right],$
(69)

and the degree of polarization in of the source can be expressed with Eq.
(37), where in terms of the components in the linear polarization basis,

$\displaystyle T^{\mathcal{S}}(k)$
$\displaystyle=\frac{1}{w^{2}}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int
d^{3}\mathbf{k}^{\prime}\left<\tilde{T}^{*}_{+}(\mathbf{k})\tilde{T}_{+}(\mathbf{k}^{\prime})+\tilde{T}^{*}_{\times}(\mathbf{k})\tilde{T}_{\times}(\mathbf{k}^{\prime})\right>,$
(70) $\displaystyle T^{\mathcal{A}}(k)$
$\displaystyle=\frac{1}{w^{2}}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int
d^{3}\mathbf{k}^{\prime}\left<\tilde{T}^{*}_{\times}(\mathbf{k})\tilde{T}_{+}(\mathbf{k}^{\prime})-\tilde{T}^{*}_{+}(\mathbf{k})\tilde{T}_{\times}(\mathbf{k}^{\prime})\right>,$
(71)

where the explicit $t$-dependence is dropped as in Eq. (39). Similarly, the
spectral energy densities of GWs is written in this basis as (again dropping
explicit $t$-dependence)

$\displaystyle E_{GW}(k)$ $\displaystyle=\frac{1}{(2\pi)^{6}}\frac{1}{16\pi
G}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int d^{3}\mathbf{k}^{\prime}$ (72)
$\displaystyle\times\left<\dot{\tilde{h}}^{*}_{+}(\mathbf{k})\dot{\tilde{h}}_{+}(\mathbf{k}^{\prime})+\dot{\tilde{h}}^{*}_{\times}(\mathbf{k})\dot{\tilde{h}}_{\times}(\mathbf{k}^{\prime})\right>,$
$\displaystyle H_{GW}(k)$ $\displaystyle=\frac{1}{(2\pi)^{6}}\frac{1}{16\pi
G}\int_{\Omega_{k}}d\Omega_{k}\,k^{2}\,\int d^{3}\mathbf{k}^{\prime}$ (73)
$\displaystyle\times\left<\dot{\tilde{h}}^{*}_{\times}(\mathbf{k})\dot{\tilde{h}}_{+}(\mathbf{k}^{\prime})-\dot{\tilde{h}}^{*}_{+}(\mathbf{k})\dot{\tilde{h}}_{\times}(\mathbf{k}^{\prime})\right>.$

##  Appendix E Polarization of the $S$ and $V$ components of the stress tensor

One can define a four-dimensional power spectrum analogous to that in Eq. (30)
from the scalar and vector components of the general stress tensor
$\lambda_{ij}$ of Eq. (1). Indeed using Eq. (5), we can write

$\displaystyle\frac{1}{w^{2}}\left<k_{\langle
i}k_{j\rangle}\tilde{\lambda}^{*S}(\mathbf{k},t)\cdot k^{\prime}_{\langle
l}k^{\prime}_{m\rangle}\tilde{\lambda}^{S}(\mathbf{k}^{\prime},t+\tau)\right>$
$\displaystyle\equiv(2\pi)^{3}\delta^{3}(\mathbf{k}-\mathbf{k}^{\prime})H^{S}_{ijlm}(\mathbf{k},t,\tau),$
(74)

and

$\displaystyle\frac{1}{w^{2}}\left<k_{(i}\tilde{\lambda}^{*V}_{j)}(\mathbf{k},t)\cdot
k^{\prime}_{(l}\tilde{\lambda}^{V}_{m)}(\mathbf{k}^{\prime},t+\tau)\right>$
(75)
$\displaystyle\equiv(2\pi)^{3}\delta^{3}(\mathbf{k}-\mathbf{k}^{\prime})H^{V}_{ijlm}(\mathbf{k},t,\tau).$

One can define antisymmetric analogs of these,

$\mathscr{H}^{S,V}_{ijlm}(\mathbf{k},t,\tau)\equiv-\epsilon_{ipq}\hat{k}_{p}H^{S,V}_{qjlm}(\mathbf{k},t,\tau),$
(76)

with respective “degrees of polarization” defined as

$\mathcal{P}^{S,V}_{T}(k,t)=\frac{\displaystyle{\int_{\Omega}d\Omega\,k^{2}\,\mathscr{H}^{S,V}_{ijij}(\mathbf{k},t,0)}}{\displaystyle{\int_{\Omega}d\Omega\,k^{2}\,H^{S,V}_{ijij}(\mathbf{k},t,0)}}.$
(77)

It follows in a straightforward manner from Eq. (E) that
$\mathscr{H}^{S}_{ijij}(\mathbf{k},t,0)=0$, and there is no polarization
associated with the scalar part of the stress tensor, as would be intuitively
expected. However, $\mathscr{H}^{V}_{ijij}(\mathbf{k},t,0)$ does not
necessarily vanish, and one can associate a degree of polarization with the
vector component of the source.

## References

## References

  * [1] M. Maggiore, “Gravitational Waves. Vol. 2: Astrophysics and Cosmology,” Oxford University Press, Oxford, UK (2018).
  * [2] R. Durrer and T. Kahniashvili, “CMB anisotropies caused by gravitational waves: A Parameter study,” Helv. Phys. Acta 71, 445 (1998).
  * [3] S. Weinberg, “Damping of tensor modes in cosmology,” Phys. Rev. D 69, 023503 (2004).
  * [4] A. Mangilli, N. Bartolo, S. Matarrese and A. Riotto, “The impact of cosmic neutrinos on the gravitational-wave background,” Phys. Rev. D 78, 083517 (2008).
  * [5] M. Lattanzi, R. Benini and G. Montani, “A possible signature of cosmic neutrino decoupling in the nHz region of the spectrum of primordial gravitational waves,” Class. Quant. Grav. 27, 194008 (2010).
  * [6] B. A. Stefanek and W. W. Repko, “Analytic description of the damping of gravitational waves by free streaming neutrinos,” Phys. Rev. D 88, 083536 (2013).
  * [7] X. J. Liu, W. Zhao, Y. Zhang and Z. H. Zhu, “Detecting Relic Gravitational Waves by Pulsar Timing Arrays: Effects of Cosmic Phase Transitions and Relativistic Free-Streaming Gases,” Phys. Rev. D 93, 024031 (2016).
  * [8] K. Saikawa and S. Shirai, “Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model,” JCAP 05, 035 (2018).
  * [9] D. V. Deryagin, D. Y. Grigoriev, V. A. Rubakov, and M. V. Sazhin, “Generation of gravitational waves by the anisotropic phases in the early universe,” Mon. Not. Roy. Astron. Soc. 229, 357 (1987).
  * [10] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal and A. Roper Pol, “Circular Polarization of Gravitational Waves from Early-Universe Helical Turbulence,” Phys. Rev. Res., in press (2021).
  * [11] E. Witten, “Cosmic Separation of Phases,” Phys. Rev. D 30, 272-285 (1984).
  * [12] C. J. Hogan, “Gravitational radiation from cosmological phase transitions,” Mon. Not. Roy. Astron. Soc. 218, 629-636 (1986).
  * [13] M. Kamionkowski, A. Kosowsky and M. S. Turner, “Gravitational radiation from first order phase transitions,” Phys. Rev. D 49, 2837.
  * [14] E. Barausse, E. Berti, T. Hertog, S. A. Hughes, P. Jetzer, P. Pani, T. P. Sotiriou, N. Tamanini, H. Witek and K. Yagi, et al. “Prospects for Fundamental Physics with LISA,” Gen. Rel. Grav. 52, no.8, 81 (2020).
  * [15] S. Dimopoulos, P. W. Graham, J. M. Hogan, M. A. Kasevich and S. Rajendran, “Gravitational Wave Detection with Atom Interferometry,” Phys. Lett. B 678, 37-40 (2009).
  * [16] C. Caprini, R. Durrer and X. Siemens, “Detection of gravitational waves from the QCD phase transition with pulsar timing arrays,” Phys. Rev. D 82, 063511 (2010).
  * [17] J. A. Ellis, F. A. Jenet and M. A. McLaughlin, “Practical Methods for Continuous Gravitational Wave Detection using Pulsar Timing Data,” Astrophys. J. 753, 96 (2012).
  * [18] S. Capozziello, M. Khodadi and G. Lambiase, “The quark chemical potential of QCD phase transition and the stochastic background of gravitational waves,” Phys. Lett. B 789, 626-633 (2019).
  * [19] Z. Arzoumanian et al. [NANOGrav], “The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background,” Astrophys. J. Lett. 905, L34 (2020).
  * [20] W. Ratzinger and P. Schwaller, “Whispers from the dark side: Confronting light new physics with NANOGrav data,” [arXiv:2009.11875 [astro-ph.CO]].
  * [21] A. Neronov, A. Roper Pol, C. Caprini and D. Semikoz, “NANOGrav signal from MHD turbulence at QCD phase transition in the early universe,” Phys. Rev. D 103, L041302 (2021).
  * [22] K. T. Abe, Y. Tada and I. Ueda, “Induced gravitational waves as a cosmological probe of the sound speed during the QCD phase transition,” [arXiv:2010.06193 [astro-ph.CO]].
  * [23] N. Kitajima, J. Soda and Y. Urakawa, “Nano-Hz gravitational wave signature from axion dark matter,” [arXiv:2010.10990 [astro-ph.CO]].
  * [24] N. Ramberg and L. Visinelli, “The QCD Axion and Gravitational Waves in light of NANOGrav results,” [arXiv:2012.06882 [astro-ph.CO]].
  * [25] S. L. Li, L. Shao, P. Wu and H. Yu, “NANOGrav Signal from First-Order Confinement/Deconfinement Phase Transition in Different QCD Matters,” [arXiv:2101.08012 [astro-ph.CO]].
  * [26] M. Gorghetto, E. Hardy and H. Nicolaescu, “Observing Invisible Axions with Gravitational Waves,” [arXiv:2101.11007 [hep-ph]].
  * [27] A. Sesana, “Systematic investigation of the expected gravitational wave signal from supermassive black hole binaries in the pulsar timing band,” Mon. Not. Roy. Astron. Soc. 433, 1 (2013).
  * [28] H. Middleton, A. Sesana, S. Chen, A. Vecchio, W. Del Pozzo and P. A. Rosado, “Massive black hole binary systems and the NANOGrav 12.5 year results,” Mon. Not. Roy. Astron. Soc. Lett. 502, 1 (2021).
  * [29] C. Caprini and D. G. Figueroa, “Cosmological Backgrounds of Gravitational Waves,” Class. Quant. Grav. 35, 163001 (2018).
  * [30] R. Durrer, P. G. Ferreira and T. Kahniashvili, “Tensor microwave anisotropies from a stochastic magnetic field,” Phys. Rev. D 61, 043001 (2000).
  * [31] A. Kosowsky, A. Mack and T. Kahniashvili, “Gravitational radiation from cosmological turbulence,” Phys. Rev. D 66, 024030 (2002).
  * [32] T. Kahniashvili, G. Gogoberidze and B. Ratra, “Polarized cosmological gravitational waves from primordial helical turbulence,” Phys. Rev. Lett. 95, 151301 (2005).
  * [33] V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, “On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe,” Phys. Lett. B 155, 36 (1985).
  * [34] M. E. Shaposhnikov, “Possible Appearance of the Baryon Asymmetry of the Universe in an Electroweak Theory,” JETP Lett. 44, 465-468 (1986)
  * [35] A. G. Cohen, D. B. Kaplan and A. E. Nelson, “Weak Scale Baryogenesis,” Phys. Lett. B 245, 561-564 (1990).
  * [36] A. G. Cohen, D. B. Kaplan and A. E. Nelson, “Baryogenesis at the weak phase transition,” Nucl. Phys. B 349, 727-742 (1991).
  * [37] D. E. Morrissey and M. J. Ramsey-Musolf, “Electroweak baryogenesis,” New J. Phys. 14, 125003 (2012).
  * [38] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, “Gravitational waves from the sound of a first order phase transition,” Phys. Rev. Lett. 112, 041301 (2014).
  * [39] M. Hindmarsh and M. Hijazi, “Gravitational waves from first order cosmological phase transitions in the Sound Shell Model,” JCAP 12, 062 (2019).
  * [40] M. Christensson, M. Hindmarsh and A. Brandenburg, “Inverse cascade in decaying 3-D magnetohydrodynamic turbulence,” Phys. Rev. E 64, 056405 (2001).
  * [41] T. Kahniashvili, A. Brandenburg, A. G. Tevzadze and B. Ratra, “Numerical simulations of the decay of primordial magnetic turbulence,” Phys. Rev. D 81, 123002 (2010).
  * [42] E. Lifshitz, “Republication of: On the gravitational stability of the expanding universe,” J. Phys. (USSR) 10, no. 2, 116 (1946) [Gen. Rel. Grav. 49, no. 2, 18 (2017)].
  * [43] L. P. Grishchuk, “Amplification of gravitational waves in an isotropic universe,” Sov. Phys. JETP 40, 409 (1975) [Zh. Eksp. Teor. Fiz. 67, 825 (1974)].
  * [44] V. F. Mukhanov, H. A. Feldman and R. H. Brandenberger, “Theory of cosmological perturbations. Part 1. Classical perturbations. Part 2. Quantum theory of perturbations. Part 3. Extensions,” Phys. Rept. 215, 203 (1992).
  * [45] R. Durrer and M. Kunz, “Cosmic microwave background anisotropies from scaling seeds: Generic properties of the correlation functions,” Phys. Rev. D 57, R3199-R3203 (1998).
  * [46] A. Brandenburg & Å. Nordlund, “Astrophysical turbulence modeling,” Rep. Prog. Phys. 74, 046901 (2011).
  * [47] A. Brandenburg and T. Kahniashvili, “Classes of hydrodynamic and magnetohydrodynamic turbulent decay,” Phys. Rev. Lett. 118, no.5, 055102 (2017).
  * [48] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky, “Numerical Simulations of Gravitational Waves from Early-Universe Turbulence,” Phys. Rev. D 102, 083512 (2020).
  * [49] A. J. Mee and A. Brandenburg, “Numerical turbulence forced through localized random expansion waves,” Mon. Not. Roy. Astron. Soc. 370, 415 (2006).
  * [50] B. B. Kadomtsev and V. I. Petviashili, “Acoustic turbulence,” Sov. Phys. Dokl. 18, 115 (1973).
  * [51] G. Gogoberidze, T. Kahniashvili and A. Kosowsky, “The Spectrum of Gravitational Radiation from Primordial Turbulence,” Phys. Rev. D 76, 083002 (2007).
  * [52] T. Kahniashvili, G. Gogoberidze and B. Ratra, “Gravitational Radiation from Primordial Helical MHD Turbulence,” Phys. Rev. Lett. 100, 231301 (2008).
  * [53] A. Brandenburg and S. Boldyrev, “The turbulent stress spectrum in the inertial and subinertial ranges,” Astrophys. J. 892, 80 (2020)
  * [54] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky and S. Mandal, “The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence,” Geophys. Astrophys. Fluid Dynamics 114, 130 (2020).
  * [55] C. Caprini, R. Durrer and T. Kahniashvili, “The Cosmic microwave background and helical magnetic fields: The Tensor mode,” Phys. Rev. D 69, 063006 (2004).
  * [56] A. S. Monin and A. M. Yaglom. Statistical Fluid mechanics mechanics of turbulence, vol. 2. MIT press, Cambridge, 1971.
  * [57] G. Falkovich, “Bottleneck phenomenon in developed turbulence,” Phys. Fluids 6, 1411 (1994).
  * [58] R. F. Stein, “Generation of Acoustic and Gravity Waves by Turbulence in an Isothermal Stratified Atmosphere,” Solar Physics 2, 385-432 (1967).
  * [59] W. Dobler, N. E. L. Haugen, T. A. Yousef, and A. Brandenburg, “Bottleneck effect in three-dimensional turbulence simulations,” Phys. Rev. E 68, 026304 (2003).
  * [60] A. Brandenburg, A. Johansen, P. A. Bourdin, W. Dobler, W. Lyra, M. Rheinhardt, S. Bingert, N. E. L. Haugen, A. Mee, F. Gent, N. Babkovskaia, C.-C. Yang, T. Heinemann, B. Dintrans, D. Mitra, S. Candelaresi, J. Warnecke, P. J. Käpylä, A. Schreiber, P. Chatterjee, M. J. Käpylä, X.-Y. Li, J. Krüger, J. R. Aarnes, G. R. Sarson, J. S. Oishi, J. Schober, R. Plasson, C. Sandin, E. Karchniwy, L. F. S. Rodrigues, A. Hubbard, G. Guerrero, A. Snodin, I. R. Losada, J. Pekkilä, and C. Qian, “The Pencil Code, a modular MPI code for partial differential equations and particles: multipurpose and multiuser-maintained,” Journal of Open Source Software 6, 2807 (2021). [arXiv:2009.08231 [astro-ph.IM]].
  * [61] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal, and A. Roper Pol, Datasets for “Circular Polarization of Gravitational Waves from Early-Universe Helical Turbulence” (v2020.11.07), https://doi.org/10.5281/zenodo.4256906; see also http://www.nordita.org/~brandenb/projects/CircPol/ for easier access.
  * [62] A. Brandenburg, E. Clarke, Y. He, T. Kahniashvili, Datasets for “Can we observe the QCD phase transition-generated gravitational waves through pulsar timing arrays?” (v2021.02.24), https://doi.org/10.5281/zenodo.4560423; see also http://www.nordita.org/~brandenb/projects/GWfromQCD/ for easier access.