\documentclass[prd,twocolumn,aps,psfig,showpacs,nofootinbib,nobibnotes,superscriptaddress,preprintnumbers,times]{revtex4}
\setlength{\topmargin}{-14mm}
\usepackage{graphicx,bm,color,braket}
\graphicspath{{./fig/}{./png/}}
\def\red{\textcolor{red}}
\def\blue{\textcolor{blue}}

%|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
%             Customized Commands
%|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
%  mathematical abbreviations
%
%
%
\newcommand{\BoldVec}[1]{\mathchoice%
  {\mbox{\boldmath $\displaystyle     #1$}}%
  {\mbox{\boldmath $\textstyle        #1$}}%
  {\mbox{\boldmath $\scriptstyle      #1$}}%
  {\mbox{\boldmath $\scriptscriptstyle#1$}}%
}
%\newcommand{\BoldVec}[1]{\bm{#1}}}
%
% math debs
\newcommand{\EQ}{\begin{equation}}
\newcommand{\EN}{\end{equation}}
\newcommand{\EQA}{\begin{eqnarray}}
\newcommand{\ENA}{\end{eqnarray}}
\newcommand{\eq}[1]{(\ref{#1})}
\newcommand{\EEq}[1]{Equation~(\ref{#1})}
\newcommand{\Eq}[1]{Eq.~(\ref{#1})}
\newcommand{\Eqs}[2]{Eqs~(\ref{#1}) and~(\ref{#2})}
\newcommand{\EEqs}[2]{Equations~(\ref{#1}) and~(\ref{#2})}
\newcommand{\eqs}[2]{(\ref{#1}) and~(\ref{#2})}
\newcommand{\Eqss}[2]{Eqs~(\ref{#1})--(\ref{#2})}
%\newcommand{\Sec}[1]{\S\,\ref{#1}}
%\newcommand{\Secs}[2]{\S\S\,\ref{#1} and~\ref{#2}}
\newcommand{\Sec}[1]{Sec.~\ref{#1}}
\newcommand{\Secs}[2]{Secs.~\ref{#1} and~\ref{#2}}
\newcommand{\App}[1]{Appendix~\ref{#1}}
\newcommand{\Fig}[1]{Fig.~\ref{#1}}
\newcommand{\FFig}[1]{Figure~\ref{#1}}
\newcommand{\Tab}[1]{Table~\ref{#1}}
\newcommand{\Figs}[2]{Figs.~\ref{#1} and \ref{#2}}
\newcommand{\Tabs}[2]{Tables~\ref{#1} and \ref{#2}}
%\newcommand{\bra}[1]{\langle #1\rangle}
\newcommand{\bbra}[1]{\left\langle #1\right\rangle}
\newcommand{\mean}[1]{\overline #1}
\newcommand{\meanB}{\overline{B}}
\newcommand{\meanC}{\overline{C}}
\newcommand{\meanU}{\overline{U}}
\newcommand{\meanW}{\overline{W}}
\newcommand{\meanPhi}{\overline{\Phi}}
\newcommand{\meanF}{\overline{\cal F}}
\newcommand{\meanR}{\overline{\cal R}}
\newcommand{\meanAA}{\overline{\bm{A}}}
\newcommand{\meanBB}{\overline{\bm{B}}}
\newcommand{\meanEE}{\overline{\bm{E}}}
\newcommand{\meanUU}{\overline{\bm{U}}}
\newcommand{\meanWW}{\overline{\bm{W}}}
\newcommand{\meanJJ}{\overline{\mbox{\boldmath $J$}}}
\newcommand{\meanuu}{\overline{\mbox{\boldmath $u$}}}
\newcommand{\meanGG}{\overline{\mbox{\boldmath $G$}}}
\newcommand{\meanAB}{\overline{\mbox{\boldmath $A\cdot B$}}}
\newcommand{\meanAoBo}{\overline{\mbox{\boldmath $A_0\cdot B_0$}}}
\newcommand{\meanApoBpo}{\overline{\mbox{\boldmath $A'_0\cdot B'_0$}}}
\newcommand{\meanApBp}{\overline{\mbox{\boldmath $A'\cdot B'$}}}
\newcommand{\meanuxB}{\overline{\mbox{\boldmath $\delta u\times \delta B$}}}
\newcommand{\meanemfs}{\overline{\cal E} {}}
\newcommand{\meanemf}{\overline{\mbox{\boldmath ${\cal E}$}} {}} %redundant
\newcommand{\meanAAAA}{\overline{\mbox{\boldmath ${\mathsf A}$}} {}}
\newcommand{\meanSSSS}{\overline{\mbox{\boldmath ${\mathsf S}$}} {}}
\newcommand{\meanAAA}{\overline{\mathsf{A}}}
\newcommand{\meanSSS}{\overline{\mathsf{S}}}
\newcommand{\meanCC}{\overline{\mbox{\boldmath ${\cal C}$}} {}}
\newcommand{\meanFF}{\overline{\mbox{\boldmath ${\cal F}$}} {}}
\newcommand{\meanRR}{\overline{\mbox{\boldmath ${\cal R}$}} {}}
\newcommand{\calFF}{\overline{\mbox{\boldmath ${\cal F}$}} {}}
\newcommand{\meanEMF}{\overline{\mbox{\boldmath ${\cal E}$}} {}}
\newcommand{\tildeFFFF}{\tilde{\mbox{\boldmath ${\cal F}$}}{}}{}
\newcommand{\hatFFFF}{\hat{\mbox{\boldmath ${\cal F}$}}{}}{}
\newcommand{\meanFFFF}{\overline{\mbox{\boldmath ${\cal F}$}}{}}{}
\newcommand{\meanFFF}{\overline{\cal F}}
\newcommand{\hatOO}{\hat{\bm{\Omega}}}
\newcommand{\hatAA}{\hat{\bm{A}}}
\newcommand{\hatBB}{\hat{\bm{B}}}
%
% tilde
%
\newcommand{\tildexx}{\tilde{\bm{x}}}
\newcommand{\tildeBB}{\tilde{\bm{B}}}
\newcommand{\tildeJJ}{\tilde{\bm{J}}}
\newcommand{\tildeA}{\tilde{A}}
\newcommand{\tildeB}{\tilde{B}}
\newcommand{\tildeJ}{\tilde{J}}
\newcommand{\tildeemf}{\tilde{\cal E}}
\newcommand{\teps}{\tilde{\epsilon} {}}
\newcommand{\tkapz}{\tilde{\kappa_0}}
\newcommand{\Oh}{\hat{\Omega}}
\newcommand{\zh}{\hat{z}}
\newcommand{\PC}{{\sc Pencil Code}~}
\newcommand{\PCS}{{\sc Pencil Code}}
%
%  unit vectors
%
\newcommand{\nullvector}{{\bf0}}
\newcommand{\nnn}{\hat{\mbox{\boldmath $n$}} {}}
\newcommand{\vvv}{\hat{\mbox{\boldmath $v$}} {}}
\newcommand{\rr}{\hat{\mbox{\boldmath $r$}} {}}
\newcommand{\xxx}{\hat{\mbox{\boldmath $x$}} {}}
\newcommand{\yyy}{\hat{\mbox{\boldmath $y$}} {}}
\newcommand{\zz}{\hat{\mbox{\boldmath $z$}} {}}
\newcommand{\pp}{\hat{\mbox{\boldmath $\phi$}} {}}
\newcommand{\ttt}{\hat{\mbox{\boldmath $\theta$}} {}}
\newcommand{\OOO}{\hat{\mbox{\boldmath $\Omega$}} {}}
\newcommand{\ooo}{\hat{\mbox{\boldmath $\omega$}} {}}
\newcommand{\BBBB}{\hat{\mbox{\boldmath $B$}} {}}
\newcommand{\kunit}{\hat{\mbox{$k$}} {}}
\newcommand{\nunit}{\hat{\mbox{$n$}} {}}
%
%  hatted quantities
%
\newcommand{\hatU}{\hat{U}}
\newcommand{\hatUU}{\hat{\bm{U}}}
%
%  vectors
%
\newcommand{\gggg}{\BoldVec{g} {}}
\newcommand{\ddd}{\BoldVec{d} {}}
\newcommand{\rrr}{\BoldVec{r} {}}
\newcommand{\xx}{\BoldVec{x}{}}
\newcommand{\yy}{\BoldVec{y} {}}
\newcommand{\zzz}{\BoldVec{z} {}}
\newcommand{\uu}{\BoldVec{u} {}}
\newcommand{\vv}{\BoldVec{v} {}}
\newcommand{\ww}{\BoldVec{w} {}}
\newcommand{\mm}{\BoldVec{m} {}}
\newcommand{\PP}{\BoldVec{P} {}}
\newcommand{\QQ}{\BoldVec{Q} {}}
\newcommand{\RR}{\BoldVec{R} {}}
\newcommand{\UU}{\BoldVec{U} {}}
\newcommand{\bb}{\BoldVec{b} {}}
\newcommand{\qq}{\BoldVec{q} {}}
\newcommand{\BB}{\BoldVec{B} {}}
\newcommand{\HH}{\BoldVec{H} {}}
\newcommand{\II}{\BoldVec{I} {}}
\newcommand{\AAA}{\BoldVec{A} {}}
\newcommand{\aaa}{\BoldVec{a} {}}
\newcommand{\aaaa}{\BoldVec{a} {}} %(convert aaa -> aaaa, compatibility problem)
\newcommand{\eee}{\BoldVec{e} {}}
\newcommand{\jj}{\BoldVec{j} {}}
\newcommand{\JJ}{\BoldVec{J} {}}
\newcommand{\nn}{\BoldVec{n} {}}
\newcommand{\ee}{\BoldVec{e} {}}
\newcommand{\ff}{\BoldVec{f} {}}
\newcommand{\hh}{\BoldVec{h} {}}
\newcommand{\EE}{\BoldVec{E} {}}
\newcommand{\FF}{\BoldVec{F} {}}
\newcommand{\TT}{\BoldVec{T} {}}
\newcommand{\CC}{\BoldVec{C} {}}
\newcommand{\KK}{\BoldVec{K} {}}
\newcommand{\MM}{\BoldVec{M} {}}
\newcommand{\GG}{\BoldVec{G} {}}
\newcommand{\kk}{\BoldVec{k} {}}
\newcommand{\SSS}{\BoldVec{S} {}}
\newcommand{\grav}{\BoldVec{g} {}}
\newcommand{\nab}{\BoldVec{\nabla} {}}
\newcommand{\OO}{\BoldVec{\Omega} {}}
\newcommand{\oo}{\BoldVec{\omega} {}}
\newcommand{\LL}{\BoldVec{\Lambda} {}}
\newcommand{\llambda}{\BoldVec{\lambda} {}}
\newcommand{\pomega}{\BoldVec{\varpi} {}}
%
%  correlation tensors
%
\newcommand{\RRRR}{\bm{\mathsf{R}}}
\newcommand{\SSSS}{\bm{\mathsf{S}}}
\newcommand{\LLLL}{\mbox{\boldmath ${\sf L}$} {}}
\newcommand{\MMMM}{\bm{\mathsf{M}}}
\newcommand{\BBB}{\mbox{\boldmath ${\cal B}$} {}}
\newcommand{\emf}{\mbox{\boldmath ${\cal E}$} {}}
\newcommand{\FFF}{\mbox{\boldmath ${\cal F}$} {}}
\newcommand{\GGG}{\mbox{\boldmath ${\cal G}$} {}}
\newcommand{\HHH}{\mbox{\boldmath ${\cal H}$} {}}
\newcommand{\QQQ}{\mbox{\boldmath ${\cal Q}$} {}}
%
%  operators  (roman)
%
\newcommand{\ii}{{\rm i}}
\newcommand{\grad}{{\rm grad} \, {}}
\newcommand{\curl}{{\rm curl} \, {}}
\newcommand{\dive}{{\rm div}  \, {}}
\newcommand{\Dive}{{\rm Div}  \, {}}
\newcommand{\diag}{{\rm diag}  \, {}}
\newcommand{\DD}{{\rm D} {}}
\newcommand{\dd}{{\rm d} {}}
\newcommand{\const}{{\rm const}  {}}
\newcommand{\crit}{{\rm crit}  {}}
\def\degr{\hbox{$^\circ$}}
\def\la{\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil
$\displaystyle##$\hfil\cr<\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\textstyle##$\hfil\cr<\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\scriptstyle##$\hfil\cr<\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\scriptscriptstyle##$\hfil\cr<\cr\sim\cr}}}}}
\def\ga{\mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil
$\displaystyle##$\hfil\cr>\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\textstyle##$\hfil\cr>\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\scriptstyle##$\hfil\cr>\cr\sim\cr}}}
{\vcenter{\offinterlineskip\halign{\hfil$\scriptscriptstyle##$\hfil\cr>\cr\sim\cr}}}}}
%
%  numbers
%
\def\Ta{\mbox{\rm Ta}}
\def\Ra{\mbox{\rm Ra}}
\def\Ma{\mbox{\rm Ma}}
\def\Sh{\mbox{\rm Sh}}
\def\Roo{\mbox{\rm Ro}^{-1}}
\def\Pra{\mbox{\rm Pr}}
\def\Pran{\mbox{\rm Pr}}
\def\Pm{\mbox{\rm Pr}_{\rm M}}
\def\Rm{\mbox{\rm Re}_{\rm M}}
\def\Rey{\mbox{\rm Re}}
\def\Imag{\mbox{\rm Im}}
\def\Pe{\mbox{\rm Pe}}
\def\epsK{\epsilon_{\rm K}}
\def\epsM{\epsilon_{\rm M}}
\def\EEi{{\cal E}_i}
\def\EEK{{\cal E}_{\rm K}}
\def\EEM{{\cal E}_{\rm M}}
\def\EEKM{{\cal E}_{\rm K/M}}
\def\EEGW{{\cal E}_{\rm GW}}
\def\OmK{{\Omega}_{\rm K}}
\def\OmM{{\Omega}_{\rm M}}
\def\OmGW{{\Omega}_{\rm GW}}
\def\hrms{{h}_{\rm rms}}
\def\EEtot{{\cal E}_{\rm tot}}
\def\EErad{{\cal E}_{\rm rad}}
\def\EElam{{\cal E}_\lambda}
\def\EEcrit{{\cal E}_{\rm crit}}
\def\HHGW{{\cal H}_{\rm GW}}
\def\HHK{{\cal H}_{\rm K}}
\def\HHM{{\cal H}_{\rm M}}
\def\EGW{E_{\rm GW}}
\def\HGW{H_{\rm GW}}
\def\EK{E_{\rm K}}
\def\EM{E_{\rm M}}
\def\HM{H_{\rm M}}
\def\hc{h_{\rm c}}
\def\cs{c_{\rm s}}
\def\xiM{\xi_{\rm M}}
\def\xiK{\xi_{\rm K}}
\def\kf{k_{\rm f}}
%\def\kf{k_\ast}
\def\vA{v_{\rm A}}
\def\urms{u_{\rm rms}}
\def\Urms{U_{\rm rms}}
\def\kappaOO{\kappa_{\Omega\Omega}}
\def\kappaO{\kappa_{\Omega}}
\def\kappat{\kappa_{\rm t}}
\def\kappatz{\kappa_{\rm t0}}
\def\nut{\nu_{\rm t}}
\def\etatz{\eta_{\rm t0}}
\def\etat{\eta_{\rm t}}
\def\etaT{\eta_{\rm T}}
\def\Beq{B_{\rm eq}}
\def\tmax{t_{\max}}
%
\newcommand{\ea}{{\em et al. }}
\newcommand{\eaa}{{\em et al. }}
\def\half{{\textstyle{1\over2}}}
\def\threehalf{{\textstyle{3\over2}}}
\def\onethird{{\textstyle{1\over3}}}
\def\twothird{{\textstyle{2\over3}}}
\def\fourthird{{\textstyle{4\over3}}}
\def\quarter{{\textstyle{1\over4}}}
%
\newcommand{\W}{\,{\rm W}}
\newcommand{\V}{\,{\rm V}}
\newcommand{\kV}{\,{\rm kV}}
\newcommand{\MeV}{\,{\rm MeV}}
\newcommand{\GeV}{\,{\rm GeV}}
\newcommand{\T}{\,{\rm T}}
\newcommand{\uG}{\,\mu{\rm G}}
\newcommand{\G}{\,{\rm G}}
\newcommand{\Hz}{\,{\rm Hz}}
\newcommand{\mHz}{\,{\rm mHz}}
\newcommand{\nHz}{\,{\rm nHz}}
\newcommand{\uHz}{\,\mu{\rm Hz}}
\newcommand{\kHz}{\,{\rm kHz}}
\newcommand{\kG}{\,{\rm kG}}
\newcommand{\K}{\,{\rm K}}
\newcommand{\g}{\,{\rm g}}
\newcommand{\s}{\,{\rm s}}
\newcommand{\ms}{\,{\rm ms}}
\newcommand{\cm}{\,{\rm cm}}
\newcommand{\m}{\,{\rm m}}
\newcommand{\km}{\,{\rm km}}
\newcommand{\kms}{\,{\rm km/s}}
\newcommand{\kg}{\,{\rm kg}}
\newcommand{\Mm}{\,{\rm Mm}}
\newcommand{\pc}{\,{\rm pc}}
\newcommand{\kpc}{\,{\rm kpc}}
\newcommand{\Mpc}{\,{\rm Mpc}}
\newcommand{\yr}{\,{\rm yr}}
\newcommand{\Myr}{\,{\rm Myr}}
\newcommand{\Gyr}{\,{\rm Gyr}}
\newcommand{\erg}{\,{\rm erg}}
\newcommand{\mol}{\,{\rm mol}}
\newcommand{\dyn}{\,{\rm dyn}}
\newcommand{\J}{\,{\rm J}}
\newcommand{\RM}{\,{\rm RM}}
\newcommand{\AU}{\,{\rm AU}}
\newcommand{\A}{\,{\rm A}}
%
\def\hX{h_\times}
\def\hT{h_+}
\def\thT{\tilde{h}_+}
\def\thX{\tilde{h}_\times}
\def\dhT{\dot{h}_+}
\def\dhX{\dot{h}_\times}
\def\dhhT{\dot{\hat{h}}_+}
\def\dhhX{\dot{\hat{h}}_\times}
\def\dhhTX{\dot{\hat{h}}_{+/\times}}
\def\dthT{\dot{\tilde{h}}_+}
\def\dthX{\dot{\tilde{h}}_\times}
\def\dthTX{\dot{\tilde{h}}_{+/\times}}
%
%  journals
%
\newcommand{\arXiv}[3]{, ``#3,'' arXiv:#2 (#1).}
\newcommand{\yjcap}[3]{, J.\ Cosmol.\ Astropart.\ Phys. {\bf #2} (#1) #3.}
\newcommand{\yjas}[3]{, J. Atmosph. Sci. {\bf #2}, #3 (#1).}
\newcommand{\yan}[3]{, Astron. Nachr. {\bf #2}, #3 (#1).}
\newcommand{\yact}[3]{, Acta Astron. {\bf #2}, #3 (#1).}
\newcommand{\yana}[3]{, Astron. Astrophys. {\bf #2}, #3 (#1).}
\newcommand{\yanas}[3]{, Astron. Astrophys. Suppl. {\bf #2}, #3 (#1).}
\newcommand{\yanal}[3]{, Astron. Astrophys. Lett. {\bf #2}, #3 (#1).}
\newcommand{\yass}[3]{, Astrophys. Spa. Sci. {\bf #2}, #3 (#1).}
\newcommand{\ysci}[3]{, Science {\bf #2}, #3 (#1).}
\newcommand{\ysph}[3]{, Solar Phys. {\bf #2}, #3 (#1).}
\newcommand{\yjetp}[3]{, Sov. Phys. JETP {\bf #2}, #3 (#1).}
\newcommand{\yspd}[3]{, Sov. Phys. Dokl. {\bf #2}, #3 (#1).}
\newcommand{\ysov}[3]{, Sov. Astron. {\bf #2}, #3 (#1).}
\newcommand{\ysovl}[3]{, Sov. Astron. Lett. {\bf #2}, #3 (#1).}
\newcommand{\ymn}[3]{, Mon.\ Not.\ R.\ Astron.\ Soc.\ {\bf #2}, #3 (#1).}
\newcommand{\ymhd}[3]{, Magnetohydrohydrodyn. {\bf #2}, #3 (#1).}
\newcommand{\yqjras}[3]{, Quart. J. Roy. Astron. Soc. {\bf #2}, #3 (#1).}
\newcommand{\ynat}[3]{, Nature {\bf #2}, #3 (#1).}
\newcommand{\yjfm}[3]{, J. Fluid Mech. {\bf #2}, #3 (#1).}
\newcommand{\pjfm}[1]{, J. Fluid Mech., in press (#1).}
\newcommand{\sjfm}[1]{, J. Fluid Mech., submitted (#1).}
\newcommand{\ypr}[3]{, Phys.\ Rev.\ {\bf #2}, #3 (#1).}
\newcommand{\yprd}[4]{, ``#4,'' Phys.\ Rev.\ D {\bf #2}, #3 (#1).}
\newcommand{\ypre}[3]{, Phys.\ Rev.\ E {\bf #2}, #3 (#1).}
\newcommand{\yprf}[4]{, ``#4,'' Phys.\ Rev.\ Fluids {\bf #2}, #3 (#1).}
\newcommand{\yprl}[3]{, Phys.\ Rev.\ Lett.\ {\bf #2}, #3 (#1).}
\newcommand{\yphl}[3]{, Phys.\ Lett.\ {\bf #2}, #3 (#1).}
\newcommand{\pprl}[1]{, Phys. Rev. Lett., in press (#1).}
\newcommand{\yepl}[3]{, Europhys. Lett. {\bf #2}, #3 (#1).}
\newcommand{\pcsf}[2]{, Chaos, Solitons \& Fractals, in press (#1).}
\newcommand{\ycsf}[3]{, Chaos, Solitons \& Fractals{\bf #2}, #3 (#1).}
\newcommand{\yprs}[3]{, Proc. Roy. Soc. Lond. {\bf #2}, #3 (#1).}
\newcommand{\yptrs}[3]{, Phil. Trans. Roy. Soc. {\bf #2}, #3 (#1).}
\newcommand{\yptrsa}[4]{, ``#4,'' Phil. Trans. Roy. Soc. Lond. A, {\bf #2}, #3 (#1).}
\newcommand{\yjcp}[3]{, J. Comp. Phys. {\bf #2}, #3 (#1).}
\newcommand{\yjgr}[3]{, J. Geophys. Res. {\bf #2}, #3 (#1).}
\newcommand{\ygrl}[3]{, Geophys. Res. Lett. {\bf #2}, #3 (#1).}
\newcommand{\yobs}[3]{, Observatory {\bf #2}, #3 (#1).}
\newcommand{\yaj}[3]{, Astronom. J. {\bf #2}, #3 (#1).}
\newcommand{\sapj}[3]{, ``#3,'' Astrophys. J., submitted, arXiv:#2  (#1).}
\newcommand{\papj}[3]{, ``#3,'' Astrophys. J., in press, arXiv:#2  (#1).}
\newcommand{\yapj}[3]{, Astrophys. J. {\bf #2}, #3 (#1).}
\newcommand{\yapjs}[3]{, Astrophys. J. Suppl. {\bf #2}, #3 (#1).}
\newcommand{\yapjl}[3]{, Astrophys. J. {\bf #2}, #3 (#1).}
\newcommand{\ycqg}[3]{, Class. Quant. Grav. {\bf #2}, #3 (#1).}
\newcommand{\ypp}[3]{, Phys. Plasmas {\bf #2}, #3 (#1).}
\newcommand{\yppcf}[3]{, Plasmas Phys. Contr. Fusion {\bf #2}, #3 (#1).}
\newcommand{\ppp}[1]{, Phys. Plasmas, in press (#1).}
\newcommand{\ypasj}[3]{, Publ. Astron. Soc. Japan {\bf #2}, #3 (#1).}
\newcommand{\ypac}[3]{, Publ. Astron. Soc. Pacific {\bf #2}, #3 (#1).}
\newcommand{\yaraa}[3]{, Ann. Rev. Astron. Astrophys. {\bf #2}, #3 (#1).}
\newcommand{\yanar}[3]{, Astron. Astrophys. Rev. {\bf #2}, #3 (#1).}
\newcommand{\yanp}[3]{, Ann. Phys. {\bf #2}, #3 (#1).}
\newcommand{\yanf}[3]{, Ann. Rev. Fluid Dyn. {\bf #2}, #3 (#1).}
\newcommand{\ypf}[4]{, ``#4,'' Phys. Fluids {\bf #2}, #3 (#1).}
\newcommand{\yphy}[3]{, Physica {\bf #2}, #3 (#1).}
\newcommand{\ygafd}[4]{, ``#4,'' Geophys. Astrophys. Fluid Dyn. {\bf #2}, #3 (#1).}
\newcommand{\yrpp}[3]{, Rep. Prog. Phys. {\bf #2}, #3 (#1).}
\newcommand{\yptp}[3]{, Progr. Theor. Phys. {\bf #2}, #3 (#1).}
\newcommand{\yjour}[5]{, ``#5,'' #2 {\bf #3}, #4 (#1).}
\newcommand{\pjour}[3]{, #2, in press (#1).}
\newcommand{\sjour}[3]{, #2, submitted (#1).}
\newcommand{\yprep}[2]{, #2, preprint (#1).}
\newcommand{\pproc}[3]{, (ed. #3), #2 (#1) (to appear).}
\newcommand{\yproc}[4]{, (ed. #4), pp. #2. #3 (#1).}
\newcommand{\ybook}[3]{, {\em #2}. #3 (#1).}

\begin{document}

\title{Can we observe the QCD phase transition-generated gravitational waves through pulsar timing arrays?}

\date{\today}
\preprint{NORDITA-2021-XXX}


\author{Axel~Brandenburg}
\email{brandenb@nordita.org}
\affiliation{Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden}
\affiliation{Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden} 
\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia}
\affiliation{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA}

\author{Emma Clarke \footnote{Corresponding author; the authors are listed alphabetically.}}
\email{emmaclar@andrew.cmu.edu}
\affiliation{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA}

\author{Yutong He}
\email{yutong.he@su.se}
\affiliation{Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden}
\affiliation{Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden}


\author{Tina~Kahniashvili}
\email{tinatin@andrew.cmu.edu}
\affiliation{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA}
\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia}
\affiliation{Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia}
\affiliation{Department of Physics, Laurentian University, Sudbury, ON P3E 2C, Canada}

\begin{abstract}
We perform numerical simulations of gravitational waves (GWs) 
induced by hydrodynamic and hydromagnetic
turbulent sources that might have been present at cosmological 
quantum chromodynamic (QCD) phase transitions. 
For turbulent energies of about 4\% of the radiation energy density,
the typical scale of such motions may have been a sizable fraction of
the Hubble scale at that time.
The resulting GWs are found to have an energy fraction of about $10^{-9}$
of the critical energy density in the nHz range today and may already
have been observed by the NANOGrav collaboration.
This is further made possible by our findings of shallower spectra 
proportional to the square root of the frequency for nonhelical
hydromagnetic turbulence.
This implies more power at low frequencies than for the steeper spectra
previously anticipated.
The behavior toward higher frequencies depends strongly on the nature of the turbulence.
For vortical hydrodynamic and hydromagnetic turbulence, there is a sharp drop of spectral GW energy
by up to five orders of magnitude in the presence of helicity, and somewhat
less in the absence of helicity.
For acoustic hydrodynamic turbulence, the sharp drop is replaced by a power law decay,
albeit with a rather steep slope.
Our study supports earlier findings of a quadratic scaling of the GW
energy with the magnetic energy of the turbulence and inverse quadratic 
scaling with the peak frequency, which leads to larger GW energies under QCD conditions.
\end{abstract}

\maketitle

\section{Introduction}
Gravitational wave (GWs) astronomy opens a new window to study the physical
processes in the early universe.
Relic GWs can be sourced by violent processes such as
cosmological phase transitions and after generation they propagate
almost freely throughout the expansion of the universe that causes
the dilution of their strain amplitude and frequency; for a review, see
\cite{Caprini:2018mtu} and references therein.
On the other hand, the detection of these relic GWs
is a challenging task due to their small amplitudes, the specific
range of the characteristic frequencies, and astrophysical foregrounds
\cite{Romano:2016dpx}. 
Despite tremendous advancements in 
GW detection techniques, the stochastic GW background 
of cosmological origin remained unobserved.

Recently, the NANOGrav collaboration reported strong evidence for a stochastic GW background \cite{Arzoumanian:2020vkk}. In addition to the possibility of GWs induced by astrophysical sources
such as supermassive black holes, the NANOGrav data can also be understood
as a possible signal from the early universe, such as inflationary GWs
\cite{Kuroyanagi:2020sfw,Tahara:2020fmn,Cai:2020qpu,Zhou:2020kkf,Vagnozzi:2020gtf},
cosmic strings and domain walls
\cite{Ellis:2020ena,Buchmuller:2020lbh,Samanta:2020cdk,
Liu:2020mru,Paul:2020wbz,Blasi:2020mfx}
phase transitions including the supercooled
phase transitions \cite{Lewicki:2020azd}, dark phase
transitions 
\cite{Nakai:2020oit,Addazi:2020zcj} and quantum
chromodynamic (QCD), with axionic string network and QCD axion
\cite{Kitajima:2020rpm,Ramberg:2020oct,Lee:2020wfn,Gorghetto:2021fsn}, 
and/or magnetic fields
\cite{Neronov:2020qrl} and turbulence \cite{Abe:2020sqb}.

In this paper we present a self-consistent study of the 
GWs from turbulent sources possibly present at QCD phase transitions.
In contrast to Ref.~\cite{Neronov:2020qrl}, we will not limit 
the turbulent sources to the presence of magnetic fields 
\cite{Quashnock:1988vs,Sigl:1996dm,Forbes:2000gr,Kisslinger:2005uy,Urban:2009sw,Buividovich:2009wi,Boyanovsky:2002wa} 
but rather consider any turbulent source possibly presented at the
QCD energy scale.
Even if primordial fields are not dynamically strong, 
turbulence can still develop at QCD energy scales 
\cite{Quashnock:1988vs,Rischke:2003mt,Kiskis:2003rd,Kisslinger:2005uy,Aoki:2006we,Blaizot:2013hx,Miniati:2017kah}; the latent 
heat they release still gives rise to pressure gradients resulting in
macroscopic plasma motions.
Given the very high Reynolds number of the primordial plasma, such
motions will inevitably decay into turbulence \cite{Quashnock:1988vs,Miniati:2017kah}.
Particularly important for our work is the earlier finding that
the separation and size of nucleation bubbles in a QCD phase
transition is a sizeable fraction of the Hubble scale, see 
Ref.~\cite{Witten:1984rs} for pioneering works and follow-up papers, 
\cite{Applegate:1985qt,Midorikawa:1985gy,Strumia:1998nf,Cottingham:1993rv,Hindmarsh:1991ay,Ignatius:1993qn,Ignatius:1997va,Strumia:1999fv,Schwarz:2003du,Tawfik:2011sh,Tawfik:2011mw}. 
Furthermore, the assumption of turbulence being driven by magnetic fields,
allows us to avoid the requirement of first order QCD phase transitions,
\cite{Miniati:2017kah}.

The paper is organized as follows.
We first review basic properties of relic GWs (\Sec{relicGWs}),
discuss then the NANOGrav observations (\Sec{NANOGravDATA}),
present our numerical approach (\Sec{ModelingGWs})
and results (\Sec{NumericalSimulations}) of our simulations,
before concluding in \Sec{Concl}.
Throughout the paper, we use natural units with  $\hbar=c=k_B=1$. 
We also set the permeability of free space to unity, i.e., $\mu_0=1$,
expressing the electromagnetic quantities in Lorentz-Heaviside units. 
The Latin indices run $i \in (1,2,3)$ and define the spatial coordinates, and the Greek indices run $\lambda \in (0,1,2,3)$. 
We choose the metric signature as $(-1,1,1,1)$. 

\section{The Early-Universe Gravitational Wave Signal}
\label{relicGWs}

GWs correspond to the tensor mode of perturbations $\delta g_{\mu\nu}$ above the spatially
flat, homogeneous, and isotropic Friedmann-Lema\^{i}tre-Robertson-Walker
(FLRW) background, in the transverse-traceless (TT) gauge\footnote{The 
TT gauge is determined by the TT projection tensor $\Lambda_{ijkl}
= P_{ik}P_{jl} - \frac{1}{2}P_{ij}P_{kl}$, where the $P_{ij}$ is a
transverse operator ($\partial_i P_{ij} =0$), defined as
$P_{ij} = \delta_{ij} - \partial_i \partial j/{\nab}^2$, where 
$\delta_{ij}$ is the Kronecker delta, $\partial_\lambda \equiv \partial/\partial x^\lambda$ denotes the partial derivative in respect of $x^\lambda$ coordinate, 
and $(\nab$ defines the vector differential operator with the components
equal to ${\nab}_i \equiv \partial_i$, i.e., ${ \nab}^2$ is the 
Laplacian in respect of spatial coordinates, 
see for more details Chapter 1 (1.2) of \cite{MaggioreText} 
} 
defined through the spatial component $h_{ij}^{\rm phys}$ with 
$a^2 h_{ij}^{\rm phys}= \Lambda_{ijlm} \delta g_{lm}$, where 
$a$ is the scale factor at the physical time $t_{\rm phys}$.
Here and below, super/subscript ``phys" denotes physical quantities.

In order to eliminate the expansion-induced dilution
from the governing hydromagnetic equations, we use rescaled 
quantities together with the conformal time $t$, defined through $dt=dt_{\rm phys}/a$, 
which reduces the metric tensor to the Minkowski form.
The background expansion of the universe during the radiation-dominated 
epoch is governed by the (dominant) radiation energy density
$\EErad = {\pi^2 g(T) T^4}/{30}$, where $g(T)$
is the effective number of relativistic degrees of freedom at temperature $T$. 
In the epoch(s) of interest, the expansion of the universe is fully governed by radiation, and 
the Hubble parameter $H \equiv a^{-1}da/dt_{\rm phys} = a^{-2}da/dt$ 
is given through $H^2 (t) = (8\pi G /3)\, \EErad (t)$,
where $G$ is Newton's gravitational constant and 
$\EErad(t)$ denotes the total energy density
of radiation (including all relativistic components).

In order to connect physical and comoving variables and to determine the
scaling of physical quantities, we compute the ratio of the scale factor 
today, $a_0=a(t_0)$ (here and below, ``0'' denotes the present moment),
to that at the time $t_*$ (at the temperature $T_*$ at which the source
becomes active and the gravitational signal is generated) corresponding
to the start of the simulation. We assume the adiabatic expansion of
the universe, such that $g_S(T)T^3a^3(T)$ is constant, where $g_S(T)$
is the number of adiabatic degrees of freedom at temperature $T$.
At high enough temperatures ($T>1\MeV$), we have $g_S(T)=g(T)$ \cite{K&T}.
Note, that our consideration below is valid for any time period 
during the radiation dominated epoch. However, we will be focused 
on the time period around the QCD energy scale ($150\MeV$). 
We also normalize the scale factor $a_* \equiv a(t_*)=1$,
which differs from the usual convention $a_0 =1$. 
Entropy conservation leads to
\begin{equation}
    \frac{a_0}{a_*} = 10^{12} \, \bigg(\frac{g_{S}(T_*)}{15}\bigg)^{1/3} 
    \bigg(\frac{T_*}{150\MeV}\bigg),
\label{Dfactor}
\end{equation} 
where we have used $T_0 = 2.73\K$ and $g_S(T_0) = 3.91$,
while at the QCD energy scale we have $g_S(T_*) \approx 15$ \cite{K&T}.
The degrees of freedom at QCD is approximate due to uncertainty in the exact temperature of the QCD transition and knowledge of the standard model (see discussions in \cite{Schwarz:2003du,Husdal:2016haj}). However, as $(a_0/a_*)\sim g_S(T_*)^{1/3}$, small deviations in $g_S(T_*)$ will not significantly impact our results. 

The GW equation in physical time and space coordinates is given by
\begin{equation}\label{eq:GWeqn1} 
    \big(\partial_{t_{\rm phys}}^2 + 3H\partial_{t_{\rm phys}}- \nab^2_\mathrm{phys}\big)  
    h^\mathrm{phys}_{ij} = 16\pi G  T^{\rm TT}_{ij,{\rm phys}},
\end{equation}
where the TT superscript denotes the TT projection 
of the stress-energy tensor 
such that $T^{\rm TT}_{ij,{\rm phys}} = \Lambda_{ijlm} T^{\rm phys}_{lm}$. 

To make the connection with observations, we define 
the characteristic strain, $h_c(t)$, which
obeys $h_c^2 (t)=\braket{(h_{ij}^{\rm phys}(\mathbf{x},t))^2}/2$, 
where angle brackets denote volume averaging in physical space, 
and the physical energy density $\EEGW^{{\rm phys}}(t)$ 
carried by the GWs is given by \cite{MaggioreText}
\begin{equation}
    \EEGW^{{\rm phys}}(t) = \frac{1}{32\pi G}\braket{(\partial_{t_\mathrm{phys}}h_{ij}^{\rm phys}(\mathbf{x},t))^2}.
\end{equation} 
It is then expressed in terms of today's
frequency $f = k / (2\pi a_0)$  that corresponds to the time Fourier transform
$Q(t)= \int_{-\infty}^\infty df \, Q(f) e^{-2\pi f t}$ 
(and $Q(f)= 2\pi \int_{-\infty}^\infty dt \, Q(t) e^{-2\pi f t}$) 
\cite{MaggioreText}.

The relic GW signal strength today is 
given through the normalized GW energy 
density parameter 
$\OmGW(f)$ reduced by the factor $(H_*/H_0)^2 (a_*/a_0)^4$, 
where $H_*$ is the Hubble parameter at $t_*$. 
This accounts for the dilution of the GW energy density parameter with
the expansion of the universe and renormalizes the GW energy density
by the critical energy density at the present time, $\EEcrit^0 =
(3H_0^2)/(8\pi G)$, where $H_0 = 100 \, h_0$ km s$^{-1}$ Mpc$^{-1}$
$\simeq 3.241 \times 10^{-18} \, h_0 \, \mathrm{s}^{-1}$ is the present
value of the Hubble parameter. 
A frequency of particular interest is 
the frequency $f_*$ corresponding to the Hubble horizon scale at $t_*$:
\begin{equation}
    f_* = \frac{a_*H_*}{a_0} \simeq (1.8 \times10^{-8}\;\mathrm{Hz})\bigg(\frac{g_*}{15}\bigg)^{1/6}\bigg(\frac{T_*}{150\MeV}\bigg).
\end{equation}

As discussed above, there is a variety of possible sources of a stochastic
GW background in the nHz frequency range, accessible to Pulsar Timing
Arrays (PTAs), see Section~\ref{NANOGravDATA} 
%TK since we discuss the limits in Sec. IV maybe should be noted here?
%AB: do you or Emma want to propose something?
for more details and
Ref.~\cite{Burke-Spolaor:2018bvk} for a review and references therein,
and these sources include a cosmic population of supermassive
black hole binaries (SMBHBs) \cite{Sesana:2008mz,Burke-Spolaor:2018bvk}, 
cosmic strings \cite{Sanidas:2012ee,Cutler:2013aja,Blanco-Pillado:2017rnf},
inflationary GWs \cite{Campeti:2020xwn,Vagnozzi:2020gtf},\footnote{The quantum
mechanical fluctuations during the inflationary epoch induces GWs 
via parametric resonance 
\cite{Grishchuk:1974ny,Rubakov:1982df,Starobinsky:1979ty}.}
and phase transitions in the early universe
\cite{Capozziello:2018qjs,Ramberg:2020oct,Caprini:2010xv,Kobakhidze:2017mru}. 
In particular, if we assume that GWs could be sourced by bubble collisions at a phase transition, 
we expect the frequency of the GWs to be related to the bubble size. 
We consider that the bubble length scale is the Hubble horizon $H_*$ 
at generation divided by the total number of phase transition bubbles $N_b$. 
Then, for the QCD phase transitions, the frequency is given by
\begin{equation}
    f_* \simeq (1.1 \times10^{-7}\;\mathrm{Hz})\bigg(\frac{g_*}{15}\bigg)^{1/6}\bigg(\frac{T_*}{150\MeV}\bigg)\bigg(\frac{N_b}{6}\bigg),
\end{equation}
where we have normalized to 6 bubbles expected at the QCD phase transition \cite{Schwarz:2003du}.

\section{NANOGrav data}
\label{NANOGravDATA}

A pulsar is a highly magnetized and rapidly rotating neutron star that
emits a beam of electromagnetic radiation along its magnetic axis \cite{Ruderman:1975ju}. 
The times of arrival (TOA) of these pulses are extremely regular and
can be predicted very accurately over long times \cite{Gold:1968}.
The presence of a GW passing between the observer and pulsar shifts
the pulse TOA proportional to the amplitude of the GW 
\cite{MaggioreText2}. 
By monitoring the fluctuations in the TOA of radio pulses from millisecond 
pulsars (see, e.g., Ref.~\cite{Goncharov:2020krd} and references therein) 
international PTA missions\footnote{The International Pulsar Timing
Array (IPTA) is a consortium of consortia, comprised of the European
Pulsar Timing Array (EPTA), the North American Nanohertz Observatory
for Gravitational Waves (NANOGrav), and the Parkes Pulsar Timing Array
(PPTA) \cite{PTA}.} aim to probe a stochastic GW background.

The maximum sensitivity of a PTA experiment is limited by the total observation time.
That is, the lowest detectable frequency is on the order of the inverse
of the time span of the data (e.g., $f \sim$ nHz for datasets spanning
$\sim$ 10 years) \cite{Hobbs:2009yy}.
Furthermore, data sampling (i.e., pulsars are usually observed on the order
of weeks \cite{MaggioreText2}) limits the maximum detectable frequency. The NANOGrav 12.5-year data is sensitive to GW frequencies between approximately
$2.5\nHz$ and $1 \uHz$ \cite{Brazier:2019mmu}.

PTA measurements typically characterize a stochastic GW background in
terms of its characteristic strain spectrum $h_c(f)$ fitted with a
power-law dependence on frequency
\cite{Arzoumanian:2020vkk},
\begin{equation}\label{eq:hfit}
    \hc(f)=A_\mathrm{CP}\bigg(\frac{f}{f_\mathrm{yr}}\bigg)^{\alpha_\mathrm{CP}},
\end{equation}
where the subscript ``CP" denotes a common-spectrum (CP) process (common to the observed pulsars), the spectral index $\alpha_\mathrm{CP}$
depends on the source of the stochastic GW background, and $A_\mathrm{CP}$
is the strain amplitude at a reference frequency of
$f_\mathrm{yr} = 1 \, \mathrm{yr}^{-1}$.
This choice of reference frequency is arbitrary and does not affect the
ability to detect a GW signal.

The energy density spectrum of the GW background today expressed in terms
of the characteristic strain spectrum is given by \cite{MaggioreText2}
\begin{equation}\label{eq:OmGWhc}
    \OmGW(f)=\frac{2\pi^2}{3H_0^2}f^2 h_c^2(f) = \OmGW^{\rm yr}\bigg(\frac{f}{f_\mathrm{yr}}\bigg)^{5-\gamma_\mathrm{CP}},
\end{equation} 
where we have used Eq.~(\ref{eq:hfit}) in the second term on the right-hand side, 
$\gamma_\mathrm{CP}=3-2\alpha_\mathrm{CP}$ and
$\OmGW^{\rm yr} \equiv 2\pi^2A_{\rm CP}^2f_{\rm yr}^2/(3H_0^2)$.
The quantity $h_0^2\OmGW(f)$ is typically considered in order to remove 
the uncertainty in the value of $H_0$.

The NANOGrav collaboration reports joint $A_\mathrm{CP}-\gamma_\mathrm{CP}$ posterior distributions \cite{Arzoumanian:2020vkk}.
Posteriors for a common-spectrum process in the NANOGrav 12.5-year data were recovered with four models:
free-spectrum, broken power law, 5-frequency power law, and 30-frequency power law. 
The fits were performed for frequencies $f \in [2.5\times10^{-9}, 7\times10^{-8}]\Hz$,
with the exception of the 5-frequency power law, which was fit to the five lowest frequency bins.
The four lowest frequency bins have the strongest response to the presence of a GW background 
(see Figure~13 of Ref.~\cite{Arzoumanian:2020vkk}).
Thus, the 5-frequency power law was fit within
the signal-dominated frequency range (approximately $f \in [2.5\times10^{-9}, 1.2\times10^{-8}]\Hz$).
Figure 1 of 
Ref.~\cite{Arzoumanian:2020vkk} shows the 1$\sigma$ and 2$\sigma$ posterior contours
for the amplitude $A_\mathrm{CP}$ and spectral slope $\gamma_\mathrm{CP}$.

Fig.~\ref{fig:NGcontours} shows the NANOGrav detection expressed
in terms of $h_0^2\Omega_\mathrm{GW}(f)$ as given by 
Eq.~(\ref{eq:OmGWhc}).
The shaded regions show the 2$\sigma$ confidence contours of the
$A_\mathrm{CP}-\gamma_\mathrm{CP}$ parameter space in terms of $f$ and
$h_0^2\Omega_\mathrm{GW}(f)$ for frequencies from 2.5--$100\nHz$
(i.e., the NANOGrav 12.5-year sensitivity range);
see Ref.~\cite{Arzoumanian:2020vkk} for more detail.

\begin{figure} 
\centering
\includegraphics[scale=0.8]{fig/NGcontours2021.eps}
\caption{NANOGrav 12.5-year data set $2\sigma$ confidence
contours for the posteriors of a common-spectrum process
(see Ref.~\cite{Arzoumanian:2020vkk} Figure 1) expressed in terms of
the GW energy density $h_0^2\Omega_\mathrm{GW}(f)$ and frequency $f$.
This is shown over the NANOGrav 12.5-year sensitivity range of 2.5--100 nHz.
The three models used to fit the process include a: broken power law
(blue), 5-frequency power law (orange), and 30-frequency power law
(green).}
\label{fig:NGcontours}
\end{figure}

\section{Gravitational Wave Generation} 
\label{ModelingGWs}

As mentioned in the introduction (see also Ref.~\cite{Burke-Spolaor:2018bvk}),
low frequency GWs can originate from various astrophysical foreground
sources (white dwarfs, SMBH mergers, etc),
and from relic sources related to inflation and cosmic strings, for
example, and in particular, from phase transition-generated turbulence
and primordial magnetic fields.
We focus here on the latter two.
Turbulence and/or magnetic fields would only be generated during a
limited amount of time before they would decay.
The decay process itself remains highly turbulent and could
affect GW production.
Let us therefore begin with some general remarks about turbulent decay. 

\subsection{Gravitational Waves from Turbulent Sources}

Using scaled quantities $h_{ij}=a h^{\rm phys}_{ij}$ and
$T_{ij}^{\rm TT}=a^4 T^{\rm TT}_{ij,{\rm phys}}$, together with 
$a(t) \propto t$
in the radiation dominated epoch, 
the GW equation takes the form
\begin{equation}\label{eq:GWeqn}
    \left(\partial_t^2 - \nabla^2\right)h_{ij} = \frac{16\pi G}{a} T^{\rm TT}_{ij}.
\end{equation}
To obtain the GW equation in Fourier (wavenumber)
space, we use the Fourier transforms and the polarization $r=(+,\times)$ decomposition of the 
tensor metric perturbations and stress energy tensor projected onto the TT gauge  (i.e.,
${\cal Q}_{ij}(\kk,t) = \sum\limits_{r=+,\times} e^r_{ij}(\hat{\kk}) Q_{ij}(\kk,t)$,
where $e^{+}_{ij}(\hat{\kk})$ and
$e^\times_{ij}(\hat{\kk})$ are the polarization tensors with 
$\hat{\kk}$ the unit vector, and $\kk = a \kk_{\rm phys}$ is the 
rescaled 
wavenumber).\footnote{
We use the spatial Fourier transform convention: 
${\cal Q}(\xx,t) = \int \frac{d^3\kk}{(2\pi)^3}e^{-i\kk\cdot\xx}{\cal Q}(\kk,t)$ 
and ${\cal Q}(\kk,t) = 
\int d^3\xx e^{-i\kk\cdot\xx}{\cal Q}(\xx,t)$. 
The transverse operator $P_{ij}$ in the Fourier
space is given $P_{ij}(\hat{\kk}) = \delta_{ij} -
\hat{k}_i\hat{j}_k$ 
and the TT projection operator is 
$\Lambda_{ijkl}(\hat{\kk}) = P_{ik}(\hat{\kk})P_{jl}(\hat{\kk})
-\frac{1}{2}P_{ij}(\hat{\kk})P_{kl}(\hat{\kk})$, correspondingly.
The polarization tensors  $e^{+}_{ij}(\hat{\kk})$ and $e^{\times}_{ij}(\hat{\kk})$ can be written as $e^+_{ij}(\kk) =
{\bf \hat{e}}^1_i{\bf \hat{e}}^1_j - {\bf \hat{e}}^2_i{\bf \hat{e}}^2_j$ and $e^\times_{ij}(\kk)
= {\bf \hat{e}}^1_i{\bf \hat{e}}^2_j + {\bf \hat{e}}^2_i{\bf \hat{e}}^1_j$, where ${\bf \hat{e}}^1$
and ${\bf \hat{e}}^2$ are unit vectors that are orthogonal to $\hat{\kk}$  
and each other; see 
Chapter 1 (1.2) of Ref.~\cite{MaggioreText}. 
and Ref.~\cite{Pol:2018pao} for further detail.} 

As in earlier work \cite{Pol:2018pao,Pol:2019yex}, we use normalized
conformal time, $\bar{t}=t/t_*$, where $t_*=H_*^{-1}$ is our starting
time, and $a_*=1$ has been chosen.
Therefore, $a=\bar{t}$.
We also use the scaled wave vector, 
$\bar{\kk}=\kk/H_*$, 
and a scaled 
normalized
stress, $\bar{T}_{+/\times}^{\rm TT} = T_{+/\times}^{\rm TT}/\EErad^*$. 
The GW equation can then be written in the form \cite{Pol:2018pao,Pol:2019yex}
\begin{equation}
\left( \partial^2_{\bar{t}} + \bar{\kk}^2\right)
h_{+/\times} (\kk, t)\, =\, {6 \over \bar{t}}
\bar{T}_{+/\times}^{\rm TT} (\kk, t),
\label{GW}
\end{equation}
but from now on we omit all overbars.

Throughout this paper, all numerical results will usually be presented
as the scaled variables introduced above.
In particular, we quote the rms strain, $\hrms=\braket{h^2}^{1/2}$,
where $h^2=h_+^2+h_\times^2=(h_{ij})^2/2$, and likewise for the
scaled GW energy, $\EEGW=\braket{\dot{h}^2}/6$, 
where $\dot{h}_{+/\times}=\partial_t h_{+/\times}$
with $\dot{h}^2\equiv\dot{h}_+^2+\dot{h}_\times^2$;
see Ref.~\cite{Pol:2018pao,Pol:2019yex} for additional
subdominant terms that are applied in the calculations. 
We sometimes also quote the (frequency dependent) characteristic
amplitude of the physical strain measured today, $h_c(f)=\hrms/a_0$;
see \Sec{relicGWs}.

\subsection{Turbulent Sources}

Turbulent flows in the early universe can be modeled by solving the
hydromagnetic equations for the density $\rho$, the velocity $\uu$,
and the magnetic field $\BB$ with $\nab\cdot\BB=0$, adopting an 
ultrarelativistic equation of state in an expanding universe
using conformal time and comoving variables
\cite{Brandenburg:1996fc,Brandenburg:2017neh} with a forcing
term $\FFF$ in the induction equation for $\BB$
\begin{eqnarray}
{\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 \\
{\partial\uu\over\partial t}&=&-\uu\cdot\nab\uu
+{\uu\over3}\left(\nab\cdot\uu+\uu\cdot\nab\ln\rho\right)
+{2\over\rho}\nab\cdot\left(\rho\nu\SSSS\right) \nonumber \\
&&-{1\over4}\nab\ln\rho -{\uu\over\rho}\left[\uu\cdot(\JJ\times\BB)+\eta\JJ^2\right] +{3\over4\rho}\JJ\times\BB,
\nonumber \\
{\partial\BB\over\partial t}&=&\nab\times(\uu\times\BB-\eta\JJ+\FFF),\quad
\JJ=\nab\times\BB.
\nonumber
\end{eqnarray}
We recall that the conformal time $t$ 
is normalized to unity at the time $t_*$
of magnetic field generation, $\rho$ is in units of the initial value,
$\uu$ is in units of the speed of light, and the magnetic energy density
$\BB^2/2$ is measured in units of the radiation density at the time
of generation.
Furthermore, 
${\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 with commas denoting
partial derivatives, $\JJ$ is the current density,
$\nu$ is the kinematic viscosity, and $\eta$ is the
magnetic diffusivity.
The electromotive force, $\FFF$, is used to model magnetic field
generation with
\begin{equation}
\FFF(\xx,t)={\rm Re}[{\cal N}\tilde{\ff}(\kk)\exp(i\kk\cdot\xx+i\varphi)],
\label{emf_model}
\end{equation}
where the wave vector $\kk(t)$ and the phase $\varphi(t)$ change randomly
from one time step to the next.
This forcing function is therefore white noise in time and consists of plane
waves with average wavenumber $\kf$ such that $|\kk|$ lies in an interval
$\kf-\delta k/2\leq|\kk|<\kf+\delta k/2$ of width $\delta k$.
Here, ${\cal N}=f_0/\delta t^{1/2}$ is a normalization factor, where
$\delta t$ is the time step and $f_0$ is varied to achieve a certain
magnetic field strength after a certain time, and
$\tilde{\ff}({\kk})=(\kk\times\ee)/[\kk^2-(\kk\cdot\ee)^2]^{1/2}$
is a nonhelical forcing function.
Here, $\ee$ is an arbitrary unit vector that is not aligned with $\kk$.
Note that $|\ff|^2=1$.
Following earlier work, the forcing is only enabled during the
time interval $1\leq t\leq 2$. 
The kinetic and magnetic energy densities are defined as
$\EEK(t) = \braket{\rho\mathbf{u}^2}/2$ and
$\EEM(t) = \braket{\mathbf{B}^2}/2$, respectively. 

\begin{figure*}[t!]\begin{center}
\includegraphics[width=.99\columnwidth]{pcomp512}
\includegraphics[width=.99\columnwidth]{pcomp512hel}
\end{center}\caption[]{
Evolution of $\EEM(t)$ and $\EEGW(t)$ for nonhelical (left) and
helical (right) cases. 
Orange, black, blue, and red are for $\kf=2$, 6, 20, and 60, respectively.
}\label{pcomp512}\end{figure*}

The vigor of turbulence is characterized by the Reynolds number,
$\Rey=\urms/\nu\kf$, where $\urms$ is the maximum rms velocity.
It can only be determined a posteriori from the velocity resulting
from the magnetic field through the Lorentz force.
For all our runs, we use $\eta=\nu$.

\subsection{Turbulent Decay Laws}

Turbulence is known to decay in power-law fashion \cite{BP56,Saf67} such that
the magnetic energy $\EEM(t)$ decays with time $t$ like $\Delta t^{-p}$
and the correlation length $\xiM(t)$ increases like $\Delta t^{q}$,
where $\Delta t=t-t_{\rm off}$ is the time interval after the forcing
has been turned off.
The exponents $p$ and $q$ are positive and depend on the physical
circumstances (magnetically or kinetically dominated turbulence),
and whether or not there is magnetic helicity.
In helical turbulence, for example, one finds $p=q=2/3$, while for
non-helical magnetically dominated turbulence one finds $p=1$ and
$q=1/2$, although other variants are sometimes possible
\cite{Brandenburg:2017rnt,Hosking21}.

In this paper, we are specifically interested in the dependence of
the decay behavior on the forcing wave number $\kf$ of the turbulence
while it was still being driven.
The parameter $\kf$ enters through the prefactor in the decay law.

Furthermore, $\Delta t^{-p}$ would become infinite for $p>0$
and $\Delta t=0$ (when $t=t_{\rm off}$).
The singularity of $\Delta t^{-p} $ at $t=0$ is a consequence of a
simplified description at the initial time moment. 
For this reason, it is convenient to express the decay laws as
\begin{equation}
\EEM(t)=\EEM^{\max}\,(1+\Delta t/\tau)^{-p},
\label{EEMt}
\end{equation}
where $\tau$ is the turnover time, which we will treat as an empirical
parameter that we expect to be of the order of $(\vA\kf)^{-1}$, where
$\vA=(3\EEM^{\max}/2)^{1/2}$ is the Alfv\'en speed, evaluated
at the time when $\EEM$ reaches its maximum value $\EEM^{\max}$.
In some simulations of purely hydrodynamic turbulence, we replace
$\EEM$ by $\EEK$ in \Eq{EEMt} and use $\tau=(\urms\kf)^{-1}$ with
$\urms=(2\EEK)^{1/2}$ as the nominal turnover time. 
We recall here that we are using nondimensional variables where the
radiation energy density is unity.

For all our simulations, we choose $t_{\rm off}=2$, i.e., turbulence
is being driven for one Hubble time during $1\leq t\leq2$.
In the following, we vary $\kf$ between 2 and 60.
For $\kf=60$, we find that $\tau$ is shorter than a Hubble time,
but in all other cases, it exceeds it by up to factors between
ten (in the nonhelical cases) and a hundred (in helical cases).

We arrange the strength of the forcing $f_0$ such that $\EEM$ is similar
for different values of $\kf$.
This allows us then to determine the resulting GW energy solely
as a function of $\kf$.
For small values of $\kf$, the turbulence may not be able to reach
a statistically steady state by the time $t_{\rm off}$, when the
driving is turned off.

It is therefore necessary to adjust $f_0$ for each value of $\kf$
separately.
Once we have two values of $\EEM^{\max}$ that are close enough
to the target strength, we determine the desired forcing
strength through linear interpolation.
We also consider the case of different values of $f_0$ for a fixed value
of $\kf$ (Runs~noh5,6 and Runs~hel5,6).

\section{Numerical simulations}
\label{NumericalSimulations}

We solve the governing equations using the {\sc Pencil Code} \cite{PC},
where the GW solver has already been implemented \cite{Pol:2018pao}. 
We consider a cubic domain of side length $2\pi/k_1$, where $k_1$ is
the smallest wave number in the domain.
We choose $k_1=\kf/6$, so that the scale separation between the initial
spectral peak and the lowest wave number in the domain is six. 
In the following, we discuss the results for different values of $\kf$.
The temporal growth of $\EEM(t)$ is similar for small values of $\kf$;
see \Fig{pcomp512}(a) and (b), where we compare the evolution of $\EEM$ and
$\EEGW$ for the nonhelical and helical cases.
The parameters of those runs are listed in \Tabs{Tsummary512}{Tsummary512hel}
(for nonhelical and helical runs).
The numerical resolution is $512^3$ mesh points, except for run noh1,
where we use $1024^3$ mesh points.
Unless specified otherwise, we use $\nu=\eta=5\times10^{-5}$.

\begin{table*}[t!]\caption{
Summary of runs with nonhelical turbulence.
}\vspace{12pt}\centerline{\begin{tabular}{ccccccccc|ccc}
Run & $\kf$ & $k_1$ & $f_0$ & $p$ & $\tau$ & $\EEM^{\max}$ & $\EEGW^{\rm sat}$ & $\hrms^{\rm sat}$ &
$B$ [$\mu$G] & $h_0^2\OmGW(f)$ & $h_c$ \\%run
\hline
noh1&$   2$&$ 0.3$&$1.9\times10^{-1}$&$ 1.0$&$  16$&$3.83\times10^{-2}$&$3.53\times10^{-4}$&$4.83\times10^{-2}$&$0.78$&$1.09\times10^{-8}$&$4.83\times10^{-14}$\\
noh2&$   6$&$   1$&$6.0\times10^{-2}$&$ 1.0$&$ 4.5$&$3.75\times10^{-2}$&$5.61\times10^{-5}$&$7.06\times10^{-3}$&$0.78$&$1.73\times10^{-9}$&$7.07\times10^{-15}$\\
noh3&$  20$&$   3$&$2.3\times10^{-2}$&$ 1.3$&$ 2.0$&$3.81\times10^{-2}$&$1.11\times10^{-5}$&$1.15\times10^{-3}$&$0.78$&$3.44\times10^{-10}$&$1.15\times10^{-15}$\\
noh4&$  60$&$  10$&$1.0\times10^{-2}$&$ 1.4$&$0.43$&$3.93\times10^{-2}$&$2.62\times10^{-6}$&$1.65\times10^{-4}$&$0.79$&$8.10\times10^{-11}$&$1.65\times10^{-16}$\\
\hline
noh5&$   2$&$ 0.3$&$1.0\times10^{-1}$&  --- &  --- &$1.06\times10^{-2}$&$2.70\times10^{-5}$&$1.40\times10^{-2}$&$0.41$&$8.37\times10^{-10}$&$1.40\times10^{-14}$\\
noh6&$   2$&$ 0.3$&$3.0\times10^{-1}$&  --- &  --- &$9.48\times10^{-2}$&$2.08\times10^{-3}$&$1.02\times10^{-1}$&$ 1.2$&$6.42\times10^{-8}$&$1.02\times10^{-13}$\\
noh7&$   6$&$   1$&$2.0\times10^{-2}$&  --- &  --- &$4.63\times10^{-3}$&$6.56\times10^{-7}$&$8.10\times10^{-4}$&$0.27$&$2.03\times10^{-11}$&$8.11\times10^{-16}$\\
\label{Tsummary512}\end{tabular}}\end{table*}

\begin{table*}[t!]\caption{
Similar to \Tab{Tsummary512}, but for helical turbulence.
}\vspace{12pt}\centerline{\begin{tabular}{ccccccccc|ccc}
Run & $\kf$ & $k_1$ & $f_0$ & $p$ & $\tau$ & $\EEM^{\max}$ & $\EEGW^{\rm sat}$ & $\hrms^{\rm sat}$ &
$B$ [$\mu$G] & $h_0^2\OmGW(f)$ & $h_c$ \\%run
\hline
hel1&$   2$&$ 0.3$&$1.9\times10^{-1}$&$0.67$&$ 100$&$3.90\times10^{-2}$&$4.85\times10^{-4}$&$4.33\times10^{-2}$&$0.79$&$1.50\times10^{-8}$&$4.33\times10^{-14}$\\
hel2&$   6$&$   1$&$5.6\times10^{-2}$&$0.67$&$  20$&$3.81\times10^{-2}$&$5.05\times10^{-5}$&$4.69\times10^{-3}$&$0.78$&$1.56\times10^{-9}$&$4.69\times10^{-15}$\\
hel3&$  20$&$   3$&$2.0\times10^{-2}$&$0.67$&$ 4.0$&$3.96\times10^{-2}$&$7.26\times10^{-6}$&$6.66\times10^{-4}$&$0.80$&$2.24\times10^{-10}$&$6.66\times10^{-16}$\\
hel4&$  60$&$  10$&$6.5\times10^{-3}$&$0.67$&$0.50$&$3.76\times10^{-2}$&$8.15\times10^{-7}$&$7.18\times10^{-5}$&$0.78$&$2.52\times10^{-11}$&$7.18\times10^{-17}$\\
\hline
hel5&$   2$&$ 0.3$&$1.0\times10^{-1}$&  --- &  --- &$1.06\times10^{-2}$&$3.61\times10^{-5}$&$1.08\times10^{-2}$&$0.41$&$1.12\times10^{-9}$&$1.08\times10^{-14}$\\
hel6&$   2$&$ 0.3$&$3.0\times10^{-1}$&  --- &  --- &$9.85\times10^{-2}$&$3.07\times10^{-3}$&$1.12\times10^{-1}$&$ 1.3$&$9.49\times10^{-8}$&$1.12\times10^{-13}$\\
hel7&$   6$&$   1$&$2.0\times10^{-2}$&  --- &  --- &$4.93\times10^{-3}$&$8.33\times10^{-7}$&$6.26\times10^{-4}$&$0.28$&$2.58\times10^{-11}$&$6.26\times10^{-16}$\\
hel8&$   6$&$   1$&$1.0\times10^{-1}$&  --- &  --- &$1.20\times10^{-1}$&$5.09\times10^{-4}$&$1.59\times10^{-2}$&$ 1.4$&$1.57\times10^{-8}$&$1.59\times10^{-14}$\\
\label{Tsummary512hel}\end{tabular}}\end{table*}

In \Tab{Tsummary512}, we have quoted the values of $\EEGW^{\rm sat}$ and
$\hrms^{\rm sat}$ obtained at the end of the simulation.
To compute the relic observable $h_0^2\OmGW$ at the present time,
we have to multiply $\EEGW^{\rm sat}$ by a factor
$(H_\ast/H_0)^2 (a_*/a_0)^4$;
see Refs.~\cite{Pol:2019yex,Pol:2018pao} for details. 
Using $g_*=15$ and $T_*=150\MeV$, we find $H_*=1.8\times10^4\s^{-1}$, 
and thus this factor is $\approx3\times10^{-5}$.
The largest value of $\EEGW^{\rm sat}$ quoted in
\Tab{Tsummary512} is $3.5\times10^{-4}$ and corresponds 
therefore to $h_0^2\OmGW\approx10^{-8}$.
Likewise, the values of $\hrms^{\rm sat}$ in \Tab{Tsummary512} 
have to be multiplied by $a_0^{-1}\approx10^{-12}$ to
obtain the observable $h_c$ at the present time; see \Eq{Dfactor}.
Again, the largest value of $\hrms^{\rm sat}=5\times10^{-2}$ 
corresponds therefore to the observable $h_c=5\times10^{-14}$.

\begin{figure*}[t!]\begin{center}
\includegraphics[width=.99\columnwidth]{pcomp512_oo}
\includegraphics[width=.99\columnwidth]{pcomp512hel_oo}
\end{center}\caption[]{
Similar to \Fig{pcomp512}(a) and (b), but in a double-logarithmic representation
where $\EEM$ is now plotted versus $\Delta t\equiv t-2$, the time after
which the electromagnetic source is turned off.
}\label{pcomp512_oo}\end{figure*}

\begin{table*}[t!]\caption{
Comparison of nonhelical magnetic turbulence (mag) with irrotational (irro)
and vortical (vort) turbulence.
}\vspace{12pt}\centerline{\begin{tabular}{lccccc|ccc}
Type & $f_0$ & $\nu$ & $\EEM^{\max}$ & $\EEGW^{\rm sat}$ & $\hrms^{\rm sat}$ &
$B$ [$\mu$G] & $h_0^2\OmGW(f)$ & $h_c$ \\%run
\hline
%force      eta       EEK         EEGW        hrms       run
magnetic&$1.9\times10^{-1}$&$5.0\times10^{-5}$&$3.83\times10^{-2}$&$3.53\times10^{-4}$&$4.83\times10^{-2}$&$0.78$&$1.09\times10^{-8}$&$4.83\times10^{-14}$\\
vortical&$3.8\times10^{-1}$&$1.0\times10^{-2}$&$4.21\times10^{-2}$&$8.81\times10^{-4}$&$8.26\times10^{-2}$&$0.82$&$2.73\times10^{-8}$&$8.27\times10^{-14}$\\
irrotational&$7.0\times10^{-1}$&$2.0\times10^{-2}$&$4.26\times10^{-2}$&$8.30\times10^{-4}$&$7.95\times10^{-2}$&$0.83$&$2.57\times10^{-8}$&$7.96\times10^{-14}$\\
\label{Tsummary_hydro}\end{tabular}}\end{table*}

To simplify the comparisons, we have arranged 
the forcing amplitude $f_0$ such that
$\EEM^{\max}$ is similar in certain cases.
The values of $\EEM^{\max}$ listed in the upper block of
\Tabs{Tsummary512}{Tsummary512hel} (for nonhelical and helical
hydromagnetic turbulence, respectively) are around 0.038 and
correspond to $0.8\uG$.
The growth phase of $\EEM(t)$ is similar, but the decay is significantly
slower when $\kf$ is smaller.
The GW energy saturates at a value $\EEGW^{\rm sat}$ some time after
$\EEM(t)$ has reached its maximum, and is smaller for larger values
of $\kf$.

It is important to realize that in all four cases, the decay of the 
magnetic energy follows an approximate power law decay, as given
by \Eq{EEMt}.
To see this, we the plot in \Fig{pcomp512_oo} the evolution of $\EEM$
versus $t-2$ in a double-logarithmic representation.
The parameters $p$ and $\tau$ describe the
decay and are also listed in \Tab{Tsummary512}.

\begin{figure*}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pres}
\includegraphics[width=\columnwidth]{EEGW_vs_EEKM}
\end{center}\caption[]{
(a) $\EEGW$ versus $\EEM/\kf$; the straight line shows
$\EEGW=5.2\times10^{-4}\,(\EEM/\kf)^{1/2}$.
(b) Positions of our runs in a diagram showing
$\EEGW^{\rm sat}$ versus $\EEM^{\max}$.
For orientation the old data points of the Ref.~\cite{Pol:2019yex}
are shown as gray symbols.
The open red (filled blue) symbols are for the helical (nonhelical) runs.
The green symbols refer to the two hydromagnetic runs of 
\Tab{Tsummary_hydro}.
}\label{pres}\end{figure*}

\begin{figure*}[t!]\begin{center}
\includegraphics[width=.99\columnwidth]{pspecm}
\includegraphics[width=.99\columnwidth]{pspecm_hel}
\end{center}\caption[]{
$h_0^2\OmGW(f)$ and $h_c(f)$ at the present time for all four runs
presented in \Tab{Tsummary512}, for the nonhelical (left) and helical
(right) runs.
The $2\sigma$ confidence contour for the 30-frequency power law of
the NANOGrav 12.5-year data set is shown in gray.
}\label{pspecm}\end{figure*}

\begin{figure*}[t!]\begin{center}
\includegraphics[width=\textwidth]{pbspecm}
\end{center}\caption[]{
Magnetic energy spectra for the (a) nonhelical and (b) helical cases.
}\label{pbspecm}\end{figure*}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pspecm_other}
\end{center}\caption[]{
Similar to \Fig{pspecm}, but comparing vortical (red) and
irrotational turbulence (blue) with MHD turbulence (black).
}\label{pspecm_other}\end{figure}

Our results confirm that the turbulence decays more slowly for
large values of $\tau$, or small values of $\kf$. 
As already found from earlier simulations \cite{Pol:2019yex},
the GW energy generally decreases with increasing $\kf$.
This is seen more clearly in a diagram of $\EEGW$ versus $\EEM/\kf$;
see \Fig{pres}(a).

For $\kf=2$, we have performed additional simulations with smaller
and larger values $f_0$, both with and without helicity.
The resulting values of $\EEGW^{\rm sat}$ obey quadratic scaling
of the form
\begin{equation}
\EEGW^{\rm sat}=\left(q\EEM^{\max}/\kf\right)^2
\label{qdef}
\end{equation}
with a coefficient $q=1.1$; see the straight line in \Fig{pres}(a).
Only the data point for $\kf=60$ is slightly above the line
represented by \Eq{qdef}.
This could be an artefact of our Reynolds numbers still not being large
enough in our simulations, especially for large value of $\kf$.

To compare with earlier work, we show in \Fig{pres}(b) the positions
of our runs in a $\EEGW^{\rm sat}$ versus $\EEM^{\max}$ diagram.
For orientation, we also show the data points from Ref.~\cite{Pol:2019yex}.
We see that the new data points are well above the older ones of
Ref.~\cite{Pol:2019yex}. 
This is mainly a consequence of using here smaller values of $\kf$
(2--60, compared to 600 in Ref.~\cite{Pol:2019yex}).
For $\kf=2$ and $k_1=0.3$, we show here the results for 
hydrodynamic runs using irrotational and vortical forcings;
see the green symbols in \Fig{pres}(b).
Those runs are listed in \Tab{Tsummary_hydro} and compared with 
the nonhelical magnetic turbulence run `noh1'.

In \Fig{pspecm}, we plot the resulting present-day GW energy and strain
spectra for our four runs with $\kf=2$, 6, 20, and 60, both without and
with helicity in the driving function $\FFF$.
The first two cases with $\kf=2$ and 6 lie well within the frequency
and amplitude range accessible to NANOGrav.
In all cases, the spectra show a sharp drop slightly above the peak
frequency.
This is a consequence of the rapid temporal growth of the spectra,
which leads to a correspondingly large growth at the peak frequency,
while at higher frequencies, the spectrum settled at values that were
determined by somewhat earlier times when the energy was still weaker.

At frequencies below the peak, we now find a spectrum that is even
shallower than the $h_0^2\Omega_{\rm GW}(f)\propto f$ spectrum found
already earlier \cite{Pol:2019yex}.
A spectrum shallower than proportional to $f$, such as the present
$f^{1/2}$ spectrum, could perhaps be explained by the finite size of the
computational domain; see Ref.~\cite{Brandenburg:2021aln}, who found even
a $f^{-1/2}$ spectrum for $\Omega_{\rm GW}(f)$. 
Alternatively, the shallower spectrum might well be physical,
or at least significantly extended over a substantial frequency
interval below the peak frequency, for example due to inverse
%cascading in helical and nonhelical cases.
%TK should not we indicate in the bracket or footnote that non-helical case -> non-helical inverse transfer?. For non-specialists there is inverse cascade in helical cases only
%AB: could just give references, like so:
cascading in helical \cite{FPLM75} and nonhelical \cite{BK17} cases.

In the absence of sources, a $\Omega_{\rm GW}(f)\propto f^\alpha$ spectrum
implies $h_c(f)\propto f^{\alpha/2-1}$ for arbitrary spectral indices $\alpha$.
For $\alpha=1/2$, we would thus expect $h_c(f)\propto f^{-3/4}$.
However, the observed strain spectrum, $h_c(f)\propto f^{-1/2}$, seems to agree
with that found previously from numerical simulations \cite{Pol:2019yex}. 
However, looking more carefully at the strain spectrum for $\kf=60$, we see
a $h_c(f)\propto f^{-3/4}$ spectrum is actually compatible with the simulation;
see the corresponding dashed-dotted line in \Fig{pspecm}(c).
This agreement is probably related to the fact that the turnover time is
shorter for the run with $\kf=60$, compared with those at smaller values 
(i.e., longer turbulence driving time will allow for
more efficient inverse cascading).

In the runs with helicity, we do find $h_0^2\Omega_{\rm GW}(f)\propto f$,
together with a slight enhancement just before reaching a maximum.
The subsequent decay for larger values of $f$ is much steeper in the
case with helicity than without.
Furthermore, in $h_c(f)$ we see a sharper drop to the right of the 
maximum than in simulations without helicity.
These differences in the spectra for helical and nonhelical cases
are surprisingly strong and might allow us to infer the presence
of magnetic helicity once such a spectrum is detected. 

It is important to note that the $h_0^2\Omega_{\rm GW}\propto f$ spectra
in Figs.~\ref{pres}(a) and (c) show an increase towards smaller $\kf$.
This is to be expected from \Eq{qdef}, but it was not included in the
sketch of Ref.~\cite{Neronov:2020qrl}; see their Fig.~(1).
By contrast, in their Eq.~(4), an effectively cubic dependence on
the magnetic energy was motivated.

The underlying magnetic energy spectrum is shown in \Fig{pbspecm}(a) for
nonhelical and in \Fig{pbspecm}(b) helical cases where $\kf$ is ranging from 2 to 60.
Those are averaged spectra obtained by averaging over the time interval
$15\leq t\leq20$.
In the nonhelical case, the amplitude of the spectrum is smaller for
larger values of $\kf$, because here the energy has decayed more rapidly.
In the helical case, the spectra have approximately the same height for
all values of $\kf$.
This is because the height of the spectrum
is related to the helicity, which is conserved.
For small values of $\kf$, the spectrum has a more extended subinertial
range.
This is because the turnover time is larger and there was not enough
time for the inverse cascade to produce energy and small values of $k$.

Finally, we compare the results for two types of purely hydrodynamic
turbulence with vortical and irrotational forcings of \Tab{Tsummary_hydro}.
The result is shown in \Fig{pspecm_other}.
All these cases are for plane wave forcings.
For irrotational forcing, we do not see the sharp drop-off of spectral
power for frequencies above the peak value as in the vortical case. 
This suggests that in the inertial range of irrotational turbulence,
there is still some power to contribute to GW driving compared with the
vortical case, where this is almost not possible at all. 
However, the spectrum in the irrotational case shows a fairly steep
spectrum proportional to $f^{-7}$, so the effect on GW production is
here also rather weak. 
Nevertheless, the spectral form of the peak might give interesting
diagnostic clues about the nature of turbulent driving at the time of
GW production.

We mention in passing that in earlier work, it was found that
irrotational turbulence is much more efficient in driving GWs
than vortical turbulence \cite{Pol:2019yex}.
Remarkably, here this is no longer the case and vortical and
irrotational turbulence have rather similar GW energies.
This could be related to the small value of $\kf$, possibly
combined with a comparatively short time of driving.
However, to clarify this further, more targeted numerical
experiments would need to be performed.

%EC putting this plot here for now:
\begin{figure*}\begin{center}
%\includegraphics[scale=1]{fig/sourcesgrid.eps}
\includegraphics[scale=1]{fig/sourcesGrid2.eps}
\end{center}
\caption{
%EC: need to modify caption, but for now:
Upper left: NANOGrav 12.5-year contours, gray shaded region corresponds to possible peak frequencies of $\OmGW$ for source described in this paper, for 1-6 bubbles at QCD; 
Upper Right: cosmic strings: four average power spectrum models described in \cite{Blanco-Pillado:2021ygr} are mono (orange), kink (green), cusp (blue), and a spectrum computed from a simulated gravitational backreaction model (`BOS') (magenta); 
%TK BOS has to be defined 
Lower Left: Inflation \cite{Vagnozzi:2020gtf}
$n_T-r$ parameter space corresponding {\it within} 
%TK Maybe for each source indicate parameter range, at least for inflation
parameter space of NANOGrav 12.5-year results;
Lower Right: GWs from SMBHBs are expected to have a spectral index of $\gamma_{\rm CP} = 13/3$, which corresponds to $A_{\rm CP}=[1.4, 2.7]\times10^{-15}$ from the NANOGrav 5-frequency and broken power laws.
}
%\label{fig:}
\end{figure*}

%EC: gray shaded 1-6 bubbles at QCD: assumes frequency of \OmGW^{sat} is located at f_GW = 2 f_B, where f_B corresponds to \OmM^max; [3.57384\times10^{-8}, 2.14431\times10^{-7}] Hz
%AB: we could also take upper and lower lines from the simulations. For kf=2, we have noh5+6 and hel5+6, but I can also try to produce kf=6 runs in that range.
%EC: figure updated with runs you sent me
%EC: Upper Left: shaded regions between noh7+8 (purple), noh5+6(yellow), hel5+6 (magenta) (?--need to verify which runs file names correspond to)
%TK is this finalized?
%AB: do we then still need the gray box?
%EC I can remove the gray. should I then plot the runs for the entire frequency range or just the range within the gray region?

%EC is it of interest to comment on how scaling of \OmGW(f) compares to NANOGrav-12.5yr results? 
%eg. ~f^{1/2} scaling implies \gamma = 4.5, which is at the edge of the 68% contours (5-frequency and broken power laws); ~f^1 implies \gamma = 4, which is within the 95% contours for 5 frequency and broken power laws, at edge of 95% contour for (2021) 30-frequency power law
%TK Yes, please draft and I will add more - the important question in the ideal case can we separate QCD and other sources mentioned here?
%EC: should I draft a paragraph on this point alone, or should we note this where we give the results of how \OmGW scales with f?
%EC: is there a paragraph drafted related to this figure? (should I draft something?)
%AB: drafting something on this sounds like a good idea.

%EC

\section{Conclusions}
\label{Concl}

In the present work, we have shown that the magnetic stress from MHD
turbulence with scales comparable to the cosmological horizon scale at the
time of the QCD phase transition can drive GWs in the range accessible
to NANOGrav, if the magnetic energy density is 3--10\% of the radiation
energy density.
The low-frequency tail below the peak frequency at $10\nHz$ or so is
shallower in the nonhelical case than in the helical one, i.e.,
$\propto f^{1/2}$ compared to $\propto f$.
Both scaling are, however, shallower than what was expected based on
earlier analytical calculations.
Also the inertial range spectrum above the peak is shallower without
helicity than with, but here, both spectra are steeper than what
is expected if the GW spectrum was a direct consequence of the MHD
turbulence spectrum.
The reason for this is primarily the relatively short time of turbulent
%driving, which we chose to be 100 Hubble times.
%TK confused to be 100 Hubble times
%AB: yes, modified
driving (one Hubble time).
This is short compared with the turnover time which,
for our runs with the smallest $\kf$ of two, is much longer:
16 (100) Hubble times for our runs without (with) helicitly.
Therefore, there was not enough time to fully establish the GW spectrum
at high wave numbers.
%AB.
%For our runs with larger values of $\kf$, this effect is less
%AB: modified
For our earlier runs with larger values of $\kf$, this effect was less
pronounced than for smaller values of $\kf$, but it is still quite noticeable,
especially in the helical case where forward cascading is weaker than
in the nonhelical case. 

Our work has led to new insights regarding the possibility of using
an observed GW spectrum for making statements about the nature of the
underlying turbulence in the early universe.
One is the already mentioned slope of the subinertial range spectrum.
Another is the position of the peak of the spectrum.
Finally, there is the strength of the drop of the spectral power for
frequencies above the peak frequency, and the subsequent slope after
the drop, which is most likely too small to be detectable.
This, however, depends on the duration of turbulent driving and could
be higher if the driving time was longer. 
The specific features of the spectrum near the peak are different
for helical and nonhelical turbulence.
This could, in principle, give information about the presence 
of parity violation, when would also lead to circularly polarized
GWs. 

\vspace{2mm}
Data availability---The source code used for the
simulations of this study, the {\sc Pencil Code},
is freely available from Ref.~\cite{PC}.
The simulation setups and the corresponding data
are freely available from Ref.~XXX.


\acknowledgements
%\vspace{3mm}
Support through the Swedish Research Council, grant 2019-04234,
and Shota Rustaveli GNSF  (grant FR/18-1462)
are gratefully acknowledged.
We acknowledge the allocation of computing resources provided by the
Swedish National Allocations Committee at the Center for Parallel
Computers at the Royal Institute of Technology in Stockholm.

%r e f
\begin{thebibliography}{99}

\bibitem{Caprini:2018mtu}
C.~Caprini and D.~G.~Figueroa,
``Cosmological Backgrounds of Gravitational Waves,''
Class. Quant. Grav. \textbf{35}, 163001 (2018) 
%doi:10.1088/1361-6382/aac608
[arXiv:1801.04268 [astro-ph.CO]].

\bibitem{Romano:2016dpx}
J.~D.~Romano and N.~J.~Cornish,
``Detection methods for stochastic gravitational-wave backgrounds: a unified treatment,''
Living Rev. Rel. \textbf{20},  2 (2017)
%doi:10.1007/s41114-017-0004-1
[arXiv:1608.06889 [gr-qc]].

\bibitem{Arzoumanian:2020vkk}
Z.~Arzoumanian \textit{et al.} [NANOGrav],
``The NANOGrav 12.5~yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background,''
Astrophys. J. Lett. \textbf{905},  L34 (2020) 
%doi:10.3847/2041-8213/abd401
[arXiv:2009.04496 [astro-ph.HE]].


\bibitem{Vagnozzi:2020gtf}
S.~Vagnozzi,
``Implications of the NANOGrav results for inflation,''
Mon. Not. Roy. Astron. Soc. \textbf{502}, L11 (2021)
doi:10.1093/mnrasl/slaa203
[arXiv:2009.13432 [astro-ph.CO]].

\bibitem{Zhou:2020kkf}
Z.~Zhou, J.~Jiang, Y.~F.~Cai, M.~Sasaki and S.~Pi,
``Primordial black holes and gravitational waves from resonant amplification during inflation,''
Phys. Rev. D \textbf{102},  103527 (2020)
%doi:10.1103/PhysRevD.102.103527
[arXiv:2010.03537 [astro-ph.CO]].

\bibitem{Tahara:2020fmn}
H.~W.~H.~Tahara and T.~Kobayashi,
``Nanohertz gravitational waves from NEC violation in the early universe,''
Phys. Rev. D \textbf{102}, 123533 (2020)
%doi:10.1103/PhysRevD.102.123533
[arXiv:2011.01605 [gr-qc]].

\bibitem{Kuroyanagi:2020sfw}
S.~Kuroyanagi, T.~Takahashi and S.~Yokoyama,
``Blue-tilted inflationary tensor spectrum and reheating in the light of NANOGrav results,''
[arXiv:2011.03323 [astro-ph.CO]].

\bibitem{Cai:2020qpu}
Y.~Cai and Y.~S.~Piao,
``Intermittent NEC violations during inflation and primordial gravitational waves,''
[arXiv:2012.11304 [gr-qc]].

\bibitem{Ellis:2020ena}
J.~Ellis and M.~Lewicki,
``Cosmic String Interpretation of NANOGrav Pulsar Timing Data,''
[arXiv:2009.06555 [astro-ph.CO]].

\bibitem{Blasi:2020mfx}
S.~Blasi, V.~Brdar and K.~Schmitz,
``Has NANOGrav found first evidence for cosmic strings?,''
Phys. Rev. Lett. \textbf{126}, 041305 (2021) 
%doi:10.1103/PhysRevLett.126.041305
[arXiv:2009.06607 [astro-ph.CO]].

\bibitem{Buchmuller:2020lbh}
W.~B\"uchm\"uller, V.~Domcke and K.~Schmitz,
``From NANOGrav to LIGO with metastable cosmic strings,''
Phys. Lett. B \textbf{811}, 135914 (2020)
%doi:10.1016/j.physletb.2020.135914
[arXiv:2009.10649 [astro-ph.CO]].

\bibitem{Samanta:2020cdk}
R.~Samanta and S.~Datta,
``Gravitational wave complementarity and impact of NANOGrav data on gravitational leptogenesis: cosmic strings,''
[arXiv:2009.13452 [hep-ph]].

\bibitem{Liu:2020mru}
J.~Liu, R.~G.~Cai and Z.~K.~Guo,
``Large anisotropies of the stochastic gravitational wave background from cosmic domain walls,''
[arXiv:2010.03225 [astro-ph.CO]].

\bibitem{Paul:2020wbz}
A.~Paul, U.~Mukhopadhyay and D.~Majumdar,
``Gravitational Wave Signatures from Domain Wall and Strong First-Order Phase Transitions in a Two Complex Scalar extension of the Standard Model,''
[arXiv:2010.03439 [hep-ph]].

\bibitem{Lewicki:2020azd}
M.~Lewicki and V.~Vaskonen,
``Gravitational waves from colliding vacuum bubbles in gauge theories,''
[arXiv:2012.07826 [astro-ph.CO]].

\bibitem{Nakai:2020oit}
Y.~Nakai, M.~Suzuki, F.~Takahashi and M.~Yamada,
``Gravitational Waves and Dark Radiation from Dark Phase Transition: Connecting NANOGrav Pulsar Timing Data and Hubble Tension,''
[arXiv:2009.09754 [astro-ph.CO]].

\bibitem{Addazi:2020zcj}
A.~Addazi, Y.~F.~Cai, Q.~Gan, A.~Marciano and K.~Zeng,
``NANOGrav results and Dark First Order Phase Transitions,''
[arXiv:2009.10327 [hep-ph]].

\bibitem{Kitajima:2020rpm}
N.~Kitajima, J.~Soda and Y.~Urakawa,
``Nano-Hz gravitational wave signature from axion dark matter,''
[arXiv:2010.10990 [astro-ph.CO]].

\bibitem{Ramberg:2020oct}
N.~Ramberg and L.~Visinelli,
``The QCD Axion and Gravitational Waves in light of NANOGrav results,''
[arXiv:2012.06882 [astro-ph.CO]].

\bibitem{Lee:2020wfn}
V.~S.~H.~Lee, A.~Mitridate, T.~Trickle and K.~M.~Zurek,
``Probing Small-Scale Power Spectra with Pulsar Timing Arrays,''
[arXiv:2012.09857 [astro-ph.CO]].

\bibitem{Gorghetto:2021fsn}
M.~Gorghetto, E.~Hardy and H.~Nicolaescu,
``Observing Invisible Axions with Gravitational Waves,''
[arXiv:2101.11007 [hep-ph]].

\bibitem{Neronov:2020qrl}
A.~Neronov, A.~Roper Pol, C.~Caprini and D.~Semikoz,
``NANOGrav signal from MHD turbulence at QCD phase transition in the early universe,''
[arXiv:2009.14174 [astro-ph.CO]].

\bibitem{Abe:2020sqb}
K.~T.~Abe, Y.~Tada and I.~Ueda,
``Induced gravitational waves as a cosmological probe of the sound speed during the QCD phase transition,''
[arXiv:2010.06193 [astro-ph.CO]].

\bibitem{Quashnock:1988vs}
J.~M.~Quashnock, A.~Loeb and D.~N.~Spergel,
``Magnetic Field Generation During the Cosmological QCD Phase Transition,''
Astrophys. J. Lett. \textbf{344}, L49 (1989)
%doi:10.1086/185528

\bibitem{Sigl:1996dm}
G.~Sigl, A.~V.~Olinto and K.~Jedamzik,
``Primordial magnetic fields from cosmological first order phase transitions,''
Phys. Rev. D \textbf{55}, 4582 (1997)
%doi:10.1103/PhysRevD.55.4582
[arXiv:astro-ph/9610201 [astro-ph]].

\bibitem{Forbes:2000gr}
M.~M.~Forbes and A.~R.~Zhitnitsky,
``Primordial galactic magnetic fields from domain walls at the QCD phase transition,''
Phys. Rev. Lett. \textbf{85}, 5268 (2000) 
%doi:10.1103/PhysRevLett.85.5268
[arXiv:hep-ph/0004051 [hep-ph]].

\bibitem{Boyanovsky:2002wa}
D.~Boyanovsky, H.~J.~de Vega and M.~Simionato,
``Large scale magnetogenesis from a nonequilibrium phase transition in the radiation dominated era,''
Phys. Rev. D \textbf{67}, 123505 (2003)
%doi:10.1103/PhysRevD.67.123505
[arXiv:hep-ph/0211022 [hep-ph]].




\bibitem{Kisslinger:2005uy}
L.~S.~Kisslinger, S.~Walawalkar and M.~B.~Johnson,
``Basic treatment of QCD phase transition bubble nucleation,''
Phys. Rev. D \textbf{71}, 065017 (2005)
%doi:10.1103/PhysRevD.71.065017
[arXiv:hep-ph/0503157 [hep-ph]].

\bibitem{Buividovich:2009wi}
P.~V.~Buividovich, M.~N.~Chernodub, E.~V.~Luschevskaya and M.~I.~Polikarpov,
``Numerical evidence of chiral magnetic effect in lattice gauge theory,''
Phys. Rev. D \textbf{80}, 054503 (2009)
%doi:10.1103/PhysRevD.80.054503
[arXiv:0907.0494 [hep-lat]].

\bibitem{Urban:2009sw}
F.~R.~Urban and A.~R.~Zhitnitsky,
``Large-Scale Magnetic Fields, Dark Energy and QCD,''
Phys. Rev. D \textbf{82}, 043524 (2010)
%doi:10.1103/PhysRevD.82.043524
[arXiv:0912.3248 [astro-ph.CO]].

\bibitem{Rischke:2003mt}
  D.~H.~Rischke,
  ``The Quark gluon plasma in equilibrium,''
  Prog.\ Part.\ Nucl.\ Phys.\  {\bf 52}, 197 (2004) 
  %doi:10.1016/j.ppnp.2003.09.002
  [arXiv: nucl-th/0305030].
  
  \bibitem{Kiskis:2003rd}
J.~Kiskis, R.~Narayanan and H.~Neuberger,
``Does the crossover from perturbative to nonperturbative physics in QCD become a phase transition at infinite N?,''
Phys. Lett. B \textbf{574}, 65 (2003)
%doi:10.1016/j.physletb.2003.08.070
[arXiv:hep-lat/0308033 [hep-lat]].

\bibitem{Aoki:2006we}
  Y.~Aoki, G.~Endrodi, Z.~Fodor, S.~D.~Katz and K.~K.~Szabo,
  ``The Order of the quantum chromodynamics transition predicted by the standard model of particle physics,''
  Nature {\bf 443}, 675 (2006) 
  %doi:10.1038/nature05120
 [arXiv:hep-lat/0611014].

\bibitem{Blaizot:2013hx}
J.~P.~Blaizot, E.~Iancu and Y.~Mehtar-Tani,
``Medium-induced QCD cascade: democratic branching and wave turbulence,''
Phys. Rev. Lett. \textbf{111}, 052001 (2013) 
%doi:10.1103/PhysRevLett.111.052001
[arXiv:1301.6102 [hep-ph]].

\bibitem{Miniati:2017kah}
  F.~Miniati, G.~Gregori, B.~Reville and S.~Sarkar,
  ``Axion-Driven Cosmic Magnetogenesis during the QCD Crossover,''
  Phys.\ Rev.\ Lett.\  {\bf 121},  021301 (2018) 
%  doi:10.1103/PhysRevLett.121.021301
  [arXiv:1708.07614 [astro-ph.CO]].

\bibitem{Witten:1984rs}
E.~Witten,
``Cosmic Separation of Phases,''
Phys. Rev. D \textbf{30}, 272 (1984). 
%doi:10.1103/PhysRevD.30.272

\bibitem{Applegate:1985qt}
J.~H.~Applegate and C.~J.~Hogan,
``Relics of Cosmic Quark Condensation,''
Phys. Rev. D \textbf{31}, 3037 (1985).
%doi:10.1103/PhysRevD.31.3037

\bibitem{Midorikawa:1985gy}
S.~Midorikawa,
``Bubble Collisions in the Cosmological Quark - Hadron Phase Transition,''
Phys. Lett. B \textbf{158}, 107 (1985). 
%doi:10.1016/0370-2693(85)91373-5

\bibitem{Hindmarsh:1991ay}
M.~Hindmarsh,
``Axions and the QCD phase transition,''
Phys. Rev. D \textbf{45}, 1130 (1992).
%doi:10.1103/PhysRevD.45.1130

\bibitem{Cottingham:1993rv}
W.~N.~Cottingham, D.~Kalafatis and R.~Vinh Mau,
``Bubble nucleation rates in first order phase transitions,''
Phys. Rev. B \textbf{48}, 6788(1993).
%doi:10.1103/PhysRevB.48.6788

\bibitem{Ignatius:1993qn}
J.~Ignatius, K.~Kajantie, H.~Kurki-Suonio and M.~Laine,
``The growth of bubbles in cosmological phase transitions,''
Phys. Rev. D \textbf{49}, 3854 (1994)
%doi:10.1103/PhysRevD.49.3854
[arXiv:astro-ph/9309059 [astro-ph]].

\bibitem{Ignatius:1997va}
J.~Ignatius,
``Early stages of growth of QCD and electroweak bubbles,''
[arXiv:hep-ph/9708383 [hep-ph]].

\bibitem{Strumia:1999fv}
A.~Strumia and N.~Tetradis,
``Bubble nucleation rates for cosmological phase transitions,''
JHEP \textbf{11}, 023 (1999)
%doi:10.1088/1126-6708/1999/11/023
[arXiv:hep-ph/9904357 [hep-ph]].

\bibitem{Strumia:1998nf}
A.~Strumia and N.~Tetradis,
``A Consistent calculation of bubble nucleation rates,''
Nucl. Phys. B \textbf{542}, 719 (1999)
%doi:10.1016/S0550-3213(98)00804-9
[arXiv:hep-ph/9806453 [hep-ph]].

\bibitem{Schwarz:2003du}
D.~J.~Schwarz,
"The first second of the universe,''
Annalen Phys. \textbf{12}, 220 (2003)
%doi:10.1002/andp.200310010
[arXiv:astro-ph/0303574 [astro-ph]].

\bibitem{Tawfik:2011mw}
A.~Tawfik,
``The Hubble parameter in the early universe with viscous QCD matter and finite cosmological constant,''
Annalen Phys. \textbf{523}, 423 (2011)
%doi:10.1002/andp.201100038
[arXiv:1102.2626 [gr-qc]].

\bibitem{Tawfik:2011sh}
A.~Tawfik and T.~Harko,
``Quark-Hadron Phase Transitions in Viscous Early Universe,''
Phys. Rev. D \textbf{85}, 084032 (2012)
%doi:10.1103/PhysRevD.85.084032
[arXiv:1108.5697 [astro-ph.CO]].

\bibitem{MaggioreText}
    M. Maggiore, "Gravitational Waves. Vol. 1: Theory and Experiments" 
    (Oxford University Press, 2007).

\bibitem{K&T}
    E. W. Kolb and M. S. Turner, ``The Early Universe," Addison-Wesley Publishing Company (1990).

\bibitem{Husdal:2016haj}
L.~Husdal,
``On Effective Degrees of Freedom in the Early Universe,''
Galaxies \textbf{4},  78 (2016) 
%doi:10.3390/galaxies4040078
[arXiv:1609.04979 [astro-ph.CO]].   

\bibitem{Burke-Spolaor:2018bvk}
S.~Burke-Spolaor, S.~R.~Taylor, M.~Charisi, T.~Dolch, J.~S.~Hazboun, A.~M.~Holgado, L.~Z.~Kelley, T.~J.~W.~Lazio, D.~R.~Madison and N.~McMann, \textit{et al.}
``The Astrophysics of Nanohertz Gravitational Waves,''
Astron. Astrophys. Rev. \textbf{27}, 5 (2019)
%doi:10.1007/s00159-019-0115-7
[arXiv:1811.08826 [astro-ph.HE]].

%\cite{Sesana:2008mz}
\bibitem{Sesana:2008mz}
A.~Sesana, A.~Vecchio and C.~N.~Colacino,
``The stochastic gravitational-wave background from massive black hole binary systems: implications for observations with Pulsar Timing Arrays,''
Mon. Not. Roy. Astron. Soc. \textbf{390}, 192 (2008)
%doi:10.1111/j.1365-2966.2008.13682.x
[arXiv:0804.4476 [astro-ph]].

\bibitem{Sanidas:2012ee}
S.~A.~Sanidas, R.~A.~Battye and B.~W.~Stappers,
``Constraints on cosmic string tension imposed by the limit on the stochastic gravitational wave background from the European Pulsar Timing Array,''
Phys. Rev. D \textbf{85}, 122003 (2012)
%doi:10.1103/PhysRevD.85.122003
[arXiv:1201.2419 [astro-ph.CO]].

\bibitem{Cutler:2013aja}
C.~Cutler, S.~Burke-Spolaor, M.~Vallisneri, J.~Lazio and W.~Majid,
``The Gravitational-Wave Discovery Space of Pulsar Timing Arrays,''
Phys. Rev. D \textbf{89}, 042003 (2014)
%doi:10.1103/PhysRevD.89.042003
[arXiv:1309.2581 [gr-qc]].

\bibitem{Blanco-Pillado:2017rnf}
J.~J.~Blanco-Pillado, K.~D.~Olum and X.~Siemens,
``New limits on cosmic strings from gravitational wave observation,''
Phys. Lett. B \textbf{778}, 392 (2018)
%doi:10.1016/j.physletb.2018.01.050
[arXiv:1709.02434 [astro-ph.CO]].

\bibitem{Campeti:2020xwn}
P.~Campeti, E.~Komatsu, D.~Poletti and C.~Baccigalupi,
``Measuring the spectrum of primordial gravitational waves with CMB, PTA and Laser Interferometers,''
JCAP \textbf{01}, 012 (2021)
%doi:10.1088/1475-7516/2021/01/012
[arXiv:2007.04241 [astro-ph.CO]].

\bibitem{Grishchuk:1974ny}
L.~P.~Grishchuk,
``Amplification of gravitational waves in an istropic universe,''
Sov. Phys. JETP \textbf{40}, 409 (1975).

\bibitem{Rubakov:1982df}
V.~A.~Rubakov, M.~V.~Sazhin and A.~V.~Veryaskin,
``Graviton Creation in the Inflationary Universe and the Grand Unification Scale,''
Phys. Lett. B \textbf{115}, 189 (1982). 
%doi:10.1016/0370-2693(82)90641-4

\bibitem{Starobinsky:1979ty}
A.~A.~Starobinsky,
``Spectrum of relict gravitational radiation and the early state of the universe,''
JETP Lett. \textbf{30}, 682(1979). 

\bibitem{Caprini:2010xv}
C.~Caprini, R.~Durrer and X.~Siemens,
``Detection of gravitational waves from the QCD phase transition with pulsar timing arrays,''
Phys. Rev. D \textbf{82}, 063511 (2010)
%doi:10.1103/PhysRevD.82.063511
[arXiv:1007.1218 [astro-ph.CO]].

\bibitem{Capozziello:2018qjs}
S.~Capozziello, M.~Khodadi and G.~Lambiase,
``The quark chemical potential of QCD phase transition and the stochastic background of gravitational waves,''
Phys. Lett. B \textbf{789}, 626 (2019)
%doi:10.1016/j.physletb.2019.01.004
[arXiv:1808.06188 [gr-qc]].

\bibitem{Kobakhidze:2017mru}
A.~Kobakhidze, C.~Lagger, A.~Manning and J.~Yue,
``Gravitational waves from a supercooled electroweak phase transition and their detection with pulsar timing arrays,''
Eur. Phys. J. C \textbf{77},  570 (2017)
%doi:10.1140/epjc/s10052-017-5132-y
[arXiv:1703.06552 [hep-ph]].

\bibitem{Ruderman:1975ju}
M.~A.~Ruderman and P.~G.~Sutherland,
``Theory of pulsars: Polar caps, sparks, and coherent microwave radiation,''
Astrophys. J. \textbf{196}, 51 (1975).
%doi:10.1086/153393

\bibitem{Gold:1968}
T. Gold, 
"Rotating Neutron Stars as the Origin of the Pulsating Radio Sources", Nature, {\bf 218}, 731 (1968).

\bibitem{MaggioreText2}
    M. Maggiore, "Gravitational Waves. Vol. 2: "Astrophysics and Cosmology" (Oxford University Press, 2017).

\bibitem{Goncharov:2020krd}
B.~Goncharov, D.~J.~Reardon, R.~M.~Shannon, X.~J.~Zhu, E.~Thrane, M.~Bailes, N.~D.~R.~Bhat, S.~Dai, G.~Hobbs and M.~Kerr, \textit{et al.}
``Identifying and mitigating noise sources in precision pulsar timing data sets,''
%doi:10.1093/mnras/staa3411
[arXiv:2010.06109 [astro-ph.HE]].

\bibitem{PTA} http:$//$ipta4gw.org$/$

\bibitem{Hobbs:2009yy}
G.~Hobbs, A.~Archibald, Z.~Arzoumanian, D.~Backer, M.~Bailes, N.~D.~R.~Bhat, M.~Burgay, S.~Burke-Spolaor, D.~Champion and I.~Cognard, \textit{et al.}
``The international pulsar timing array project: using pulsars as a gravitational wave detector,''
Class. Quant. Grav. \textbf{27}, 084013 (2010)
%doi:10.1088/0264-9381/27/8/084013
[arXiv:0911.5206 [astro-ph.SR]].

\bibitem{Brazier:2019mmu}
A.~Brazier, S.~Chatterjee, T.~Cohen, J.~M.~Cordes, M.~E.~DeCesar, P.~B.~Demorest, J.~S.~Hazboun, M.~T.~Lam, R.~S.~Lynch and M.~A.~McLaughlin, \textit{et al.}
``The NANOGrav Program for Gravitational Waves and Fundamental Physics,''
[arXiv:1908.05356 [astro-ph.IM]].
    
\bibitem{Brandenburg:1996fc}
  A.~Brandenburg, K.~Enqvist and P.~Olesen,
  ``Large scale magnetic fields from hydromagnetic turbulence in the very early universe,''
  Phys.\ Rev.\ D {\bf 54}, 1291 (1996)
%doi:10.1103/PhysRevD.54.1291
[arXiv:astro-ph/9602031 [astro-ph]]. 

\bibitem{Brandenburg:2017neh}
  A.~Brandenburg, T.~Kahniashvili, S.~Mandal, A.~R.~Pol, A.~G.~Tevzadze and T.~Vachaspati,
  ``Evolution of hydromagnetic turbulence from the electroweak phase transition,''
  Phys.\ Rev.\ D {\bf 96}, 123528 (2017).
  %doi:10.1103/PhysRevD.96.123528
[arXiv:1711.03804 [astro-ph.CO]].

\bibitem{Pol:2018pao}
A.~Roper Pol, A.~Brandenburg, T.~Kahniashvili, A.~Kosowsky and S.~Mandal,
``The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence,''
Geophys. Astrophys. Fluid Dynamics \textbf{114}, 130 (2020)
%doi:10.1080/03091929.2019.1653460
[arXiv:1807.05479 [physics.flu-dyn]].

\bibitem{Pol:2019yex}
  A.~Roper~Pol, S.~Mandal, A.~Brandenburg, T.~Kahniashvili and A.~Kosowsky,
  ``Numerical Simulations of Gravitational Waves from Early-Universe Turbulence,''
  Phys. Rev. D {\bf 102}, 083512 (2020) 
  %doi:10.1103/PhysRevD.102.083512 
[arXiv:1903.08585 [astro-ph.CO]].

\bibitem{BP56}
G. K. Batchelor and I. Proudman\yptrsa{1956}{248}{369} 
{The large-scale structure of homogeneous turbulence}

\bibitem{Saf67}
P. G. Saffman\ypf{1967}{10}{1349} 
{Note on decay of homogeneous turbulence}

\bibitem{Brandenburg:2017rnt}
A.~Brandenburg, T.~Kahniashvili, S.~Mandal, A.~R.~Pol, A.~G.~Tevzadze and T.~Vachaspati,
``The dynamo effect in decaying helical turbulence,''
Phys. Rev. Fluids. \textbf{4}, 024608 (2019)
%doi:10.1103/PhysRevFluids.4.024608
[arXiv:1710.01628 [physics.flu-dyn]].

\bibitem{Hosking21}
D. N. Hosking and A. A. Schekochihin\arXiv{2021}
{2012.01393}
{Reconnection-controlled decay of magnetohydrodynamic turbulence and the role of invariants}

\bibitem{PC} 
A. Brandenburg, A. Johansen, P. A. Bourdin, W. Dobler, W. Lyra, M. Rheinhardt, S. Bingert, N. E. L. Haugen, A. Mee, F. Gent, N. Babkovskaia, C.-C. Yang, T. Heinemann, B. Dintrans, D. Mitra, S. Candelaresi, J. Warnecke, P. J. K\"apyl\"a, A. Schreiber, P. Chatterjee, M. J. K\"apyl\"a, X.-Y. Li, J. Kr\"uger, J. R. Aarnes, G. R. Sarson, J. S. Oishi, J. Schober, R. Plasson, C. Sandin, E. Karchniwy, L. F. S. Rodrigues, A. Hubbard, G. Guerrero, A. Snodin, I. R. Losada, J. Pekkil\"a, and C. Qian\yjour{2021}{Journal of Open Source Software}{6}{2807}
{The Pencil Code, a modular MPI code for partial differential equations and particles: multipurpose and multiuser-maintained}

\bibitem{Brandenburg:2021aln}
A. Brandenburg, Y. He, T. Kahniashvili, M. Rheinhardt, and J. Schober\papj{2021}
{2101.08178}
{Gravitational waves from the chiral magnetic effect}

\bibitem{FPLM75}
U. Frisch, A. Pouquet, J. L\'eorat, and A. Mazure\yjfm{1975}{68}{769}
{778}{Possibility of an inverse cascade of magnetic helicity in hydrodynamic turbulence}

\bibitem{BK17}
A. Brandenburg and T. Kahniashvili\yprl{2017}{118}{055102}
{Classes of hydrodynamic and magnetohydrodynamic turbulent decay}

\bibitem{Blanco-Pillado:2021ygr}
J.~J.~Blanco-Pillado, K.~D.~Olum and J.~M.~Wachter,
``Comparison of cosmic string and superstring models to NANOGrav 12.5-year results,''
[arXiv:2102.08194 [astro-ph.CO]].

\end{thebibliography}
\end{document}
