Low frequency tail of gravitational wave spectra from hydromagnetic turbulence

Ramkishor Sharma^{1,2}^{∗} and Axel Brandenburg^{1,2,3,4}^{†}

1Nordita, KTH Royal Institute of Technology and Stockholm University,

Hannes Alfvéns väg 12, 10691 Stockholm, Sweden

2The Oskar Klein Centre, Department of Astronomy, Stockholm University, 10691 Stockholm, Sweden 3McWilliams Center for Cosmology & Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA and

4School of Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili Avenue, 0194 Tbilisi, Georgia

Hydrodynamic and magnetohydrodynamic (MHD) turbulence in the early Universe can drive gravitational waves (GWs) and imprint their spectrum onto that of GWs, which might still be observable today. We study the production of the GW background from freely decaying MHD turbulence from helical and nonhelical initial magnetic fields. To understand the produced GW spectra, we develop a simple model on the basis of the evolution of the magnetic stress tensor. We find that the GW spectra obtained in this model reproduce those obtained in numerical simulations if we consider the detailed time evolution of the low frequency tail of the stress spectrum from numerical simulations. We also show that the shapes of the produced GW frequency spectra are different for helical and nonhelical cases for the same initial magnetic energy spectra. Such differences can help distinguish helical and nonhelical initial magnetic fields from a polarized background of GWs – especially when the expected circular polarization cannot be detected directly.

I. INTRODUCTION

Magnetohydrodynamic (MHD) turbulence in the early Universe can be a powerful source of gravitational waves (GWs) that could be observable as a stochastic background today [1–5]. The frequency spectrum of these waves is related to the spectrum of the underlying turbulence. Such turbulence could be induced during the various epochs^{1} in the early Universe [8–11] or the possible presence of primordial magnetic fields [12–19]. These GWs produced by turbulence at an epoch of the electroweak phase transition lie in the sensitivity range of the proposed Laser Interferometer Space Antenna and pulsar timing arrays for the turbulence induced around an epoch of the quantum chromodynamics (QCD) phase transition. Recently, various pulsar timing arrays [20–23] have reported evidence for the presence of a common spectrum process across analyzed pulsars in the search of the presence of an isotropic stochastic GW background. This evidence has been used to constrain the strength and correlation length of magnetic fields generated at the QCD epoch [24–26]. However, the presence of a quadrupolar spatial correlation [27], a characteristics of a GW background, is yet to be claimed.

Numerical simulations have confirmed that there is indeed a direct connection between the slopes of the turbulence and GW spectra [28], except that at low frequencies, below the peak of the spectrum, the GW spectrum was found to be shallower in the simulations than what was previously expected from analytical calculations. We

∗ ramkishor.sharma@su.se
† brandenb@nordita.org

1 The electroweak and QCD epochs are accompanied by crossovers in the standard model [6, 7]. However, many extensions of the standard model can lead to a first order phase transition.

call this part the low frequency tail of the GW spectrum. However, there is the worry that this shallow tail could be caused by unknown numerical artifacts such as the finite size of the computational domain and the way the turbulence is initiated in the simulations.

To understand the origin of the low frequency tail, the authors of Ref. [26] have recently compared numerical MHD simulations with an analytic model, where the stress is assumed constant for a certain interval of time. Their model predicts a flat spectrum whose extent depends on the duration over which the stress is held constant. In this way, it was possible to determine an effective duration for a given numerical simulation. This duration was found to be different for different simulations. Their model is therefore descriptive rather than predictive. Furthermore, in the numerical solutions the stress was not actually constant for any duration of time.

In another recent approach, the authors of Ref. [29] have focused on the importance of unequal time correlation functions of the Fourier components of the velocity field for purely hydrodynamic turbulence. While the authors acknowledge the potential importance of the initial growth phase of the turbulence, they also have no inverse cascade in their simulations. This is different from MHD turbulence, which can display inverse cascading even in the absence of net magnetic helicity [30, 31]. This will be crucial to the approach discussed in the present paper.

To address the problem of a limited computational domain, it is important to use large enough computational domains so that its minimum wave number is as small as possible. In this paper, we discuss two MHD simulations where the wave number corresponding to the peak of the GW spectrum and the wave numbers below that corresponding to the horizon size at the initial time are well resolved. Since the stress appears explicitly in the linearized GW equation, we also analyze for these simulations the evolution of the stress spectrum along with the magnetic and GW spectra. Such a detailed comparison between the simulated stress and resulting GWs is an important new aspect of the present work. Second, we develop a simple model, motivated by the stress evolution seen in the present simulations, to explain the GW spectrum obtained. In this model, our main focus is to understand the nature of the GW spectrum below the wave number corresponding to the peak of the spectrum. Our simulations are similar to those of Ref. [26], but our interpretation and corresponding modeling of the stress is not. There is no time interval during which the stress is constant. Our results are therefore not characterized by the duration of such a time interval. In addition, we determine and analyze spectral differences between runs with and without magnetic helicity, which were not noticed previously. We also emphasize that the Hubble horizon wave number poses an ultimate cutoff for the flat spectrum toward low wave numbers.

This paper is organized as follows. In Section II, we discuss the evolution of the magnetic field, stress, and GW spectrum in our new runs. In this section, we also discuss how the stress spectrum evolves when inverse transfer and inverse cascade of the turbulence correspond to the evolution for nonhelical and helical magnetic fields in the early Universe. In Section III, we discuss the model to explain the low frequency tail of the GW spectrum. Further, in Section IV, we compare the GW spectrum obtained from our numerical simulations and our model. We conclude in Section V.

II. NONHELICAL AND HELICAL CASCADES

Various phenomena such as primordial magnetic fields and phase transitions can lead to the generation of turbulence in the early Universe. The stress associated with magnetic fields and turbulence lead to the production of GWs. This has been studied in the literature both analytically [4, 32–37] and numerically [26, 28, 29, 38–40]. In the present paper, we perform new simulations of decaying MHD turbulence, where we resolve the scales which are smaller than the Hubble horizon size at the initial time. Before explaining the simulations in detail, let us begin by summarizing the basic equations.

A. GWs from MHD turbulence

We follow here the formalism of Ref. [28, 41], where conformal time is normalized to unity at the initial time. One could associate this with the electroweak phase transition, for example. The velocity u is normalized to the speed of light. The magnetic field B = ∇ × A is written in terms of the magnetic vector potential A, and the current density is written as J = ∇×B. Following Ref. [42], the energy density ρ includes the restmass density, so its evolution equation obeys a continuity equation that also includes magnetic energy terms. As in [41], ρ is normalized to the critical energy density for a flat Universe. We solve for the Fourier transformed plus and cross polarizations of the gravitational strain, h̃_{+} and h̃_{×}, which are driven by the corresponding projections of the stress, which, in turn, is composed of kinetic and magnetic contributions,

T_{ij} = \frac{4}{3}γ_{Lor}^{2}ρu_{i}u_{j} − B_{i}B_{j} + ...,  (1)

where γ_{Lor} = (1 − u^{2})^{−1/2} is the Lorentz factor, and the ellipsis denotes terms proportional to δ_{ij}, which do not contribute to the projected source T̃_{+/×}.

Assuming the Universe to be conformally flat, its expansion can be scaled out by working with conformal time t and comoving variables [42]. We use the fact that in the radiation-dominated era, the scale factor grows linearly with conformal time. The only explicit occurrence of conformal time is then in the GW equation, where a 6/t factor occurs in the source term [41]. The full set of equations is therefore

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

\frac{Du}{Dt} = \frac{1}{ρ}∇ · (2ρνS) − \frac{1}{4}∇ ln ρ + \frac{u}{3} (∇ · u + u · ∇ ln ρ)

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

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

\frac{∂^{2}}{∂t^{2}} h̃_{+/×}(k, t) + k^{2}h̃_{+/×}(k, t) = \frac{6}{t}T̃_{+/×}(k, t),  (4)

where D/Dt ≡ ∂/∂t + u · ∇ is the advective derivative, η is the magnetic diffusivity, ν is the kinematic viscosity, S_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) − \frac{1}{3}δ_{ij}∇ · u are the components of the rate-of-strain tensor S with commas denoting partial derivatives. Fourier transformation in space is denoted by a tilde. In all cases studied in this paper, the initial conditions are such that B consists of a weak Gaussiandistributed seed magnetic field, u = 0, ρ = 1.

We work with spectra that are defined as integrals over concentric shells in wave number space k with k = |k|. They are normalized such that their integrals over k give the mean square of the corresponding quantity, i.e., ∫ Sp(B) dk = ⟨B^{2}⟩, where Sp(B) = Sp(B_{x}) + Sp(B_{y}) + Sp(B_{z}). Similarly, Sp(h) = Sp(h_{+})+Sp(h_{×}) is defined as the sum over the two polarization modes. Of particular interest will also be the stress spectrum Sp(T), which is defined analogously through Sp(T) = Sp(T_{+}) + Sp(T_{×}). To study the evolution of the stress at selected Fourier modes, we compute |T̃(k, t)| ≡ \sqrt{Sp(T)/4πk^{2}}, which scales the same way as |T̃_{+}(k, t)| and |T̃_{×}(k, t)|.

B. Evolution of the stress and strain spectra

To put our results into perspective and compare with earlier work, we study cases of suddenly initiated tur-

TABLE I. Summary of simulation parameters

Run  E_{M}^{i}  k_{p} E_{GW}^{sat}  Ω_{GW}^{sat}
HEL 5.4 × 10^{−3} 10 3.7 × 10^{−7} 5.9 × 10^{−12}
NHEL 5.5 × 10^{−3} 10 3.5 × 10^{−7} 5.6 × 10^{−12}

bulence. We perform simulations similar to those of Ref. [28] by using as initial condition for the magnetic field a random Gaussian-distributed magnetic field with a k^{4} spectrum for k < k_{p} and a k^{−5/3} spectrum for k > k_{p}. For details of such a magnetic field, see Ref. [43]. As initial condition for the GW field, we assume that h and ḣ vanish. The strength of the GW field is then strongly determined by the sudden initialization of a fully developed turbulence spectrum. The details of the simulations are given in Table I. In this table, the first column represents the name of the runs (HEL and NHEL), E_{EM}^{i} is the initial value of the magnetic energy density compared to the background energy density, k_{p} is the wave number at which the magnetic energy spectrum peaks and it is normalized by the wave number corresponding to the Hubble horizon size at the initial time, E_{GW}^{sat} is the value of the GW energy density after saturation compared to the background energy density, and Ω_{GW}^{sat} is the density parameter of GWs, representing the ratio of the GW energy density compared to the critical energy density at present. Ω_{GW}^{sat} has been calculated considering the production of GWs around the electroweak phase transition.

We consider two runs where the initial magnetic field is either fully helical (HEL) or nonhelical (NHEL). In Figs. 1 and 2, we show for HEL and NHEL, respectively, the time evolution of the spectra of the magnetic field, the TT-projected stress, the strain derivative, and the strain. Inverse cascading is seen in the magnetic energy spectra, which leads to the expected increase of the spectral stress at small k; see Figs. 1(a) and (b) for the HEL run. By contrast, in NHEL, the stress spectrum always decreases at small k. We also see in Fig. 1(c) that the GW energy spectrum has a maximum at k ≈ 20, which is not present in NHEL; cf. Fig. 2(c). Their spectra fall off toward smaller k proportional to k and k^{1.5} for HEL and NHEL, respectively.

In Figs. 3 and 4, we show the time evolution of the modulus and phase of stress and strain derivative for HEL and NHEL for wave vectors k = (k, 0, 0) and five different values of k, which are all below k_{p}. The purpose of this is to see how representative these individual wave vectors are compared to the collective effect of all others of similar length, and whether there is any important effect resulting from the phases of the stress.

Broadly speaking, the time evolution of the modulus of the stress |T̃| for any of the five wave vectors does not seems to reflect the expectation from the evolution of the shell-integrated stress spectrum, which is increasing for HEL and decreasing for NHEL, as was see in Figs. 1 and 2. This is an important observation and may be due to the fact that in Figs. 3 and 4, we have shown stress and strain derivatives only for particular values of k.

From Figs. 3 and 4, it is evident that arg(T̃) remains constant for some time and starts evolving more rapidly after that. It is also interesting to note that the amplitude of |˙h̃| increases up to the time until which arg(˙h̃) is roughly constant. After this time, |˙h̃| enters an oscillatory regime and its amplitude does not change much. Other wave vectors of the same length show a similar behavior of arg(T̃) and arg(˙h̃), as is shown in the figures of the Supplemental Material provided along with the data in Ref. [44]. This conclusion applies for both runs shown in Figs. 3 and 4 and it leads us to develop a simple model to understand the GW spectrum in these cases. In this model we replace T̃ by its wave vector-averaged magnitude, |T̃|, as discussed in Section III.

Further, to understand the role of the phases of the stress tensor in the production of GWs, we run two new simulations analogous to runs HEL and NHEL, where we replace T̃(k, t) by its modulus at each time step. The final GW spectrum in these modified runs turns out to be virtually the same as in the original HEL and NHEL runs; see Fig. 5. The comparisons of HEL and NHEL are shown in panels (a) and (b) of this figure, respectively. Dashed red and gray curves at times t = 1.5 and t = 37, respectively, are for the cases when T̃(k, t) has been replaced by its modulus. It is evident from the figure that there is hardly any difference in the actual Sp(ḣ) after replacing the stress with its modulus. On the basis of this observation, we develop a model to obtain the GW spectrum from the time evolution of the spectrum of the stress tensor.

A striking difference between HEL and NHEL is the more pronounced peak in the spectral GW energy in the former. As we show in Appendix A, this is due to the fact that the stress spectrum for NHEL is different from that of HEL due to the presence of additional helical contributions to the two-point correlation of the magnetic field vectors. This difference is shown in Fig. 6 and the details are explained next.

C. Overall behavior of the stress

At the most minimalistic level, we can say that the magnetic field shows an approximately self-similar evolution at late times, where for HEL, the peak value of E_{M}(k, t) is unchanged, but the position of the peak k_{p} goes to progressively smaller values as k_{p} ∼ t^{−2/3}.

To understand the consequences for the evolution of the stress, let us now consider an idealized model, where E_{M}(k) ≡ Sp(B)/2 has a k^{4} subinertial range, where k < k_{p}, with k_{p}(t) being the peak wave number, and a k^{−5/3} inertial range spectrum for k > k_{p}. The spectrum of the transverse traceless part of the stress, Sp(T)/2, can

FIG. 1. Spectra of the magnetic field, the TT-projected stress, the strain derivative, and the strain for suddenly initiated turbulence with magnetic helicity.

FIG. 2. Same as Fig. 1, but for the nonhelical case.

FIG. 3. Modulus and phase of T̃(k, t) and ˙ h̃(k, t) for the helical case for k = (k, 0, 0) with k = 0.3 (orange), 0.4 (red), 0.5 (green), 0.6 (blue), and 0.7 (black). The inset shows the phase with a linear abscissa.

FIG. 4. Same as Fig. 3, but for the nonhelical case.

FIG. 5. Sp(ḣ)(k, t) vs k: (a) The solid black and blue curves represent Sp(ḣ)(k, t) at times t = 1.5 and t = 37 for run hel. The dashed red and gray curves show Sp(ḣ)(k, t) for the case when the stress spectrum has been replaced by its modulus in the GW evolution equation. The black curve coincides with the dashed red curve and solid blue curve coincides with dashed gray curve. (b) Same as (a), but for run nonhel.

FIG. 6. In this figure, magnetic field energy spectrum, E_{M}(k) (Dashed curves) and Sp(T) (Solid curves) for the helical and non helical case. The blue and red curves are for nonhelical and helical case respectively.

be computed analytically using the expressions given in Appendix A (see Refs. [37, 45] for details) and is shown in Figs. 7 and 8 for the helical and nonhelical cases, respectively. In these figures, we take three instances where the magnetic peaks are at wave numbers k_{p} = 1, 0.3, and 0.1.

For HEL, the position of the peak of E_{M}(k) is unchanged. We see that, in agreement with earlier work [28], the positions of the peak of Sp(T) are always at 2k_{p}. However, even though the peak values of E_{M}(k) are unchanged, except for the factor of two, those of Sp(T) are not and decay. Nevertheless, at small k, Sp(T) still increases proportional to k_{p}^{−1}. If k_{p} ∝ t^{−2/3}, as expected for helical turbulence [46–48], we find that Sp(T) ∝ t^{2/3} for small k.

For NHEL, as shown in Ref. [48], the peak of the spectrum decreases with decreasing values of k_{p} proportional

TABLE II. Time evolution of Sp(T) and |T̃| from theory

helical  nonhelical
Sp(T) vs k_{p} at small k ∝ 1/k_{p}  ∝ k_{p}
k_{p} vs t  k_{p} ∝ t^{−2/3} k_{p} ∝ t^{−1/2}
Sp(T) vs t  Sp(T) ∝ t^{2/3} Sp(T) ∝ t^{−1/2}
|T̃| vs t  |T̃| ∝ t^{1/3} |T̃| ∝ t^{−1/4}

to k_{p}^{β}, where β is an exponent that can be between one and four. In Fig. 8, we present the case with β = 1 and find that now Sp(T)(k) ∝ k_{p} for small k and ∝ k_{p}^{14/3} for large k. If k_{p} ∝ t^{−1/2}, as expected for the nonhelical case for β = 1, Sp(T) ∝ t^{−1/2} for small k. We have summarized the behavior of Sp(T) with time in Table II for helical and nonhelical cases.

Recently, it has been found that the Hosking integral [49], a Saffman-like helicity integral, is well conserved in nonhelical magnetically dominated decaying turbulence [50, 51], which implies β = 1.5. For the general case, we write Sp(T) ∝ k_{p}^{2β−1} (for k < k_{p}), which implies Sp(T) ∝ t^{−8/9} for β = 1.5 and k_{p} ∝ t^{−4/9}.

It is also interesting to note that, for a given spectrum of the magnetic field, Sp(T) is different for HEL and NHEL; see the blue and red curves in Fig. 6), respectively. For the helical case, Sp(T) has smaller values compared to the nonhelical case at wave numbers below k_{p}. However, it has large values for wave numbers around the peak and above. Such a feature of the stress spectrum also translates to the GW spectrum and that is why we see a difference in the final GW spectrum produced from helical and nonhelical cases discussed in the previous section.

FIG. 7. Left: Solutions for Sp(T) (red) for different E_{M}(k) (blue) for three values of k_{p}. Right: solutions for Sp(T) scaled by k_{p} (blue) and k_{p}^{−8/3} (red), to see its scalings in the subinertial and inertial ranges, respectively.

FIG. 8. Similar to Fig. 7, but for a case with β = 1. On the right, the solutions for Sp(T) are scaled by k_{p} (blue) and k_{p}^{14/3} (red), to see its scalings in the subinertial and inertial ranges, respectively. Violet indicates the two overlap for k_{p} = 1.

III. PREDICTIONS FROM ALGEBRAICALLY GROWING STRESS

With the detailed information above, we are now in a position to compare with the predictions from a simple time-dependent model. In this section, we compute GW spectra by considering a simple model for the time evolution of the stress. It is assumed to increase algebraically as a power law characterized by a power law index p during the time interval from t = 1 to t_{e}.

A. The model

We model the + and × polarizations of the Fouriertransformed stress, T̃(k, t), as

T̃(k, t) = \{ \begin{matrix}|T̃_{0}(k)| t^{p}, & 1 ≤ t ≤ t_{e}, \\ 0, & t > t_{e}, \end{matrix}  (5)

where |T̃_{0}(k)| represents |T̃(k, t)| at the initial time and is obtained for given energy and helicity spectra of the magnetic field; see Appendix A for details. We note that the authors of Ref. [26] have developed an analytical model for the GW spectrum on the basis of the time evolution of the stress, which they assumed constant during a certain interval – unlike our case. The authors explain the location of certain breaks in their GW spectrum as a consequence of the finite duration over which the stress is constant. This duration is an empirical input parameter. In our model, by contrast, the stress evolves as a power law with an index that is in principle known from MHD theory, although we can get even better agreement with the simulations when we take the actual power-law index that is realized in the simulations.

To obtain the GW spectrum for our model, we first solve Eq. (4) for a case when the source is active during the interval 1 < t < t_{e} and thus obtain h̃(k, t) and ˙ h̃(k, t).

FIG. 9. (a) Sp(ḣ(k, t) at different times. Here, we assume t = t_{e}, E_{M} = c(k/k_{p})^{4}/(1+(k/k_{p})^{17/3}), where k_{p} = 10 and c = 10^{−4} and p = −1/4. The red, blue, and black curves are for t_{e} = 2, 4, and 10 respectively. The two black vertical lines correspond to k_{∗} and 2k_{p}. (b) Sp(ḣ) at times t_{e} = 10, 20, and 30.

The solution for t ≥ t_{e} is given by

h̃(k, t) = ∫_{1}^{t} \frac{sin k(t − t′)}{k} \frac{6T̃(k, t′)}{t′} dt′,  (6)

˙h̃(k, t) = ∫_{1}^{t} cos k(t − t′) \frac{6T̃(k, t′)}{t′} dt′.  (7)

Using Eq. (7) and our model for T̃(k, t), we obtain

˙h̃(k, t) = \frac{−3|T̃_{0}(k)|}{(kt_{0})^{p}}\{e^{i(kt−pπ/2)}[Γ(p, ikt_{e}) − Γ(p, ikt_{0})]
+ e^{−i(kt−pπ/2)}[Γ(p, −ikt_{e}) − Γ(p, −ikt_{0})]}. (8)

In the above expression, t_{0} = 1 represents the initial time. In Fig. 9(a), we show Sp(ḣ) at different times for this model with p = −1/4. The red, blue, and black curves represent Sp(ḣ) at t = 2, 4, and 10, respectively. It is evident from this figure that Sp(ḣ) is almost flat for 1 < ∼ k < ∼ 2k_{p} and declines as ∝ k^{−11/3} for k > 2k_{p}. Sp(ḣ) is proportional to k^{2} for k < k_{H}, where k_{H} represents the wave number corresponding to the Hubble horizon size at t = t_{e}. Further, as time increases, Sp(ḣ) at low wave numbers (k_{H} < k < 1) grows and saturates, as is evident from Fig. 9(b).

To understand the role of the power-law index in the algebraically growing part of the stress on Sp(ḣ), we calculate Sp(ḣ) for different values of p. Those are shown in the right-hand panel of Fig. 10. Here, Sp(ḣ) is rapidly oscillating, so we plot in this figure only its envelope. From this figure, we conclude that Sp(ḣ) can be divided into three regime. We begin discussing first the high wave number regime (k > k_{0}, regime I), where k_{0} represents the wave number corresponding to the Hubble horizon at the initial time. Sp(ḣ) is flat and changes to k^{−11/3} for k > 2k_{p}. For very low wave numbers corresponding to the superhorizon range (k < k_{H}, regime III) Sp(ḣ) is proportional to k^{2}. In the intermediate regime [k_{H} < ∼ k < ∼ (1 − p)/t_{0}, regime II], Sp(ḣ) changes from a

FIG. 10. Sp(ḣ) for different values of p. Here, the left, middle, and right black vertical lines represent the wave numbers corresponding to the horizon size at the final time t, the initial time t_{0} = 1, and 2k_{p}, respectively.

flat spectrum to a k^{2} spectrum as the wave number decreases. Note that, as the wave numbers decrease, the transition from a flat spectrum to a k^{2} spectrum is faster for the case when p = −1/4 than when it is p = 1/3. The wave number at which this transition occurs depends on the value of p and can be understood as follows. In the algebraically growing phase, the typical time scale over which T̃/t decays, is δt_{T} ∼ t/(1 − p) and the typical time scale for sourcing GWs at a given wave number k just after t = 1 is δt_{GW} ∼ 1/k, as can be inferred from the cosine function in Eq. (7). The value of T̃/t does not change much when δt_{GW}/δt_{T} ≤ 1. This implies that, for k > (1 − p)/t_{0}, there will be a finite interval during which T̃/t can be assumed constant. However, for k < (1 − p)/t_{0}, T̃/t always changes. The wave number k ∼ (1 − p)/t_{0} corresponds to the wave number where Sp(ḣ) starts changing from a flat spectrum.

The nature of Sp(ḣ) can also be understood by writing the expression of ˙ h̃(k, t), given in Eq. (8), for different limits depending on the values of kt_{0} and kt_{e}. For t_{e} ≫ t_{0}, which is indeed the case, and p < 1, Eq. (8) reduces to

\frac{˙h̃(k, t)}{6|T̃_{0}(k)|} ≈ \{ \begin{matrix}\frac{sin k(t−t_{0})}{kt_{0}} (I), \\ \frac{Γ[p]}{(kt_{0})^{p}} cos (kt − \frac{pπ}{2} ) − \frac{cos kt}{p} (II), \\ \frac{cos kt}{p} [(\frac{t_{e}}{t_{0}} )^{p} − 1] (III). \end{matrix} (9)

Using this, we calculate the spectrum of ˙ h̃; it is given by

\frac{Sp(ḣ)(k, t)}{36|T̃_{0}(k)|^{2}} ≈ \{ \begin{matrix}[\frac{sin k(t−t_{0})}{t_{0}}]^{2} (I), \\ k^{2}[ − \frac{cos kt}{p} + \frac{Γ[p]}{(kt_{0})^{p}} cos(kt − \frac{pπ}{2} )]^{2} (II), \\ k^{2}\{\frac{cos kt}{p} [(\frac{t_{e}}{t_{0}} )^{p} − 1] \}^{2} (III). \end{matrix}

(10)

From the above expression, we conclude that the break points for the different slopes of Sp(ḣ) is decided by |T̃_{0}(k)|^{2} for k > t_{0}^{−1}. Here, |T̃_{0}(k)|^{2} is flat for k < 2k_{p} and proportional to k^{−11/3} for k > 2k_{p}. For the superhorizon modes, i.e., k < 1/t_{e}, Sp(ḣ) is proportional to k^{2}, and for wave numbers t_{0}^{−1} < k < t_{e}^{−1}, Sp(ḣ) changes from a flat spectrum to k^{2}, as shown as the blue curves in Fig. 10.

In this model, we take the same algebraic evolution with a constant power-law index for all wave numbers. In general, however, the time evolution of T̃(k, t) is different for wave numbers below and above the peak of Sp(T)(k, t). For the case of helical magnetic fields discussed in Fig. 7, the value of Sp(T) for a particular k < 2k_{p} at the initial time, grows as t^{2/3} until the time for which Sp(T) peaks at this particular wave number. After this time, the value of Sp(T) for this particular k starts decreasing as t^{−16/9}. This implies that we would have assumed T̃(k) ∝ t^{1/3} for k < k_{p} and T̃(k) ∝ t^{−8/9} for k > k_{p}. For the nonhelical magnetic field shown in Fig. 8, Sp(T) is always decreasing. For k < 2k_{p} at the initial time, it first decreases as t^{−1/2} and later switches to t^{−7/3}. We study Sp(ḣ) by incorporating this and find that there is no difference in the final Sp(ḣ) compared to the one obtained in our model. This is why we did not consider the aforementioned evolution of the stress to keep the model simple.

IV. COMPARISON OF THE ANALYTICAL MODEL WITH SIMULATION RESULTS

In the above section, we discussed Sp(ḣ) in a model inspired by the fact that the GW spectrum does not change if we change the stress tensor by its modulus for decaying MHD turbulence in the early Universe. In this model, we approximate the stress tensor by |T̃(k, t)| ≡ \sqrt{Sp(T)/4πk^{2}} and its time evolution is parameterized as a power law with index p; see Fig. 11. In this section, we provide a comparison for Sp(ḣ) obtained in this model with the simulation results discussed in Section II B. To compare, we show Sp(ḣ) obtained from the model and different simulations together. In Figs. 12(a) and (b), we plot the spectra for the runs shown in Figs. 1 and 2 and discussed in Section II B. For the dotted-dashed black curve, |T̃_{0}(k)| is obtained by using Eq. (A6) of Appendix A, where we take E_{M} = c(k/k_{p})^{4}/[1 + (k/k_{p})^{17/3}] and the value of the constant c is determined such that the obtained stress spectrum matches with that of the simulation at t = 1. For this case, the stress spectrum evolution is modeled as a single power law and the value of p is 1/3 and −1/4 for the helical and nonhelical cases, respectively. The value of p is decided from the time evolution of the low wave number tail of Sp(T), as discussed in Section II C. From this figure, it is concluded that the spectral nature of Sp(ḣ) matches well with the prediction from the model. However, there is a difference at small wave numbers, especially for the nonhelical case. This is due to the fact that the modelling of Sp(T) by a single power does not provide a better fit to the evolution obtained in NHEL.

A double power law of the form t^{−1/3}/[1+(t−1)^{n}]^{5/7n}, where n regulates the transition, here with n = 10, provides a better fit to Sp(T) for the low frequency tail for NHEL; see Fig. 11. In this figure, we plot T̃(k, t) obtained from the simulation in solid and dotted blue for the wave numbers k = 0.3 and 0.5, respectively and the double power law fit to the blue curves is in dashed red color. The double power law, which fits T̃(k, t) for HEL, is 1/[1 + (t − 0.2)^{n}]^{5/24n}, where n = 20. After considering such a time evolution, the obtained Sp(ḣ) is shown as the dotted red curve in Fig. 12(b). For the dashed green curve, we consider |T̃_{0}(k)| = \sqrt{Sp(T)/4πk^{2}} and Sp(T) is obtained from the simulation at t = 1. These different forms of the time evolution of T̃(k, t) are given in Table III. The spectra in Figs. 12(a) and (b) are plotted at a time when the value of the Sp(ḣ) for each wave number has reached approximately a constant value. The actual time for HEL is t = 175.5 and for NHEL it is t = 99.5. For earlier times, the mean value of Sp(ḣ) is in reasonably good agreement with the values obtained from the model.

We notice that the nature of the GW spectra in HEL and NHEL are different. There is large power in HEL compared to NHEL around the peak of the GW spectrum for the same strength of the initial magnetic field. This is due to the presence of an additional term due to the helicity spectrum in the stress spectrum; see Appendix A.

V. CONCLUSIONS

In this work, we have suggested a simple model to understand the GW spectrum obtained for decaying MHD turbulence in the early Universe. The Fouriertransformed stress is taken to be |T̃(k, t)|, i.e., we ignore changes in the phase, and its time evolution is parameter-

FIG. 11. T̃(k, t) vs. t. (a) The solid and dotted blue curves correspond to the time evolution of T̃(k, t) vs t obtained from the simulation for k = 0.3 and 0.5, respectively, for HEL and the red curve corresponds to a broken power law fit to the blue curves. (b) Same as (a), but for NHEL.

FIG. 12. Sp(ḣ) vs k: (a) The solid blue curve represents Sp(ḣ) corresponding to the run shown in Fig. 1, respectively. The dotted red and dashed green curves represent the spectra obtained in the model for the time evolution of T̃(k, t) given in column 2 of Table III. The dot dashed black curve represents Sp(ḣ) for the time evolution of T̃(k, t) given in column 1 of Table III. The red curves are for the case when |T̃_{0}(k)| is obtained using Eq. (A6) of Appendix A. For the green curves, |T̃_{0}(k)| is taken from the simulation. (b) Same as (a), but for nonhelical case.

TABLE III. Time dependence of T̃(k, t) taken in our analysis

run from theory from simulation
hel  ( \frac{t}{t_{0}} )^{1/3}  \frac{1}{(1+(t−0.2)^{n})^{5/24n}}
nonhel ( \frac{t}{t_{0}} )^{−1/4}  \frac{t^{−1/6}}{(1+(t−1)^{n})^{5/14n}}

ized by a power law. Such a time evolution of the stress is motivated by the simulations for the decaying MHD turbulence at low wave numbers discussed in Section II B. We find that the spectral nature of the GW spectrum is well represented by this simple model. In this work, we also show that the nature of the GW spectra in the helical case are different from those in the nonhelical case. Apart from the polarization of GW, this spectral difference may also be important in distinguishing the helical and nonhelical nature of the primordial magnetic field.

In this work, we have developed a model to understand the low frequency tail of the GW spectrum in the cases where turbulence is initiated suddenly. However, it will now also be interesting to study cases where the magnetic field is generated selfconsistently, such as through the chiral magnetic effect in the early Universe [52, 53]. It would be interesting to see if a model such as the one discussed in this paper can also explain the GW spectra obtained through the chiral magnetic effect. This, we hope to report in a future study.

ACKNOWLEDGMENTS

RS would like to thank Hongzhe Zhou for helping him analyzing the output data with the Mathematica routines of the Pencil Code. This work was supported by the Swedish Research Council (Vetenskapsradet, 2019-04234). Nordita is sponsored by Nordforsk. 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 and Linköping.

Data availability—The source code used for the simulations of this study, the Pencil Code, is freely available from Refs. [54]. Supplemental Material with plots similar to Figs. 3 and 4, but for other wave vectors of the same length, along with the simulation setups and the corresponding data are freely available from Ref. [44].

Appendix A: Sp(T) in terms of magnetic spectrum

Here, we provide the expressions for the stress spectrum in terms of the magnetic spectrum for helical and nonhelical cases; see Refs. [45] and [37] for the derivation. The two-point correlation for the nonhelical magnetic field in Fourier space is given by,

⟨B̃_{i}(k)B̃_{j}^{∗}(k′)⟩ = (2π)^{3}δ(k − k′)(δ_{ij} − k̂_{i}k̂_{j})P_{SM}(k) (A1)

Assuming a Gaussian nature of the magnetic field fluctuations, the stress spectrum is given by

Sp(T) ≡4πk^{2}⟨T̃_{ij}^{TT}(k)T̃_{∗ij}^{TT}(k)⟩ = 4πk^{2} ∫ d^{3}q[P_{SM}(q)
P_{SM}(|k − q|)(1 + γ^{2} + β^{2} + γ^{2}β^{2})].  (A2)

Here, T̃_{ij}^{TT}(k) = −B_{i}B_{j} + \frac{1}{2}δ_{ij}B_{k}B_{k}. In terms of the energy spectrum, E_{M}(k) ≡ 4πk^{2}P_{SM}(k), the above expression reduces to

Sp(T) = \frac{1}{4π} ∫ d^{3}q \frac{k^{2}}{q^{2}|k − q|^{2}} [E_{M}(q)
E_{M}(|k − q|)(1 + γ^{2} + β^{2} + γ^{2}β^{2})].  (A3)

In the above expressions γ = k̂ · q̂ and β = k̂ · k̂ − q. For helical magnetic fields, there is an additional antisymmetric contribution to the two-point correlation and is given by,

⟨B̃_{i}(k)B̃_{j}^{∗}(k′)⟩ = (2π)^{3}δ(k − k′)((δ_{ij} − k̂_{i}k̂_{j})P_{SM}(k)
+ iϵ_{ijm}k̂_{m}P_{AM}(k)).  (A4)

The stress spectrum for this case is given by

Sp(T) =4πk^{2} ∫ d^{3}q[P_{SM}(q)P_{SM}(|k − q|)(1 + γ^{2} + β^{2} + γ^{2}β^{2})
+ 4γβP_{AM}(q)P_{AM}(|k − q|)].  (A5)

In terms of the energy spectrum, E_{M}(k), and the helicity spectrum, H_{M}(k) ≡ 4πk^{2}P_{AM}(k), the above expression reduces to

Sp(T) =\frac{1}{4π} ∫ d^{3}q \frac{k^{2}}{q^{2}|k − q|^{2}} [E_{M}(q)E_{M}(|k − q|)
(1 + γ^{2} + β^{2} + γ^{2}β^{2}) + 4γβH_{M}(q)H_{M}(|k − q|)]. (A6)

[1] M. Kamionkowski,  A. Kosowsky,  and
M.  S.  Turner,  Phys. Rev. D 49, 2837 (1994),

arXiv:astro-ph/9310044 [astro-ph].

[2] C. J. Hogan, Phys. Rev. Lett. 85, 2044 (2000).

[3] A. Kosowsky, A. Mack, and T. Kahniashvili, Phys. Rev. D 66, 024030 (2002), arXiv:astro-ph/0111483 [astro-ph].

[4] A. D. Dolgov, D. Grasso, and A. Nicolis, Phys. Rev. D 66, 103505 (2002), arXiv:astro-ph/0206461 [astro-ph].

[5] C. Caprini and D. G. Figueroa, CQGra 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO].

[6] K. Kajantie, M. Laine, K. Rummukainen, and M. Shaposhnikov, Phys. Rev. Lett. 77, 2887 (1996).

[7] Y. Aoki, G. EndrHodi, Z. Fodor, S. D. Katz, and K. K. Szabó, Nature (London) 443, 675 (2006), arXiv:hep-lat/0611014 [hep-lat].

[8] E. Witten, Phys. Rev. D 30, 272 (1984).

[9] C. J. Hogan, Mon. Not. Roy. Astron. Soc. 218, 629 (1986).

[10] A. Mazumdar and G. White, Reports on Progress in Physics 82, 076901 (2019).

[11] M. Hindmarsh, M. Lben, J. Lumma, and M. Pauly, SciPost Phys. Lect. Notes , 24 (2021).

[12] M. S. Turner and L. M. Widrow, Phys. Rev. D 37, 2743 (1988).

[13] T. Vachaspati, PhLB 265, 258 (1991).

[14] B. Ratra, Astrophys. J. Lett. 391, L1 (1992).

[15] R. Durrer and A. Neronov, Astron. Astrophys. Rev. 21, 62 (2013), arXiv:1303.7121 [astro-ph.CO].

[16] K. Subramanian, RPPh 79, 076901 (2016).

[17] K. Subramanian, Galaxies 7, 47 (2019), arXiv:1903.03744 [astro-ph.CO].

[18] T. Vachaspati, Rep. Prog. Phys. (2021).

[19] A.  Brandenburg,  Y.  He,  and
R.  Sharma,  Astrophys. J. 922, 192 (2021),

arXiv:2107.12333 [astro-ph.CO].

[20] Z. Arzoumanian et al. (NANOGrav), Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE].

[21] B. Goncharov et al., Astrophys. J. Lett. 917, L19 (2021), arXiv:2107.12112 [astro-ph.HE].

[22] S. Chen et al., Mon. Not. Roy. Astron. Soc. 508, 4970 (2021), arXiv:2110.13184 [astro-ph.HE].

[23] J. Antoniadis et al., Mon. Not. Roy. Astron. Soc. 510, 4873 (2022), arXiv:2201.03980 [astro-ph.HE].

[24] A. Neronov, A. R. Pol, C. Caprini, and D. Semikoz, Phys. Rev. D 103, L041302 (2021), arXiv:2009.14174 [astro-ph.CO].

[25] R. Sharma, Phys. Rev. D 105, L041302 (2022), arXiv:2102.09358 [astro-ph.CO].

[26] A. Roper Pol, C. Caprini, A. Neronov, and D. Semikoz, Phys. Rev. D 105, 123502 (2022), arXiv:2201.05630 [astro-ph.CO].

[27] R. W. Hellings and G. S. Downs, Astrophys. J. Lett. 265, L39 (1983).

[28] A. Roper Pol,  S. Mandal,  A. Bran-
denburg,  T.  Kahniashvili,  and
A.  Kosowsky,  Phys. Rev. D 102, 083512 (2020),

arXiv:1903.08585 [astro-ph.CO].

[29] P. Auclair, C. Caprini, D. Cutting, M. Hindmarsh, K. Rummukainen, D. A. Steer, and D. J. Weir, (2022), arXiv:2205.02588 [astro-ph.CO].

[30] A. Brandenburg, T. Kahniashvili, and A. G. Tevzadze, PhRvL 114, 075001 (2015), arXiv:1404.2238.

[31] J. Zrake, Astrophys. J. Lett. 794, L26 (2014), arXiv:1407.5626 [astro-ph.HE].

[32] A. Kosowsky, A. Mack, and T. Kahniashvili, PhRvD 66, 024030 (2002).

[33] G. Gogoberidze,  T. Kahniashvili,  and
A.  Kosowsky,  Phys. Rev. D 76, 083002 (2007),

arXiv:0705.1733 [astro-ph].

[34] T. Kahniashvili, A. Kosowsky, G. Gogoberidze, and Y. Maravin, PhRvD 78, 043003 (2008).

[35] C. Caprini, R. Durrer, and G. Servant, J. Cosmol. Astropart. Phys. 2009, 024 (2009), arXiv:0909.0622 [astro-ph.CO].

[36] P. Niksa, M. Schlederer, and G. Sigl, Class. Quant. Grav. 35, 144001 (2018), arXiv:1803.02271 [astro-ph.CO].

[37] R. Sharma, K. Subramanian, and T. R. Seshadri, Phys. Rev. D 101, 103526 (2020), arXiv:1912.12089 [astro-ph.CO].

[38] A. Roper Pol, S. Mandal, A. Brandenburg, and T. Kahniashvili, J. Cosmol. Astropart. Phys. 2022, 019 (2022), arXiv:2107.05356 [gr-qc].

[39] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal, and A. R. Pol, PhRvR 3, 013193 (2021), arXiv:2011.05556 [astro-ph.CO].

[40] J. Dahl, M. Hindmarsh, K. Rummukainen, and D. Weir, arXiv e-prints , arXiv:2112.12013 (2021), arXiv:2112.12013 [gr-qc].

[41] A. Roper Pol, A. Brandenburg, T. Kahniashvili, A. Kosowsky, and S. Mandal, GApFD 114, 130 (2020).

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

arXiv:astro-ph/9602031 [astro-ph].

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

[44] R. Sharma and A. Brandenburg, Supplemental Material and Datasets for Low frequency tail of gravitational wave spectra from hydromagnetic turbulence, doi:10.5281/zenodo.6596525 (v2022.05.31); see also http://norlx65.nordita.org/~brandenb/projects/LowFreqTail/ for easier access 10.5281/zenodo.6596525.

[45] C. Caprini, R. Durrer, and T. Kahniashvili, Phys. Rev. D 69, 063006 (2004), arXiv:astro-ph/0304556 [astro-ph].

[46] T. Hatori, JPSJ 53, 2539 (1984).

[47] D. Biskamp and W.-C. Müller, Phys. Rev. Lett. 83, 2195 (1999), arXiv:physics/9903028 [physics.flu-dyn].

[48] A. Brandenburg and T. Kahniashvili, PhRvL 118, 055102 (2017), arXiv:1607.01360 [physics.flu-dyn].

[49] A. A. Schekochihin, arXiv e-prints , arXiv:2010.00699 (2020), arXiv:2010.00699 [physics.plasm-ph].

[50] D. N. Hosking and A. A. Schekochihin, Phys. Rev. X 11, 041005 (2021).

[51] H. Zhou, R. Sharma, and A. Brandenburg, arXiv e-prints , arXiv:2206.07513 (2022), arXiv:2206.07513 [physics.plasm-ph].

[52] I. Rogachevskii, O. Ruchayskiy, A. Boyarsky, J. Fröhlich, N. Kleeorin, A. Brandenburg, and J. Schober, Astrophys. J. 846, 153 (2017), arXiv:1705.00378 [physics.plasm-ph].

[53] J. Schober, I. Rogachevskii, A. Brandenburg, A. Boyarsky, J. Fröhlich, O. Ruchayskiy, and N. Kleeorin, Astrophys. J. 858, 124 (2018).

[54] Pencil Code Collaboration, A. Brandenburg, A. Johansen, P. Bourdin, W. Dobler, W. Lyra, M. Rheinhardt, S. Bingert, N. Haugen, A. Mee, F. Gent, N. Babkovskaia, C.-C. Yang, T. Heinemann, B. Dintrans, D. Mitra, S. Candelaresi, J. Warnecke, P. Käpylä, A. Schreiber, P. Chatterjee, M. Käpylä, X.-Y. Li, J. Krüger, J. Aarnes, G. Sarson, J. Oishi, J. Schober, R. Plasson, C. Sandin, E. Karchniwy, L. Rodrigues, A. Hubbard, G. Guerrero, A. Snodin, I. Losada, J. Pekkilä, and C. Qian, JOSS 6, 2807 (2021).