The Turbulent Stress Spectrum in the Inertial and Subinertial Ranges

In aeroacoustics, the stress tensor of the turbulent velocity field plays an important role in sound generation. Its theory goes back to the work of Lighthill (1952a, 1952b), whose equation is also used in astrophysics to describe the heating of stellar coronae by pressure waves excited in the outer convection zones of stars (Stein 1967). Similarly, in the early universe, the velocity stress and also the combined stress of velocity and magnetic fields can be responsible for driving primordial gravitational waves (Kamionkowski et al. 1994; Durrer et al. 2000). In that case, it is important to relate spectra of the turbulence to the spectra of the kinetic and magnetic stresses in order to compute the spectrum of the gravitational waves (Gogoberidze et al. 2007; Roper Pol et al. 2019).

Empirically, it was known that a velocity or magnetic field with a Kolmogorov-type power-law spectrum produces a similar spectrum for the stress, except that in the subinertial range, where the spectral energy increases with wavenumber k, the spectral slope of the stress never increases with k faster than for white noise (Roper Pol et al. 2019), even if the turbulence has a blue spectrum. This has important implications for understanding the gravitational wave production at very low frequencies from primordial magnetic fields. Such magnetic fields can be generated at the electroweak phase transition (see Subramanian 2016, for a review), but their spectrum would be steeper than that of white noise (Durrer & Caprini 2003) and could not readily explain the shallower white noise spectrum of the stress.

There are different conventions for expressing energy spectra. In this paper, we always present the energy per uniform (linear as opposed to logarithmic) wavenumber interval, so the mean energy density is therefore ${\int }_{0}^{\infty }E(k){\text{}}{dk}$. In three dimensions, a white noise spectrum is then proportional to k². At some wavenumber k_*, the spectral energy begins to decline again. The value of k_* determines the scale where most of the energy resides. At an even higher wavenumber k_D, dissipation becomes important and the spectral energy falls off exponentially. The spectral range from k_* to k_D is called the inertial range. Its spectral slope is determined by the nature of turbulence. For Kolmogorov turbulence, it would be proportional to k^−5/3. The spectral range below k_* is called the subinertial range. Here, the flow tends to be completely uncorrelated, and this is what determines its spectral slope.

In the early universe, when it was just ${10}^{-11}\,{\rm{s}}$ old, magnetic fields are believed to have been produced with a blue subinertial range spectrum proportional to k⁴ (Durrer & Caprini 2003). This is because the magnetic field is divergence free, so the magnetic field itself does not have a white noise spectrum, but it must be the magnetic vector potential that does. Since the magnetic field is the curl of the vector potential, the spectrum of magnetic energy has an extra k² factor as compared to white noise, which is the reason why the magnetic energy spectrum is steeper than that of white noise.

There are other applications where the knowledge of the spectrum of a squared function is important. An example is the magnetic pressure, which can lead to a modulation of the gas pressure and the gas density in the interstellar medium and hence to interstellar scintillation (Lithwick & Goldreich 2001). Similarly, the square of the magnetic field perpendicular to the line of sight affects dust polarization as well as synchrotron radiation. Both dust and synchrotron emission, as well as interstellar scintillation can provide useful turbulence diagnostics in astrophysics, provided we understand the relationship between the spectra of the magnetic field and its square.

The purpose of the present paper is to derive the relationship between the spectrum of the turbulence and that of the resulting stress. Our calculations are independent of the physical model of the turbulence and apply equally to fluid and magnetohydrodynamic turbulence. With the help of several examples, we illustrate the detailed crossover behavior between different power laws. In all cases, we ignore the temporal evolution of the fluctuations. The temporal correlations are important for the radiation produced by turbulence, e.g., the gravitational waves (Gogoberidze et al. 2007), where the turbulent stress tensor enters as a source in the wave equation. Studying this in detail will be the subject of a separate investigation. Here we focus instead on the specific relationships between the spectra of a field and that of its stress found in the numerical simulations of Roper Pol et al. (2019). To illustrate the nature of the problem, it is useful to begin with a simple example of a one-dimensional scalar field and turn then to three-dimensional cases for scalar and vector fields, with and without helicity. The calculations are relatively straightforward, but we are not aware of earlier work addressing this question.

Let us consider the fluctuations of a scalar field (e.g., temperature, chemical concentration, etc.) θ(x) as a function of position x. We write θ(x) in terms of its Fourier transform as

$$\begin{eqnarray}&&\theta (x)=\int \tilde{\theta }(k){e}^{{ikx}}\displaystyle \frac{{\text{}}{dk}}{2\pi }.\end{eqnarray} \tag{ 1 }$$

Due to spatial homogeneity, the correlation function of the field can be written as

$$\begin{eqnarray}&&\langle \tilde{\theta }(k){\tilde{\theta }}^{* }(k^{\prime} )\rangle =2\pi E(k)\delta (k-k^{\prime} ),\end{eqnarray} \tag{ 2 }$$

where E(k) is the energy spectrum of θ. Its Fourier transform yields the two-point correlation function,

$$\begin{eqnarray}&&\langle \theta (x)\theta (x^{\prime} )\rangle =\int E(k){e}^{{ik}(x-x^{\prime} )}\displaystyle \frac{{\text{}}{dk}}{2\pi },\end{eqnarray} \tag{ 3 }$$

and therefore

$$\begin{eqnarray}&&\langle {\theta }^{2}(x)\rangle =\int E(k)\displaystyle \frac{{\text{}}{dk}}{2\pi }.\end{eqnarray} \tag{ 4 }$$

Consider now the fluctuations of the squared field

$$\begin{eqnarray}&&\phi (x)={\theta }^{2}(x).\end{eqnarray} \tag{ 5 }$$

We are interested in the two-point correlation function of ϕ(x). We now make an important simplifying assumption (for which a physical justification will be provided later) that the four-point correlation function of θ can be split into two-point correlation functions analogously to the Gaussian rule. We then obtain

$$\begin{eqnarray}&&\langle \phi (x)\phi (x^{\prime} )\rangle =\langle {\theta }^{2}{\rangle }^{2}+2\langle \theta (x)\theta (x^{\prime} ){\rangle }^{2}.\end{eqnarray} \tag{ 6 }$$

In order to find the energy spectrum of ϕ, we Fourier transform Equation (6) to obtain

$$\begin{eqnarray}&&\langle \phi (x)\phi (x^{\prime} )\rangle =\int F(k){e}^{{ik}(x-x^{\prime} )}\displaystyle \frac{{\text{}}{dk}}{2\pi },\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&F(k)=2\pi \langle {\theta }^{2}{\rangle }^{2}\delta (k)+2\int E(k-k^{\prime} )E(k^{\prime} )\displaystyle \frac{{\text{}}{dk}^{\prime} }{2\pi }.\end{eqnarray} \tag{ 8 }$$

The first term could be removed by subtracting the average of $\langle {\theta }^{2}{\rangle }^{2}$.

Let us assume we know the spectrum E(k). Our question concerns the resulting spectrum F(k). Specifically, we may think of a piecewise power law of the form E(k) ∝ k^α, where α is positive for 0 < k < k_*, and negative for ${k}_{* }\leqslant k\leqslant {k}_{{\rm{D}}}$, so that the energy is contained mostly at the scale ${k}_{* }^{-1}$, which is the outer scale of fluctuations. We expect F(k) to be asymptotically also of piecewise power-law form, F(k) ∝ k^β within a certain k-range. For k > 0 we have

$$\begin{eqnarray}&&F(k)=2{\int }_{-\infty }^{\infty }E(k^{\prime} )E(k-k^{\prime} )\displaystyle \frac{{\text{}}{dk}^{\prime} }{2\pi },\end{eqnarray} \tag{ 9 }$$

where we have highlighted the fact that the integration over k' goes from $-\infty$ to $+\infty$.

At small wavenumbers k ≪ k_*, the integral in Equation (9) is dominated by the scales $k^{\prime}$ comparable to the outer scale k_*, so we may expand

$$\begin{eqnarray}&&E(k^{\prime} -k)\approx E(k^{\prime} )-k\displaystyle \frac{{\text{}}{dE}(k^{\prime} )}{{\text{}}{dk}^{\prime} }+\displaystyle \frac{1}{2}{k}^{2}\displaystyle \frac{{\text{}}{d}^{2}E(k^{\prime} )}{{\text{}}{{dk}}^{{\prime} 2}}.\end{eqnarray} \tag{ 10 }$$

We then obtain from Equation (9) the asymptotic behavior of F(k) at small wavenumbers as $F(k)\approx {c}_{1}-{c}_{2}{k}^{2}$, where c₁ and c₂ are positive constants. This means that the spectrum F(k) is flat at small k, that is, β = 0.

In order to find the asymptotic behavior at large wavenumbers, k ≫ k_*, we note that, if the energy spectrum in this interval is E(k) ∝ k^α, and −3 < α < −1, then the correlation function of θ behaves at small scales as

$$\begin{eqnarray}&&\displaystyle \frac{\langle \theta (x)\theta (x^{\prime} )\rangle }{\langle {\theta }^{2}\rangle }\approx 1-{\left|\displaystyle \frac{x-x^{\prime} }{L}\right|}^{-\alpha -1},\end{eqnarray} \tag{ 11 }$$

where L ∼ 1/k_* is a scale comparable to the outer scale of the fluctuations. The square of this correlation function then scales as

$$\begin{eqnarray}&&\displaystyle \frac{\langle \theta (x)\theta (x^{\prime} ){\rangle }^{2}}{\langle {\theta }^{2}{\rangle }^{2}}\approx 1-2{\left|\displaystyle \frac{x-x^{\prime} }{L}\right|}^{-\alpha -1},\end{eqnarray} \tag{ 12 }$$

where we have expanded the right-hand side in the small parameter $| x-x^{\prime} | /L$. Therefore, asymptotically at large k, the spectrum F(k) ∼ k^β should scale with the same scaling exponent as the original energy spectrum E(k), that is, β = α.

Expressions (11) and (12) allow us to provide a physical motivation for splitting the fourth-order correlation functions of θ in the pair-wise ones in formula (6). For that, consider the Fourier component of the ϕ field,

$$\begin{eqnarray}&&\phi (k)=\int \tilde{\theta }(k^{\prime} )\tilde{\theta }(k-k^{\prime} )\displaystyle \frac{{\text{}}{dk}^{\prime} }{2\pi }.\end{eqnarray} \tag{ 13 }$$

One can ask what typical wavenumbers $k^{\prime}$ and $k-k^{\prime}$ contribute to this integral. The first possibility would be to have both wavenumbers of the same order, $k^{\prime} \sim k-k^{\prime} \sim k/2$. The second possibility is to have one of these numbers much larger than the other one, say $k^{\prime} \approx k$ and $k-k^{\prime} \approx 0$. Since for the Kolmogorov spectrum, the intensity of fluctuations declines rapidly with wavenumber, the dominant contribution is expected to come from the second possibility. This means that the fluctuating fields $\tilde{\theta }(k)$ contributing to the integral (13) have rather disparate wavenumbers. A similar argument also applies to the integral (8). Since the assumption of locality of Kolmogorov turbulence implies that the small-scale fluctuations are uncorrelated from the large-scale ones, we may average the $\tilde{\theta }$ fields in ϕ independently, which formally leads to the "Gaussian splitting rule" resulting in Equation (6).

To illustrate the resulting slope of F(k) for given E(k), we consider as a first example

$$\begin{eqnarray}E(k)=\left\{\begin{array}{ll}{k}^{-2} & {\rm{for}}\ 1\,\leqslant \,k\,\leqslant \,100\\ 0 & {\rm{otherwise}}\end{array}\right..\end{eqnarray} \tag{ 14 }$$

We compute the convolution in Equation (9) through multiplication of the Fourier transformed quantities, i.e., through their autocorrelation function,

$$\begin{eqnarray}&&\tilde{F}(x)={\left|\tilde{E}(x)\right|}^{2},\end{eqnarray} \tag{ 15 }$$

where $\tilde{E}(x)=\int {e}^{{ikx}}E(k){\text{}}{dk}/2\pi$, and likewise for $\tilde{F}(x)$. Note that the Fourier integral is carried out from $-\infty$ to $+\infty$ and that E(k) is symmetric about k = 0. We evaluate the integral in Equation (9) numerically. The energy-carrying wavenumber in our example is ${k}_{* }\equiv 1$.

The result is plotted in Figure 1. We see that in the range 5 < k < 100, we have F(k) ≈ E(k). At k = 2, the profile of F(k) has a sharp dip, which results in a local maximum at k ≈ 3, before approaching E(k). At small values of k, however, F(k) always has a flat spectrum.

Zoom In Zoom Out Reset image size

In our second example, we define the spectrum as

$$\begin{eqnarray}&&E(k)=\displaystyle \frac{| k{| }^{{\alpha }_{1}}}{{\left[1+| k{| }^{({\alpha }_{1}-{\alpha }_{2})q}\right]}^{1/q}}{e}^{-{\left(| k| /{k}_{{\rm{D}}}\right)}^{2}},\end{eqnarray} \tag{ 16 }$$

where we allow for an exponential cutoff at the dissipation wavenumber k_D = 100 and a power-law subinertial range of the form ${k}^{{\alpha }_{1}}$ with α₁ = 10, 4, 2, 1, and 0. Here we have included the exponents q and 1/q with q = 4 to sharpen the crossover between the subinertial range and the inertial range power laws. Furthermore, α₂ = −5/3 is the exponent chosen for the Kolmogorov-type inertial range.

The result is plotted in Figure 2. In all these cases, we see that F(k) has a flat spectrum for k ≲ 0.5. We see that the dip at k = 2 has disappeared for α₁ = 0, but it becomes stronger when α₁ is large and β₁ develops a white noise spectrum for small k.

Zoom In Zoom Out Reset image size

Next, we consider fully three-dimensional examples. In the case of a scalar field, $\theta ({\boldsymbol{r}})$, the derivation is similar to the one-dimensional case. The Fourier transform of the field is defined as

$$\begin{eqnarray}&&\theta ({\boldsymbol{r}})=\int \tilde{\theta }({\boldsymbol{k}}){e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{r}}}\displaystyle \frac{{\text{}}{d}^{3}k}{{\left(2\pi \right)}^{3}}.\end{eqnarray} \tag{ 17 }$$

Then, given the correlation function of the fields

$$\begin{eqnarray}&&\langle \tilde{\theta }({\boldsymbol{k}}){\tilde{\theta }}^{* }({\boldsymbol{k}}^{\prime} )\rangle ={\left(2\pi \right)}^{3}I({\boldsymbol{k}})\delta ({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ),\end{eqnarray} \tag{ 18 }$$

one derives the correlation function of the quadratic field $\phi ({\boldsymbol{r}})={\theta }^{2}({\boldsymbol{r}})$ in the form

$$\begin{eqnarray}&&\langle \phi ({\boldsymbol{r}})\phi ({\boldsymbol{r}}^{\prime} )\rangle =\int H({\boldsymbol{k}})\displaystyle \frac{{\text{}}{d}^{3}k}{{\left(2\pi \right)}^{3}},\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}&&H({\boldsymbol{k}})={\left(2\pi \right)}^{3}\langle {\theta }^{2}{\rangle }^{2}\delta ({\boldsymbol{k}})+2\int I({\boldsymbol{k}}^{\prime} )I({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )\displaystyle \frac{{\text{}}{d}^{3}k^{\prime} }{{\left(2\pi \right)}^{3}}.\quad \end{eqnarray} \tag{ 20 }$$

3.1. Nonhelical Vector Fields

The situation is qualitatively similar for a vector field. Let us consider an incompressible vector field ${\boldsymbol{u}}({\boldsymbol{r}})$, representing a velocity or magnetic field. Its Fourier transform is defined as

$$\begin{eqnarray}&&{\boldsymbol{u}}({\boldsymbol{r}})=\int \tilde{{\boldsymbol{u}}}({\boldsymbol{k}}){e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{r}}}\displaystyle \frac{{\text{}}{d}^{3}k}{{\left(2\pi \right)}^{3}}.\end{eqnarray} \tag{ 21 }$$

We assume that the distribution of this field is homogeneous and isotropic, so that its correlation function is given by

$$\begin{eqnarray}&&\langle {\tilde{u}}^{i}({\boldsymbol{k}}){\tilde{u}}^{* j}({\boldsymbol{k}}^{\prime} )\rangle ={\left(2\pi \right)}^{3}I({\boldsymbol{k}})\delta ({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ){P}_{{ij}}({\boldsymbol{k}}),\end{eqnarray} \tag{ 22 }$$

where we have denoted ${P}_{{ij}}({\boldsymbol{k}})={\delta }_{{ij}}-{k}_{i}{k}_{j}/{k}^{2}$. The energy of this field then satisfies

$$\begin{eqnarray}&&\langle {u}^{2}({\boldsymbol{r}})\rangle =\int 2I({\boldsymbol{k}})\displaystyle \frac{{\text{}}{d}^{3}k}{(2{\pi }^{3})}.\end{eqnarray} \tag{ 23 }$$

Similarly to the one-dimensional case, we are interested in the correlation function of the quadratic field ${\phi }^{{ij}}({\boldsymbol{r}})={u}^{i}({\boldsymbol{r}}){u}^{j}({\boldsymbol{r}})$. Assuming that the four-point correlation functions of the u-field can be split into the two-point ones by using the Gaussian rule, we get

$$\begin{eqnarray}&&\begin{array}{l}\langle {\phi }^{{ij}}({\boldsymbol{k}}){\phi }^{* {lm}}(\tilde{{\boldsymbol{k}}})\rangle \\ =\,\delta ({\boldsymbol{k}})\delta (\tilde{{\boldsymbol{k}}})\displaystyle \int I({\boldsymbol{k}}^{\prime} )I({\boldsymbol{k}}^{\prime\prime} ){P}_{{ij}}({\boldsymbol{k}}^{\prime} ){P}_{{lm}}({\boldsymbol{k}}^{\prime\prime} ){\text{}}{d}^{3}k^{\prime} {\text{}}{d}^{3}k^{\prime\prime} \\ \,\,+\,\delta ({\boldsymbol{k}}-\tilde{{\boldsymbol{k}}})\displaystyle \int I({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )I({\boldsymbol{k}}^{\prime} ){P}_{{il}}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ){P}_{{jm}}({\boldsymbol{k}}^{\prime} ){\text{}}{d}^{3}k^{\prime} \\ \,\,+\,\delta ({\boldsymbol{k}}-\tilde{{\boldsymbol{k}}})\displaystyle \int I({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )I({\boldsymbol{k}}^{\prime} ){P}_{{im}}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ){P}_{{jl}}({\boldsymbol{k}}^{\prime} ){\text{}}{d}^{3}k^{\prime} .\end{array}\end{eqnarray} \tag{ 24 }$$

As an example, consider the correlation function of energy fluctuations,

$$\begin{eqnarray}&&\displaystyle \frac{\langle {\phi }^{{ii}}({\boldsymbol{k}}){\phi }^{* {ll}}(\tilde{{\boldsymbol{k}}})\rangle }{{\left(2\pi \right)}^{3}}=\delta ({\boldsymbol{k}})\langle u{\left({\boldsymbol{r}}\right)}^{2}{\rangle }^{2}+2\delta ({\boldsymbol{k}}-\tilde{{\boldsymbol{k}}})H({\boldsymbol{k}}),\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&H({\boldsymbol{k}})=\int I({\boldsymbol{k}}^{\prime} )I({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )\left[1+\displaystyle \frac{{\left[{\boldsymbol{k}}^{\prime} \cdot ({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )\right]}^{2}}{{k}^{{\prime} 2}{\left({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} \right)}^{2}}\right]\displaystyle \frac{{\text{}}{d}^{3}k^{\prime} }{{\left(2\pi \right)}^{3}}.\end{eqnarray} \tag{ 26 }$$

In a statistically isotropic case, instead of the power spectrum I(k), it is convenient to use the energy spectrum $E(k)=4\pi {k}^{2}I(k)$, and similarly $F(k)=4\pi {k}^{2}H(k)$. The above equation then becomes

$$\begin{eqnarray}&&F(k)=\int E(k^{\prime} )E(\kappa )\displaystyle \frac{{k}^{2}}{{\kappa }^{2}}\left[1+\displaystyle \frac{{\left(k\mu -k^{\prime} \right)}^{2}}{{\kappa }^{2}}\right]\displaystyle \frac{{\text{}}{dk}^{\prime} \,{\text{}}d\mu }{{\left(2\pi \right)}^{3}},\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&\kappa =| {\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} | =\sqrt{{k}^{2}+{k}^{{\prime} 2}-2{kk}^{\prime} \mu },\end{eqnarray} \tag{ 28 }$$

and $\mu ={\boldsymbol{k}}\cdot {\boldsymbol{k}}^{\prime} /{kk}^{\prime}$ is the cosine of the angle between ${\boldsymbol{k}}$ and ${\boldsymbol{k}}^{\prime}$.

3.2. Helical Vector Fields

For completeness, we also consider a more general case, when the system is not mirror invariant. In this case, the field correlation function has an extra term,

$$\begin{eqnarray}\begin{array}{rcl}\langle {\tilde{u}}^{i}({\boldsymbol{k}}){\tilde{u}}^{* j}({\boldsymbol{k}}^{\prime} )\rangle & = & {\left(2\pi \right)}^{3}I(k)\delta ({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ){P}_{{ij}}({\boldsymbol{k}})\\ & - & {\left(2\pi \right)}^{3}J(k)\delta ({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ){\epsilon }_{{ijl}}{{ik}}^{l}/k.\end{array}\end{eqnarray} \tag{ 29 }$$

The last term in this expression is responsible for the helicity of the field. In particular, it enters the helicity integral

$$\begin{eqnarray}&&\int \langle ({\boldsymbol{\nabla }}\times {\boldsymbol{u}})\cdot {\boldsymbol{u}}\rangle \,{\text{}}{d}^{3}x=2\int k\,J(k)\displaystyle \frac{{\text{}}{d}^{3}k}{{\left(2\pi \right)}^{3}},\end{eqnarray} \tag{ 30 }$$

which would be zero in a mirror-invariant case. The helicity spectral function J(k) is not necessarily positive, but it has to satisfy the realizability condition $| J(k)| \leqslant I(k)$. The generalization of our results to the helical case is straightforward; it is achieved by replacing the nonhelical terms $I(k){P}_{{ij}}({\boldsymbol{k}})$ in Equation (24) by their helical counterparts $I(k){P}_{{ij}}({\boldsymbol{k}})-J(k){\epsilon }_{{ijl}}{{ik}}^{l}/k$. For example, in the correlation function of the energy fluctuations (25) we will need to replace the $H({\boldsymbol{k}})$ function (26) by the general expression

$$\begin{eqnarray}\begin{array}{rcl}H({\boldsymbol{k}}) & = & \displaystyle \int I({\boldsymbol{k}}^{\prime} )I({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )\left[1+\displaystyle \frac{{\left[{\boldsymbol{k}}^{\prime} \cdot ({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )\right]}^{2}}{{k}^{{\prime} 2}{\left({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} \right)}^{2}}\right]\displaystyle \frac{{\text{}}{d}^{3}k^{\prime} }{{\left(2\pi \right)}^{3}}\\ & - & 2\displaystyle \int J({\boldsymbol{k}}^{\prime} )J({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )\left[\displaystyle \frac{{\boldsymbol{k}}^{\prime} \cdot ({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )}{k^{\prime} | {\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} | }\right]\displaystyle \frac{{\text{}}{d}^{3}k^{\prime} }{{\left(2\pi \right)}^{3}}.\end{array}\end{eqnarray} \tag{ 31 }$$

Again, in a statistically isotropic case, instead of the power spectrum J(k), it is convenient to introduce the helicity spectrum defined as $G(k)=4\pi {k}^{2}J(k)$. The above equation then becomes

$$\begin{eqnarray}\begin{array}{rcl}F(k) & = & \displaystyle \int E(k^{\prime} )E(\kappa )\displaystyle \frac{{k}^{2}}{{\kappa }^{2}}\left[1+\displaystyle \frac{{\left(k\mu -k^{\prime} \right)}^{2}}{{\kappa }^{2}}\right]\displaystyle \frac{{\text{}}{dk}^{\prime} \,{\text{}}d\mu }{{\left(2\pi \right)}^{3}}\\ & - & 2\displaystyle \int G(k^{\prime} )G(\kappa )\displaystyle \frac{{k}^{2}}{{\kappa }^{2}}\,\displaystyle \frac{k\mu -k^{\prime} }{\kappa }\,\displaystyle \frac{{\text{}}{dk}^{\prime} \,{\text{}}d\mu }{{\left(2\pi \right)}^{3}},\end{array}\end{eqnarray} \tag{ 32 }$$

which can readily be evaluated using numerical integration.

3.3. Examples

In Figure 3, we show the results for the case of a single power law, as in Equation (12). We consider two values for the slope α (−2 and −3) in the range 1 ≤ k ≤ 100. We see that, for k > 2, we obtain for F(k) a power law, F(k) ∝ k^β, with β = α, as in the one-dimensional case. In the range 1 < k < 2, F(k) is still increasing with k, but the slope is slightly less steep than two. We emphasize that this behavior is different from that in the one-dimensional case, where we saw instead a marked dip in F(k).

Zoom In Zoom Out Reset image size

In Figure 3, we also plot F(k) for the case of a fully helical field using Equation (32), where G(k) = E(k) is assumed. We see that the basic features of F(k) are rather similar to the case without helicity, but there is now slightly more power in the subinertial range, where F(k) appears to be elevated by about a factor of two. In the inertial range, on the other hand, F(k) is not affected by the presence of helicity.

Next, we consider a spectrum with two different slopes, α₁ and α₂, along with an exponential cutoff, just as in Equation (16). Again, we denote the corresponding slopes of F(k) as β₁ and β₂, respectively. Here, we always assume a Kolmogorov inertial range spectrum for E(k), i.e., α₂ = −5/3, and we vary α₁ from 0 to 10. Physically relevant are the Saffman (α₁ = 2) and Batchelor (α₁ = 4) asymptotic scalings for k → 0 (e.g., Davidson 2015). In our finite simulation domain it is, however, interesting to consider arbitrary values of α₁.

Figure 4 confirms the statement of Roper Pol et al. (2019) that β = 2 is obtained even if E(k) has a blue spectrum, i.e., α ≥ 2. We also see from Figure 4 that for α₁ = 2, the crossover from the k² scaling for small k to the k^−5/3 scaling for large k extends now over more than one decade (0.2 < k < 5). This shows that we may expect slight differences when approximate scalings are based on the inspection of spectra over a limited dynamical range.

Zoom In Zoom Out Reset image size

To study in more detail the crossover from β = 2 for α ≥ 2 to β = α for α of around and below −5/3, for example, let us now consider single power-law spectra within a more extended range 1 ≤ k ≤ 1000 using intermediate values α = −1, 0, and 1. No distinction between α₁ and α₂ will therefore be made. The result is shown in Figure 5. We see that in this range of α, β is always larger than α. We see that already for the scale-invariant k⁻¹ spectrum, we have β = −0.86, so the β = α relation is only approximately obeyed.

Zoom In Zoom Out Reset image size

To determine the relation between α and β in the intermediate regime, we now compute β using the same numerical setup as before, but we now consider power-law scalings in a range that is 100 times larger, 1 ≤ k ≤ 10⁵. The result is shown in Figure 6. Here we also compare with the corresponding results in one dimension. We now see that in the range −1 < α < 1, the value of β deviates markedly from the β = α relation, and that we have β ≈ 2 already for α = 1.

Zoom In Zoom Out Reset image size

Likewise, in one dimension, the β = α relation is only true for α < −1. Again, there is a small intermediate interval, −1 < α < 0, where β is somewhat larger than α, but this departure is by far not as dramatic as in three dimensions; see Figure 6.

To address the assumption of Gaussianity, we now present the results of numerical simulations of the hydromagnetic equations for a weakly compressible gas in a cubic domain of size L³ using 1024³ mesh points. We employ the Pencil Code,⁷ which uses sixth order accurate finite differences for the spatial discretization and a third order time stepping scheme. We first consider decaying nonhelical turbulence. Our simulations are similar to Run A of Brandenburg et al. (2017), where α₁ ≈ 4 and α₂ ≈ −5/3. The turbulence is magnetically dominated, so the velocity is almost entirely the result of the Lorentz force.

In Figure 7, we present the results for the magnetic energy spectrum ${E}_{{\rm{M}}}(k)$ after about 100 Alfvén times. Initially, the peak of ${E}_{{\rm{M}}}(k)$ was at $k={k}_{* }$ with ${k}_{* }L/2\pi =60$, but, owing to an effect similar to the inverse cascade—here without helicity—the peak has moved to about k_* L/2π = 15 by the end of the simulation; see Brandenburg et al. (2015) for similar results. In Figure 7 we also compare with the corresponding spectrum of ${{\boldsymbol{B}}}^{2}$, referred to as F_i(k), where i = M stands for the magnetic stress. We see that the spectral slopes of ${E}_{{\rm{M}}}$ and ${F}_{{\rm{M}}}$ agree in the inertial range. In the subinertial range, the slopes of ${E}_{{\rm{M}}}$ and ${F}_{{\rm{M}}}$ are, as expected, different from each other. However, the values of the slopes, 3.5 and 1.5, respectively, are below the expectations of 4 and 2, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To characterize the departure from Gaussianity, we have computed the kurtosis of the magnetic field separately for all three components and then take the average, which is denoted by

$$\begin{eqnarray}&&{\rm{kurt}}{\boldsymbol{B}}=-3+\displaystyle \frac{1}{3}\sum _{i=1}^{3}\langle {B}_{i}^{4}\rangle /\langle {B}_{i}^{2}{\rangle }^{2}.\end{eqnarray} \tag{ 33 }$$

We find a rather small value of less than 0.1. Thus, the field is close to Gaussian and our results are qualitatively in close agreement with those of the present paper.

To compare with a field where the assumption of Gaussianity cannot be justified, we also show the results for the normalized current density ${\boldsymbol{J}}={\boldsymbol{\nabla }}\times {\boldsymbol{B}}$. We denote the corresponding spectra for current density by E_i(k) and F_i(k) with i = C. The kurtosis of ${\boldsymbol{J}}$, defined analogous to ${\rm{kurt}}{\boldsymbol{B}}$, is about 4.5.

In the inertial range, ${E}_{{\rm{M}}}$ has a slope of about −2.2, which is steeper than that of the Kolmogorov spectrum. The current density spectrum has a slope of about −0.2. The slopes were initially closer to the Kolmogorov values, but they became steeper with time. What is important, however, is that in the inertial range, the spectral slopes of ${F}_{{\rm{M}}}$ and ${F}_{{\rm{C}}}$ agree with those of ${E}_{{\rm{M}}}$ and ${E}_{{\rm{C}}}$, respectively, i.e., both have slopes of −2.2 and −0.2, respectively. In the subinertial range, the slopes are again somewhat different from the expectation. We find ${\alpha }_{1}=3.5$ and 5.5 for the spectra of ${E}_{{\rm{M}}}$ and ${E}_{{\rm{C}}}$, respectively, but 1.5 for the slopes of both ${F}_{{\rm{M}}}$ and ${F}_{{\rm{C}}}$. It should be noted that, even under the assumption of isotropy, the stress tensor contains different contributions (scalar, vector, and tensor modes) that might behave differently. However, the resulting differences are also sensitive to the nature of the turbulence, whose study is beyond the scope of the present paper.

Next, we compare with two runs of forced turbulence. In these two examples, we consider the forcing wavenumbers k_* = 30 and 6, respectively. In the former case (Figure 8), the subinertial range is more developed. The magnetic field and current density are now closer to being Gaussian ($| {\rm{kurt}}{\boldsymbol{B}}| \lesssim 0.1$ and ${\rm{kurt}}{\boldsymbol{J}}\approx 1.4$). We see that in the subinertial range, the slope of ${F}_{{\rm{C}}}$ is now even more shallow (${\beta }_{1}=1$), while that of ${F}_{{\rm{M}}}$ is slightly steeper (${\beta }_{1}=1.7$), but still not quite as steep as what is expected (β₁ = 2).

In the latter run with k_* = 6 (Figure 9), the inertial range is more developed and we see a clear k^−5/3 spectrum in ${E}_{{\rm{M}}}$. The magnetic field and current density are now further away from being Gaussian (${\rm{kurt}}{\boldsymbol{B}}=0.1...0.2$ and ${\rm{kurt}}{\boldsymbol{J}}\approx 5$). In the inertial range, the slope for ${E}_{{\rm{C}}}$ agrees with what is expected for Kolmogorov-type turbulence (1/3 for ${F}_{{\rm{M}}}$). However, we also see departures from the $\beta =\alpha$ relation, where the β for ${F}_{{\rm{M}}}$ is slightly larger, while that of ${F}_{{\rm{C}}}$ is now smaller and even negative.

We have derived the general formula that allows us to compute a spectrum F(k) of the square of a fluctuating field whose spectrum, in turn, is E(k). Our results are independent of whether the spectrum is that of a scalar or that of a vector field. We have seen that in the inertial range with $E(k)\propto {k}^{\alpha }$ and α ≲ −1, we find a spectrum F(k) ∝ k^β with β ≈ α if we are sufficiently far away from the boundaries of the validity range of where the power law applies. In the subinertial range, where α ≳ 1, we find β ≈ 2.

A possible application of our work concerns the generation of gravitational waves from hydrodynamic and hydromagnetic turbulence with a known energy spectrum E(k). The resulting stress T_ij sources a wave equation of the form $\square {h}_{{ij}}={T}_{{ij}}$, where, except for normalization factors, h_ij is the linearized strain field, $\square ={{\rm{\nabla }}}^{2}-{c}^{-2}{\partial }^{2}/\partial {t}^{2}$ is the d'Alembertian wave operator, and c is the speed of light. The actual gravitational wave fields are the transverse and traceless projections of h_ij. The nature of the wave operator can lead to a complicated relation between the spectra of h_ij and T_ij (e.g., Roper Pol et al. 2020). If, however, the time delay in the wave equation can be neglected, the spectrum of h_ij is proportional to $F(k)/{k}^{4}$.

The simulations of Roper Pol et al. (2020) show that most of the wave generation occurs at the time when the stress has reached maximum amplitude. Subsequent changes of the source hardly contribute to wave production. It may be for this reason that the assumption of no time delay is a reasonable one. Under this assumption, we expect that the gravitational wave energy, which is proportional to ${(\partial {h}_{{ij}}/\partial t)}^{2}$, should be proportional to $F(k)/{k}^{2};$ see Roper Pol et al. (2020) for details. The extent of the empirically determined departures from this simplistic way of estimating the gravitational wave energy spectrum are not yet fully understood and would need to be determined numerically or analytically, similarly to earlier work using the aeroacoustic approximation of Lighthill (1952a), as already done by Gogoberidze et al. (2007) in the context of primordial gravitational waves.

Cosmological magnetic fields may well be helical (Tashiro et al. 2014). They could be generated by the chiral magnetic effect (Joyce & Shaposhnikov 1997; Boyarsky et al. 2012; Yamamoto 2016; Anand et al. 2017). However, even under the most optimistic conditions, this effect can only produce about ${10}^{-18}\,{\rm{G}}$ on a scale of $1\,\mathrm{Mpc}$ at the present time (Brandenburg et al. 2017). Our work now shows that the presence of magnetic helicity enhances the stress spectrum by up to a factor of two in the subinertial range, while leaving that in the initial range unchanged. This implies a small shift in the peak of the stress spectrum, which itself would hardly be a distinguishing feature. However, helical magnetic fields lead to circular polarization of gravitational waves (Kahniashvili et al. 2005), which may be detectable with the Laser Interferometer Space Antenna if there is a sufficiently strong dipolar anisotropy in the signal (Domcke et al. 2019).

Another potentially important application concerns the spectrum of the parity even and parity odd linear polarization modes. Those depend quadratically on the magnetic field components perpendicular to the line of sight (Caldwell et al. 2017; Kandel et al. 2017; Brandenburg et al. 2019). Our work now suggests that, measuring the polarization spectrum, one can only infer the spectrum of the underlying turbulence if $\alpha \lt -1$.

We acknowledge initial discussions with Ethan Vishniac during the program on the Turbulent life of Cosmic Baryons at the Aspen Center for Physics, which is supported by the National Science Foundation grant PHY-1607611. We also thank the anonymous referee for encouraging us to study the case with helicity. The work of A.B. was supported by the National Science Foundation under the grant AAG-1615100, and the work of S.B. was supported by the National Science Foundation under the grant PHY-1707272 and by NASA under the grant 80NSSC18K0646. S.B. was also partly supported by the DOE grant DE-SC0018266. 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.