%\documentclass[twocolumn,twocolappendix]{aastex631} \documentclass[twocolumn,aps,prd,nofootinbib]{revtex4-1} \usepackage{graphicx} \usepackage{hyperref} \usepackage{xcolor} \def\kf{k_{\rm f}} \graphicspath{{./fig/}{./png/}} \usepackage{amsmath} \usepackage{macros} \begin{document} %\author[0000-0002-2549-6861]{Ramkishor Sharma} %\author[0000-0002-7304-021X]{Axel Brandenburg} %\email{ramkishor.sharma@su.se} %\email{brandenb@nordita.org} \title{Low frequency tail of gravitational wave spectra from hydromagnetic turbulence} \author{Ramkishor Sharma$^{1,2}$} \email{ramkishor.sharma@su.se} \author{Axel Brandenburg$^{1,2,3,4}$} \email{brandenb@nordita.org} \affiliation{$^{1}$Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfv\'ens v\"ag 12, 10691 Stockholm, Sweden} \affiliation{$^{2}$The Oskar Klein Centre, Department of Astronomy, Stockholm University, 10691 Stockholm, Sweden} \affiliation{$^{3}$McWilliams Center for Cosmology \& Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA} \affiliation{$^{4}$School of Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili Avenue, 0194 Tbilisi, Georgia} \begin{abstract} 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. \end{abstract} \keywords{gravitational waves---early Universe---turbulence---magnetic fields---MHD} \maketitle \section{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 \citep{1994PhRvD..49.2837K,PhysRevLett.85.2044,2002PhRvD..66b4030K,Dolgov+02,2018CQGra..35p3001C}. The frequency spectrum of these waves is related to the spectrum of the underlying turbulence. Such turbulence could be induced during the various \blue{epochs}\footnote{\blue{The electroweak and QCD epochs are accompanied by crossovers in the standard model \cite{PhysRevLett.77.2887,2006Natur.443..675A}. However, many extensions of the standard model can lead to a first order phase transition.}} in the early Universe \citep{PhysRevD.30.272, PhysRevD.30.272,10.1093/mnras/218.4.629, Mazumdar_2019,10.21468/SciPostPhysLectNotes.24} or the possible presence of primordial magnetic fields \citep{1988PhRvD..37.2743T, tanmay1991, 1992ApJ...391L...1R,durrer2013, subramanian2016, kandu2019, tanmay2021, Bran+He+Shar21}. 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 \citep{NANOGrav2020, Goncharov_2021, Chen2021, Antoniadis2022} 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 \citep{2021PhRvD.103L1302N, Sharma21, RoperPol+22}. However, the presence of a quadrupolar spatial correlation \citep{Hellings&Downs}, 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 \citep{RoperPol+20}, 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 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. \blue{To understand the origin of the low frequency tail,} the authors of Ref.~\cite{RoperPol+22} 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. \blue{Furthermore, in the numerical solutions the stress was not actually constant for any duration of time.} In another recent approach, the authors of Ref.~\cite{Auclair:2022jod} 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 \cite{2015PhRvL.114g5001B,2014ApJ...794L..26Z}. This will be crucial to the approach discussed in the present paper. To address the problem of a limited computational domain, it is important \blue{ 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 \blue{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.~\cite{RoperPol+22}, 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 \Sec{TheModel}, 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 \Sec{simplemodel}, we discuss the model to explain the low frequency tail of the GW spectrum. Further, in \Sec{comparison}, we compare the GW spectrum obtained from our numerical simulations and our model. We conclude in \Sec{conclusion}. \section{Nonhelical and helical cascades} \label{TheModel} 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 \citep{Dolgov+02,kosowsky2002,Gogo+07,tina2008,Chiara2009,sigl2018,Sharma+20} and numerically \citep{RoperPol+20,RoperPol+21,Kahniashvili+21,Jani2021,RoperPol+22,Auclair:2022jod}. 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. \subsection{GWs from MHD turbulence} \label{GWfromCME} We follow here the formalism of Ref.~\cite{RoperPol+20b, RoperPol+20}, where conformal time is normalized to unity at the initial time. One could associate this with the electroweak phase transition, for example. The velocity $\uu$ is normalized to the speed of light. The magnetic field $\BB=\nab\times\AAA$ is written in terms of the magnetic vector potential $\AAA$, and the current density is written as $\JJ=\nab\times\BB$. Following Ref.~\cite{BEO96}, the energy density $\rho$ includes the restmass density, so its evolution equation obeys a continuity equation that also includes magnetic energy terms. As in \cite{RoperPol+20b}, $\rho$ 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, $\tilde{h}_+$ and $\tilde{h}_\times$, which are driven by the corresponding projections of the stress, which, in turn, is composed of kinetic and magnetic contributions, \EQ \TTT_{ij} =\frac{4}{3}\gamma_{\rm Lor}^2\rho u_i u_j-B_i B_j+..., \EN where $\gamma_{\rm Lor}=(1-\uu^2)^{-1/2}$ is the Lorentz factor, and the ellipsis denotes terms proportional to $\delta_{ij}$, which do not contribute to the projected source $\tilde{T}_{+/\times}$. Assuming the Universe to be conformally flat, its expansion can be scaled out by working with conformal time $t$ and comoving variables \citep{BEO96}. 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 \citep{RoperPol+20b}. The full set of equations is therefore \begin{eqnarray} &&\frac{\partial\BB}{\partial t}= \nab\times(\uu\times\BB-\eta\nab\times\BB),\label{dAdt}\\ &&{\DD\uu\over\DD t}= {1\over\rho}\nab\cdot\left(2\rho\nu\SSSS\right)-{1\over4}\nab\ln\rho +{\uu\over3}\left(\nab\cdot\uu+\uu\cdot\nab\ln\rho\right) \nonumber \\ &&\qquad\quad-{\uu\over\rho} \left[\uu\cdot(\JJ\times\BB)+\eta \JJ^2\right] +{3\over4\rho}\JJ\times\BB, \label{dudt} \\ &&{\partial\ln\rho\over\partial t} =-\frac{4}{3}\left(\nab\cdot\uu+\uu\cdot\nab\ln\rho\right) +{1\over\rho}\left[\uu\cdot(\JJ\times\BB)+\eta \JJ^2\right]\!, \nonumber \\ &&\frac{\partial^2}{\partial t^2} \tilde{h}_{+/\times} (\kk, t) +k^2\tilde{h}_{+/\times} (\kk, t) = {6\over t} \tilde{T}_{+/\times}(\kk,t), \label{GW4} \end{eqnarray} where $\DD/\DD t\equiv\partial/\partial t+\uu\cdot\nab$ is the advective derivative, $\eta$ is the magnetic diffusivity, $\nu$ is the kinematic viscosity, ${\sf S}_{ij}=\half(u_{i,j}+u_{j,i})-\onethird\delta_{ij}\nab\cdot\uu$ are the components of the rate-of-strain tensor $\SSSS$ 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 $\BB$ consists of a weak Gaussian-distributed seed magnetic field, $\uu=0$, $\rho=1$. We work with spectra that are defined as integrals over concentric shells in wave number space $\kk$ with $k=|\kk|$. They are normalized such that their integrals over $k$ give the mean square of the corresponding quantity, i.e., $\int\Sp(\BB)\,\dd k=\bra{\BB^2}$, where $\Sp(\BB)=\Sp(B_x)+\Sp(B_y)+\Sp(B_z)$. Similarly, $\Sp(\hh)=\Sp(h_+)+\Sp(h_\times)$ is defined as the sum over the two polarization modes. Of particular interest will also be the stress spectrum $\Sp(\TT)$, which is defined analogously through $\Sp(\TT)=\Sp(T_+)+\Sp(T_\times)$. To study the evolution of the stress at selected Fourier modes, we compute $|\tildeT(k,t)|\equiv\sqrt{\Sp(\TT)/4\pi k^2}$, which scales the same way as $|\tilde{T}_+(k,t)|$ and $|\tilde{T}_\times(k,t)|$. \subsection{Evolution of the stress and strain spectra} \label{EvolutionOfStress} To put our results into perspective and compare with earlier work, we study cases of suddenly initiated turbulence. We perform simulations similar to those of Ref.~\cite{RoperPol+20} by using as initial condition for the magnetic field a random Gaussian-distributed magnetic field with a $k^4$ spectrum for $kk_{\rm p}$. For details of such a magnetic field, see Ref.~\cite{Bran+17}. As initial condition for the GW field, we assume that $h$ and $\dot{h}$ 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 \Tab{table1}. In this table, the first column represents the name of the runs (HEL and NHEL), $\EEEM^{i}$ is the initial value of the magnetic energy density compared to the background energy density, $\kp$ 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, $\EEGW^{\rm sat}$ is the value of the GW energy density after saturation compared to the background energy density, and $\OmGW^{\rm sat}$ is the density parameter of GWs, representing the ratio of the GW energy density compared to the critical energy density at present. $\OmGW^{\rm sat}$ has been calculated considering the production of GWs around the electroweak phase transition. %Fig.1 \begin{figure*}\begin{center} %sayan/GW/LowFreq1024sig1_k01_kf10b_rep2 %/cfs/klemming/projects/snic/snic2020-4-12/brandenb/sayan/GW/LowFreq1024sig1_k01_kf10b_rep2 \includegraphics[width=\textwidth]{pstress_etc_LowFreq_sig1} \end{center}\caption{ Spectra of the magnetic field, the TT-projected stress, the strain derivative, and the strain for suddenly initiated turbulence with magnetic helicity. }\label{pstress_etc_LowFreq_sig1}\end{figure*} %Fig.2 \begin{figure*}\begin{center} \includegraphics[width=\textwidth]{pstress_etc_LowFreq_sig0} %sayan/GW/LowFreq1024sig0_k01_kf10b2_rep2 %/cfs/klemming/projects/snic/snic2020-4-12/brandenb/sayan/GW/LowFreq1024sig0_k01_kf10b2_rep2 \end{center}\caption{ Same as \Fig{pstress_etc_LowFreq_sig1}, but for the nonhelical case. }\label{pstress_etc_LowFreq_sig0}\end{figure*} \blue{We consider two runs where the initial magnetic field is either fully helical ($\runh$) or nonhelical ($\runnh$).} In \Figs{pstress_etc_LowFreq_sig1}{pstress_etc_LowFreq_sig0}, we show \blue{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 \Figsp{pstress_etc_LowFreq_sig1}{a}{b} for the $\runh$ run. \blue{By contrast, in $\runnh$, the stress spectrum always decreases at small $k$.} We also see in \Figp{pstress_etc_LowFreq_sig1}{c} that the GW energy spectrum has a maximum at $k\approx20$, which is not present in $\runnh$; cf.\ \Figp{pstress_etc_LowFreq_sig0}{c}. Their spectra fall off toward smaller $k$ proportional to $k$ and $k^{1.5}$ for $\runh$ and $\runnh$, respectively. \begin{table}\caption{ Summary of simulation parameters }\begin{center} \begin{tabular}{ccccccc} \hline Run & $\EEM^i$ &$\kp$ & $\EEGW^{\rm sat}$ & $\OmGW^{\rm sat}$\\ \hline \runh & $5.4\times10^{-3}$ &$10$& $3.7\times10^{-7}$ & $5.9\times10^{-12}$\\ %LowFreq1024sig1_k01_kf10b_rep2 \hline \runnh & $5.5 \times 10^{-3}$ &$10$& $3.5\times10^{-7}$ & $5.6\times10^{-12}$\\%LowFreq1024sig0_k01_kf10b2_rep2 \hline \label{table1}\end{tabular} \end{center} \end{table} %Fig.3 \begin{figure*}\begin{center} % sayan/GW/LowFreq1024sig1_k01_kf10b_rep2 \includegraphics[width=\textwidth]{rslice_stress_plot_LowFreq} \end{center}\caption{ Modulus and phase of $\tilde{T}(k,t)$ and $\dot{\tilde{h}}(k,t)$ for the helical case for $\kk=(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. }\label{rslice_stress_plot_LowFreq}\end{figure*} %Fig.4 \begin{figure*}\begin{center} \includegraphics[width=\textwidth]{rslice_stress_plot_LowFreq_nohel} \end{center}\caption{ Same as \Fig{rslice_stress_plot_LowFreq}, but for the nonhelical case. }\label{rslice_stress_plot_LowFreq_nohel}\end{figure*} %Fig.5 \begin{figure*}\begin{center} %\includegraphics[width=\textwidth]{pcascade} \includegraphics[scale=0.7]{helical_with_and_wo_phase.eps} \includegraphics[scale=0.7]{nonhelical_with_and_wo_phase.eps} \end{center}\caption{$\Sp(\dot{\hh})(k,t)$ vs $k$: (a) The solid black and blue curves represent $\Sp(\dot{\hh})(k,t)$ at times $t=1.5$ and $t=37$ for run hel. The dashed red and gray curves show $\Sp(\dot{\hh})(k,t)$ for the case when the stress spectrum has been replaced by its modulus in the GW evolution equation. \blue{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.} \label{with_and_wo_phase} \end{figure*} In \Figs{rslice_stress_plot_LowFreq}{rslice_stress_plot_LowFreq_nohel}, we \blue{show the time evolution of the modulus and phase of stress and strain derivative for $\runh$ and $\runnh$ for wave vectors $\kk=(k,0,0)$ and five different values of $k$, which are all below $k_{\rm 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. } \blue{ Broadly speaking, the time evolution of the modulus of the stress $|\tildeT|$ 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 $\runh$ and decreasing for $\runnh$, as was see in \Figs{pstress_etc_LowFreq_sig1}{pstress_etc_LowFreq_sig0}. This is an important observation and may be due to the fact that in \Figs{rslice_stress_plot_LowFreq}{rslice_stress_plot_LowFreq_nohel}, we have shown stress and strain derivatives only for particular values of $k$. } \blue{From \Figs{rslice_stress_plot_LowFreq}{rslice_stress_plot_LowFreq_nohel},} it is evident that $\arg(\tildeT)$ remains constant for some time and starts evolving more rapidly after that. It is also interesting to note that the amplitude of $|\dot{\tildeh}|$ increases up to the time until which $\arg(\dot{\tildeh})$ is roughly constant. After this time, $|\dot{\tildeh}|$ enters an oscillatory regime and its amplitude does not change much. \blue{Other wave vectors of the same length show a similar behavior of $\arg(\tildeT)$ and $\arg(\dot{\tildeh})$, as is shown in the figures of the Supplemental Material provided along with the data in Ref.~\cite{DATA}. } This conclusion applies for both runs shown in \Figs{rslice_stress_plot_LowFreq}{rslice_stress_plot_LowFreq_nohel} and it leads us to develop a simple model to understand the GW spectrum in these cases. In this model we replace $\tildeT$ by its wave vector-averaged magnitude, $|\tildeT|$, as discussed in \Sec{simplemodel}. 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 $\runh$ and $\runnh$, where we replace $\tildeT(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 $\runh$ and $\runnh$ runs; see \Fig{with_and_wo_phase}. 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 $\tildeT(\kk,t)$ has been replaced by its modulus. It is evident from the figure that there is hardly any difference in the actual $\Sp(\dot{\hh})$ 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 $\runh$ and $\runnh$ is the more pronounced peak in the spectral GW energy in the former. As we show in \App{appendixa}, 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{stress_spectrum_helvsnonhel} and the details are explained next. %Fig.6 \begin{figure} \begin{center} \includegraphics[scale=0.6]{helical_vs_nonhelical.eps} \end{center}\caption{ In this figure, magnetic field energy spectrum, $\EM(k)$ (Dashed curves) and $\Sp(\TT)$ (Solid curves) for the helical and non helical case. The blue and red curves are for nonhelical and helical case respectively. }\label{stress_spectrum_helvsnonhel} \end{figure} %Fig.7 \begin{figure*}\begin{center} \includegraphics[scale=0.7]{helical1.eps} \includegraphics[scale=0.7]{helical2.eps} \end{center}\caption{ Left: Solutions for $\Sp(\TT)$ (red) for different $\EM(k)$ (blue) for three values of $k_{\rm p}$. Right: solutions for $\Sp(\TT)$ scaled by $k_{\rm p}$ (blue) and $k_{\rm p}^{-8/3}$ (red), to see its scalings in the subinertial and inertial ranges, respectively. }\label{pcascade}\end{figure*} %Fig.8 \begin{figure*}\begin{center} \includegraphics[scale=0.7]{nonhelical1.eps} \includegraphics[scale=0.7]{nonhelical2.eps} \end{center}\caption{ Similar to \Fig{pcascade}, but for a case with $\beta=1$. On the right, the solutions for $\Sp(\TT)$ are scaled by $k_{\rm p}$ (blue) and $k_{\rm p}^{14/3}$ (red), to see its scalings in the subinertial and inertial ranges, respectively. Violet indicates the two overlap for $k_{\rm p}=1$. }\label{pcascade_beta2}\end{figure*} \subsection{Overall behavior of the stress} \label{stress_evolution} At the most minimalistic level, we can say that the magnetic field shows an approximately self-similar evolution at late times, where for $\runh$, the peak value of $\EM(k,t)$ is unchanged, but the position of the peak $k_{\rm p}$ goes to progressively smaller values as $k_{\rm p}\sim t^{-2/3}$. To understand the consequences for the evolution of the stress, let us now consider an idealized model, where $\EM(k)\equiv\Sp(\BB)/2$ has a $k^4$ subinertial range, where $kk_{\rm p}$. The spectrum of the transverse traceless part of the stress, $\Sp(\TT)/2$, can be computed analytically using the expressions given in \App{appendixa} (see Refs.~\cite{CDK04,Sharma+20} for details) and is shown in \Figs{pcascade}{pcascade_beta2} for the helical and nonhelical cases, respectively. In these figures, we take three instances where the magnetic peaks are at wave numbers $\kp=1$, 0.3, and $0.1$. For HEL, the position of the peak of $\EM(k)$ is unchanged. We see that, in agreement with earlier work \citep{RoperPol+20}, the positions of the peak of $\Sp(\TT)$ are always at $2k_{\rm p}$. However, even though the peak values of $\EM(k)$ are unchanged, except for the factor of two, those of $\Sp(\TT)$ are not and decay. Nevertheless, at small $k$, $\Sp(\TT)$ still increases proportional to $k_{\rm p}^{-1}$. If $k_{\rm p}\propto t^{-2/3}$, as expected for helical turbulence \citep{Hat84,BM99,BK17}, we find that $\Sp(\TT)\propto t^{2/3}$ for small $k$. For NHEL, as shown in Ref.~\cite{BK17}, the peak of the spectrum decreases with decreasing values of $k_{\rm p}$ proportional to $k_{\rm p}^\beta$, where $\beta$ is an exponent that can be between one and four. In \Fig{pcascade_beta2}, we present the case with $\beta=1$ and find that now $\Sp(\TT)(k)\propto k_{\rm p}$ for small $k$ and $\propto k_{\rm p}^{14/3}$ for large $k$. If $\kp\propto t^{-1/2}$, as expected for the nonhelical case for $\beta=1$, $\Sp(\TT)\propto t^{-1/2}$ for small $k$. \blue{We have summarized the behavior of $\Sp(\TT)$ with time in \Tab{tildeT_vs_t} for helical and nonhelical cases.} \begin{table}\caption{ Time evolution of $\Sp(\TT)$ and $|\tilde{T}|$ from theory }\begin{center} \begin{tabular}{c|c|c} &\multicolumn{1}{c}{helical}& \multicolumn{1}{c}{nonhelical}\\ \hline $\Sp(\TT)$ vs $k_{\rm p}$ at small $k$& $\propto 1/k_{\rm p}$& $\propto k_{\rm p}$ \\ \hline $k_{\rm p}$ vs $t$ &$k_{\rm p} \propto t^{-2/3}$&$k_{\rm p} \propto t^{-1/2}$\\ \hline $\Sp(\TT)$ vs $t$ & $\Sp(\TT) \propto t^{2/3}$& $\Sp(\TT)\propto t^{-1/2}$ \\ \hline $|\tilde{T}|$ vs $t$ & $|\tilde{T}| \propto t^{1/3}$& $|\tilde{T}|\propto t^{-1/4}$ \\ \hline \label{tildeT_vs_t}\end{tabular} \end{center} \end{table} Recently, it has been found that the Hosking integral \cite{Sch20}, a Saffman-like helicity integral, is well conserved in nonhelical magnetically dominated decaying turbulence \citep{hosking20,ZBS22}, which implies $\beta=1.5$. For the general case, we write $\Sp(\TT)\propto \kp^{2\beta-1}$ (for $k<\kp$), which implies $\Sp(\TT)\propto t^{-8/9}$ for $\beta=1.5$ and $\kp\propto t^{-4/9}$. It is also interesting to note that, for a given spectrum of the magnetic field, $\Sp(\TT)$ is different for HEL and NHEL; see the blue and red curves in \Fig{stress_spectrum_helvsnonhel}), respectively. For the helical case, $\Sp(\TT)$ has smaller values compared to the nonhelical case at wave numbers below $\kp$. 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. \section{Predictions from algebraically growing stress}\label{simplemodel} 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 $\power$ during the time interval from $t=1$ to $\te$. \subsection{The model} We model the $+$ and $\times$ polarizations of the Fourier-transformed stress, $\tildeT(k,t)$, as \begin{equation}\label{stresscase2} \tildeT(k,t)= \left\{\begin{array}{ll} |\tildeT_0(k)| \, t^{\power} , \quad & 1\le t \le \te, \\ 0, & t > \te, \end{array}\right. \end{equation} where $|\tildeT_0(k)|$ represents $|\tildeT(k,t)|$ at the initial time and is obtained for given energy and helicity spectra of the magnetic field; see \App{appendixa} for details. We note that the authors of Ref.~\cite{RoperPol+22} 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. %Fig.9 \begin{figure*} \includegraphics[scale=0.7]{model.eps} \includegraphics[scale=0.7]{model_fig2.eps} \caption{ (a) $\Sp(\dot{\hh}(k,t)$ at different times. Here, we assume $t=\te$, $\EM=c (k/\kp)^4/(1+(k/\kp)^{17/3})$, where $\kp=10$ and $c=10^{-4}$ and $\power=-1/4$. The red, blue, and black curves are for $\te=2$, $4$, and $10$ respectively. The two black vertical lines correspond to $k_{*}$ and $2 \kp$. (b) $\Sp(\dot{\hh})$ at times $\te=10$, $20$, and $30$. }\label{results_from_model} \end{figure*} To obtain the GW spectrum for our model, we first solve \Eq{GW4} for a case when the source is active during the interval $1 2\kp$. $\Sp(\dot{\hh})$ is proportional to $k^2$ for $k<\kH$, where $\kH$ represents the wave number corresponding to the Hubble horizon size at $t=\te$. Further, as time increases, $\Sp(\dot{\hh})$ at low wave numbers ($\kH < k < 1$) grows and saturates, as is evident from \Fig{results_from_model}(b). To understand the role of the power-law index in the algebraically growing part of the stress on $\Sp(\dot{\hh})$, we calculate $\Sp(\dot{\hh})$ for different values of $\power$. Those are shown in the right-hand panel of \Fig{sphdot_for_different_gamma}. Here, $\Sp(\dot{\hh})$ is rapidly oscillating, so we plot in this figure only its envelope. From this figure, we conclude that $\Sp(\dot{\hh})$ 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(\dot{\hh})$ is flat and changes to $k^{-11/3}$ for $k>2\kp$. For very low wave numbers corresponding to the superhorizon range ($k<\kH$, regime~III) $\Sp(\dot{\hh})$ is proportional to $k^2$. In the intermediate regime [$\kH\lesssim k\lesssim (1-\power)/t_0$, regime~II], $\Sp(\dot{\hh})$ changes from a 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 $\power=-1/4$ than when it is $\power=1/3$. The wave number at which this transition occurs depends on the value of $\power$ and can be understood as follows. In the algebraically growing phase, the typical time scale over which $\tildeT/t$ decays, is $\delta t_T\sim t/(1-\power)$ and the typical time scale for sourcing GWs at a given wave number $k$ just after $t=1$ is $\delta t_{\rm GW}\sim 1/k$, as can be inferred from the cosine function in \Eq{hdot}. The value of $\tildeT/t$ does not change much when $\delta t_{\rm GW}/\delta t_T\le 1$. This implies that, for $k>(1-\power)/t_0$, there will be a finite interval during which $\tildeT/t$ can be assumed constant. However, for $k<(1-\power)/t_0$, $\tildeT/t$ always changes. The wave number $k\sim (1-\power)/t_0$ corresponds to the wave number where $\Sp(\dot{\hh})$ starts changing from a flat spectrum. \begin{figure} \begin{center} \includegraphics[scale=0.8]{sp_hdot_for_different_p.eps} \caption{$\Sp(\dot{\hh})$ for different values of $\power$. 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 $2\kp$, respectively.} \label{sphdot_for_different_gamma} \end{center} \end{figure} The nature of $\Sp(\dot{\hh})$ can also be understood by writing the expression of $\dot{\tildeh}(k,t)$, given in \Eq{full_solution}, for different limits depending on the values of $kt_0$ and $k\te$. For $\te\gg t_0$, which is indeed the case, and $\power<1$, \Eq{full_solution} reduces to \begin{equation} \frac{\dot{\tildeh}(k,t)}{6 |\tildeT_0(k)|}\approx \left\{\begin{array}{lr} \frac{\sin{k(t-t_0)}}{kt_0} \quad \text{(I)},\\[5pt] %k>t_0^{-1} \frac{\Gamma[\power]}{(kt_0)^\power}\cos\left({kt-\frac{\power\pi}{2}}\right)-\frac{\cos{kt}}{\power} \quad \text{(II)},\\[5pt] %t_e^{-1}t_0^{-1}$. Here, $|\tildeT_0(k)|^2$ is flat for $k<2 \kp$ and proportional to $k^{-11/3}$ for $k>2 \kp$. For the superhorizon modes, i.e., $k<1/\te$, $\Sp(\dot{\hh})$ is proportional to $k^2$, and for wave numbers $t_0^{-1}\kp$. For the nonhelical magnetic field shown in \Fig{pcascade_beta2}, $\Sp(\TT)$ is always decreasing. For $k<2 \kp$ at the initial time, it first decreases as $t^{-1/2}$ and later switches to $t^{-7/3}$. We study $\Sp(\dot{\hh})$ by incorporating this and find that there is no difference in the final $\Sp(\dot{\hh})$ 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. \section{Comparison of the analytical model with simulation results} \label{comparison} In the above section, we discussed $\Sp(\dot{\hh})$ 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 $|\tildeT(k,t)|\equiv\sqrt{\Sp(\TT)/4\pi k^2}$ and its time evolution is parameterized as a power law with index $\power$; see \Fig{t_evolution}. In this section, we provide a comparison for $\Sp(\dot{\hh})$ obtained in this model with the simulation results discussed in \Sec{EvolutionOfStress}. To compare, we show $\Sp(\dot{\hh})$ obtained from the model and different simulations together. In \Figsp{helical_comparison}{a}{b}, we plot the spectra for the runs shown in \Figs{pstress_etc_LowFreq_sig1}{pstress_etc_LowFreq_sig0} and discussed in \Sec{EvolutionOfStress}. For the dotted-dashed black curve, $|\tildeT_0(k)|$ is obtained by using \Eq{tijh} of \App{appendixa}, where we take $\EM=c (k/\kp)^4/[1+(k/\kp)^{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(\TT)$, as discussed in \Sec{stress_evolution}. From this figure, it is concluded that the spectral nature of $\Sp(\dot{\hh})$ 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(\TT)$ by a single power does not provide a better fit to the evolution obtained in $\runnh$. 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(\TT)$ for the low frequency tail for NHEL; see \Fig{t_evolution}. In this figure, we plot $\tildeT(k,t)$ obtained from the simulation \blue{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 $\tildeT(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(\dot{\hh})$ is shown as the dotted red curve in \Fig{helical_comparison}(b). For the dashed green curve, we consider $|\tildeT_0(k)|=\sqrt{\Sp(\TT)/4\pi k^2}$ and $\Sp(\TT)$ is obtained from the simulation at $t=1$. These different forms of the time evolution of $\tildeT(k,t)$ are given in \Tab{table2}. \blue{The spectra in \Figsp{helical_comparison}{a}{b} are plotted at a time when the value of the $\Sp(\dot{\hh})$ for each wave number has reached approximately a constant value. The actual time for $\runh$ is $t=175.5$ and for $\runnh$ it is $t=99.5$. For earlier times, the mean value of $\Sp(\dot{\hh})$ is in reasonably good agreement with the values obtained from the model.} We notice that the nature of the GW spectra in $\runh$ and $\runnh$ are different. There is large power in $\runh$ compared to $\runnh$ 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 \App{appendixa}. \begin{table}\caption{ Time dependence of $\tildeT(k,t)$ taken in our analysis }\begin{center} \begin{tabular}{c|c|c} run & from theory & from simulation\\ \hline hel& $\left(\frac{t}{t_0}\right)^{1/3}$& $\frac{1}{(1+(t-0.2)^n)^{5/24n}}$ \\ \hline nonhel& $\left(\frac{t}{t_0}\right)^{-1/4}$& $\frac{t^{-1/6}}{(1+(t-1)^n)^{5/14n}}$ \\ \hline \label{table2}\end{tabular} \end{center} \end{table} \section{Conclusions}\label{conclusion} In this work, we have suggested a simple model to understand the GW spectrum obtained for decaying MHD turbulence in the early Universe. The Fourier-transformed stress is taken to be $|\tilde{T}(k,t)|$, i.e., we ignore changes in the phase, and its time evolution is parameterized 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 \Sec{EvolutionOfStress}. 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 \citep{Roga_etal17,Schober+18}. 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. \begin{acknowledgements} RS would like to thank Hongzhe Zhou for helping him analyzing the output data with the Mathematica routines of the {\sc Pencil Code}. This work was supported by the Swedish Research Council %(Vetenskapsr{\aa}det, 2019-04234). (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\"oping. \end{acknowledgements} \vspace{2mm} {\bf Data availability}---The source code used for the simulations of this study, the {\sc Pencil Code}, is freely available from Refs.~\cite{JOSS}. \blue{ Supplemental Material with plots similar to \Figs{rslice_stress_plot_LowFreq}{rslice_stress_plot_LowFreq_nohel}, but for other wave vectors of the same length, along with the} simulation setups and the corresponding data are freely available from Ref.~\cite{DATA}. \appendix \section{$\Sp(\TT)$ in terms of magnetic spectrum} \label{appendixa} Here, we provide the expressions for the stress spectrum in terms of the magnetic spectrum for helical and nonhelical cases; see Refs.~\cite{CDK04} and \cite{Sharma+20} for the derivation. The two-point correlation for the nonhelical magnetic field in Fourier space is given by, \begin{align} \langle \tilde{B}_{i}(\kk)\tilde{B}^*_{j}(\kk') \rangle &= (2 \pi)^3 \delta(\kk-\kk')(\delta_{ij}-\hat{\kk}_i \hat{\kk}_j)P_{\rm SM}(k) \end{align} Assuming a Gaussian nature of the magnetic field fluctuations, the stress spectrum is given by \begin{align} \Sp(\TT)\equiv &4 \pi k^2\langle\tildeT_{ij}^{\rm TT}(\kk)\tildeT_{*ij}^{\rm TT}(\kk)\rangle=4\pi k^2\int d^3 q \Big[ P_{\rm SM}(q)\nonumber\\ &P_{\rm SM}({|\kk-\qq|}) (1+\gamma^2+\beta^2+\gamma^2 \beta^2)\Big]. \end{align} Here, $\tildeT_{ij}^{\rm TT}(\kk)=-B_i B_j+\frac{1}{2}\delta_{ij}B_kB_k$. In terms of the energy spectrum, $\EM(k)\equiv 4\pi k^2 P_{\rm SM}(k)$, the above expression reduces to \begin{align}\label{tijnh} \Sp(\TT)&=\frac{1}{4\pi}\int d^3 q \frac{k^2}{q^2 |\kk-\qq|^2} \Big[ \EM(q)\nonumber\\ &\EM({|\kk-\qq|}) (1+\gamma^2+\beta^2+\gamma^2 \beta^2)\Big]. \end{align} In the above expressions $\gamma=\hat{\kk}\cdot \hat{\qq}$ and $\beta=\hat{\kk}\cdot \widehat{\kk-\qq}$. For helical magnetic fields, there is an additional antisymmetric contribution to the two-point correlation and is given by, \begin{align} \langle \tilde{B}_{i}(\kk)\tilde{B}^*_{j}(\kk') \rangle&= (2 \pi)^3 \delta(\kk-\kk')\Big((\delta_{ij}-\hat{\kk}_i \hat{\kk}_j)P_{\rm SM}(k)\nonumber\\ &+i \epsilon_{ijm}\hat{\kk}_m P_{AM}(k) \Big). \end{align} The stress spectrum for this case is given by \begin{align} \Sp(\TT)=&4\pi k^2 \int d^3 q \Big[ P_{\rm SM}(\qq)P_{\rm SM}({|\kk-\qq|})(1+\gamma^2+\beta^2+\gamma^2 \beta^2)\nonumber\\ &+4 \gamma \beta P_{\rm AM}(q)P_{\rm AM}({|\kk-\qq|})\Big]. \end{align} In terms of the energy spectrum, $\EM(k)$, and the helicity spectrum, $H_M(k)\equiv 4\pi k^2 P_{\rm AM}(k)$, the above expression reduces to \begin{align}\label{tijh} \Sp(\TT)=&\frac{1}{4\pi}\int d^3 q \frac{k^2}{q^2 |\kk-\qq|^2} \Big[ \EM(\qq)\EM({|\kk-\qq|})\nonumber\\&(1+\gamma^2+\beta^2+\gamma^2 \beta^2)+4 \gamma \beta H_M(q)H_M({|\kk-\qq|})\Big]. \end{align} \bibliography{ref} %\bibliographystyle{apsrev4-1} %\bibliographystyle{unsrt} \end{document}