\documentclass[prd,twocolumn,aps,psfig,showpacs,nofootinbib,nobibnotes,superscriptaddress,preprintnumbers,times]{revtex4}
\setlength{\topmargin}{-14mm}
\usepackage{graphicx,bm,color}
\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 defs
\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{\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\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{\GeV}{\,{\rm GeV}}
\newcommand{\T}{\,{\rm T}}
\newcommand{\G}{\,{\rm G}}
\newcommand{\Hz}{\,{\rm Hz}}
\newcommand{\mHz}{\,{\rm mHz}}
\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}[2]{, 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}[3]{, Phys.\ Rev.\ D {\bf #2}, #3 (#1).}
\newcommand{\ypre}[3]{, Phys.\ Rev.\ E {\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}[3]{, 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}[2]{, Astrophys. J., submitted, 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}[3]{, 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}[4]{, #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}
\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}
%
%TK for now the authors are listed alphabetically

\begin{abstract}
  We perform the numerical simulations of gravitational waves (GWs) induced by the turbulent sources possibly present at cosmological quark chromodynamic (QCD) phase transitions. 
  
\end{abstract}

\maketitle

\section{Introduction}
Gravitational waves astronomy opens a new window to study the physical
processes in the early universe.
Relic gravitational waves 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, see for a review
\cite{Caprini:2018mtu} and references therein.
On the other hand the detection of these relic gravitational waves
is a challenging task due to their small amplitudes, the specific
range of the characteristic frequencies, and astrophysical foregrounds
\cite{Romano:2016dpx}.
Despite of tremendous advancements in gravitational
waves detection techniques, the stochastic gravitational waves background
of cosmological origin remained unobserved yet.

Recently NanoGrav collaboration reported the strong evidence of the stochastic gravitational waves background \cite{Arzoumanian:2020vkk} that was interpreted as a signal possibly from the early universe, that was treated a signal possibly from the early universe,  such as inflationary gravitational waves \cite{Kuroyanagi:2020sfw,Tahara:2020fmn,Cai:2020qpu,Zhou:2020kkf}, cosmic strings and domain walls \cite{Ellis:2020ena,Buchmuller:2020lbh,Samanta:2020cdk,Liu:2020mru,Paul:2020wbz}, 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} and/or magnetic fields \cite{Neronov:2020qrl} and turbulence \cite{Abe:2020sqb}. 

In this paper we present a self-consistent study of the gravitational
waves from the turbulent sources possibly presented at QCD phase transitions.
In contrast to Ref.~\cite{Neronov:2020qrl} we will not limit
the turbulence sources assuming the magnetic fields presence, but
rather considering {\it any} turbulent source possibly presented at
QCD energy scale.
Even if the primordial fields are not dynamically strong,
turbulence can still develop during at QCD energy scale
Refs.~\cite{Rischke:2003mt,Kiskis:2003rd,Aoki:2006we}; and the latent
%AB: Tina, could yu insert the Kiskis:2003rd reference?
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{Miniati:2017kah}.
This approach will allow us to avoid the requirement of the first order
QCD phase transitions.

\section{Modeling}

\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$, 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 $\emf$ 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+\emf),\quad
\JJ=\nab\times\BB.
\nonumber
\end{eqnarray}
Here 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$ 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, $\emf$, is used to model the magnetic field
generation with
\begin{equation}
\emf(\xx,t)={\rm Re}\{{\cal N}\tilde{\ff}(\kk(t))\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}$ is a normalization factor 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 an arbitrarily
short time interval $1\leq t\leq 2$.

\subsection{Gravitational waves from turbulent sources}

We solve the non-dimensional GW equation in the radiation era for
the scaled strains $\tilde{h}_{+/\times}$, where the subscripts $+$
and $\times$ denote the plus and cross polarization modes, $\kk$ is
the comoving wave vector, normalized by the corresponding scale
$1/(c t_\ast)$ \cite{Pol:2018pao,Pol:2019yex}:
\begin{equation}
\left( \partial^2_t + \kk^2\right)
\tilde{h}_{+/\times} (\kk, t)\, =\, {6 \over t}
\tilde{T}_{+/\times}^{\rm TT} (\kk, t),
\label{GW}
\end{equation}
where $\tilde{T}_{+/\times}^{\rm TT}$ is the comoving stress tensor,
projected into the traceless and transverse (TT) gauge, described by
the linear polarization modes, and normalized by the energy density
at $t_\ast$.

The vigor of turbulence is characterized by the Reynolds number,
\EQ
\Rey=\urms/\nu\kf,
\EN
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.



\section{NanoGrav data}
Pulsar timing array (PTA) experiments typically characterize a stochastic GW background in terms of its characteristic strain spectrum $h_c(f)$ fitted with a power-law dependence on frequency:
\begin{equation}\label{eq:hfit}
    h_c(f)=A_\mathrm{CP}\bigg(\frac{f}{f_\mathrm{yr}}\bigg)^{\alpha_\mathrm{CP}}
\end{equation}
where $f_\mathrm{yr} = 1 \mathrm{yr}^{-1}$ is the reference frequency and the spectral index $\alpha_\mathrm{CP}$ depends on the source of the stochastic GW background.

The energy density spectrum of the GW background expressed in terms of the characteristic strain is
\begin{equation}
    \Omega_\mathrm{GW}(f)=\frac{2\pi^2}{3H_0^2}f^2 h_c^2(f).
\end{equation}
 
The GW density parameter may then be expressed with the power-law dependence on frequency:
\begin{equation}\label{eq:OmegaGWf}
    \Omega_\mathrm{GW}(f)=\frac{2\pi^2}{3H_0^2}A^2_\mathrm{CP}f_\mathrm{yr}^2\bigg(\frac{f}{f_\mathrm{yr}}\bigg)^{5-\gamma_\mathrm{CP}}
\end{equation}
where $\alpha_\mathrm{CP}=(3-\gamma_\mathrm{CP})/2$.

PTA results are reported via joint $A_\mathrm{CP}-\gamma_\mathrm{CP}$ posterior distributions. Posteriors for a common-spectrum process in the NANOGrav 12.5-yr 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 Ref. \cite{Arzoumanian:2020vkk} Figure 13), and so 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 \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 $\Omega_\mathrm{GW}(f)$ as given by equation (\ref{eq:OmegaGWf}) with $H_0 = 70$ km/s/Mpc. The shaded regions show the 2$\sigma$ confidence contours of the $A_\mathrm{CP}-\gamma_\mathrm{CP}$ parameter space in terms of $f$ and $\Omega_\mathrm{GW}(f)$ for frequencies from 2.5--100 nHz.%roughly corresponding to the NANOGrav sensitivity range 
\begin{figure} 
\centering
\includegraphics[scale=0.8]{fig/NGContours.eps}
\caption{NANOGrav 12.5-year data set $2\sigma$ confidence contours for the posteriors of a common-spectrum (CP) process (see Ref. \cite{Arzoumanian:2020vkk} Figure 1) expressed in terms of the gravitational wave energy density parameter $\Omega_\mathrm{GW}(f)$ and frequency $f$. This is shown over the frequency range of 2.5--100 nHz. The three models used to fit the CP process include a: broken power law (blue), 5-frequency power law (orange), and 30-frequency power law (green).}
\label{fig:NGcontours}
\end{figure}







\section{Numerical simulations}

We solve the governing equations using the {\sc Pencil Code} \cite{PC},
where the GW solver has already been implemented \cite{Pol:2018pao}.
In the following, we discuss the results for different values of $\kf$.
At sufficiently high Reynolds numbers, the temporal growth of $\EEM(t)$
is similar for different values of $\kf$; see \Fig{pcomp512}(b).
However, for smaller values of $\kf$, the turbulence decay is slower.
Furthermore, the GW energy increases with increasing $\kf$.
This is seen more clearly in a diagram of $\EEGW$ versus $\EEM/\kf$.

It is important to realize that in all four cases, the decay of the
magnetic energy follows an approximate powerlaw decay with
\EQ
\left.\EEM\approx\EEM^{\max}\right/\left(1+\delta t/\tau\right)^p,
\EN
where $\Delta t=t-2$ is the time after which the
electromagnetic source is turned off, $\tau$ is the decay time scale,
and $p$ is the exponent describing the decay.
The values of $\tau$ and $p$ are also listed in \Tab{Tsummary512}.

\begin{table}[htb]\caption{
All runs with $\eta=5\times10^{-5}$ and $512^3$ meshpoints.
}\vspace{12pt}\centerline{\begin{tabular}{cccccccc}
$\kf$ & $k_1$ & $f_0$ & $n$ & $\tau$ & $\EEM^{\max}$ & $\EEGW^{\rm sat}$ & $\hrms^{\rm sat}$ \\%run
\hline
 2 & 0.3& 1.9e-01 & 1.0 & 0.43 & 3.832e-02 & 3.529e-04 & 4.827e-02 \\%M1024sig0_k6_wav03_1p9em1b
 6 &  1 & 6.0e-02 & 1.0 & 2.0  & 3.748e-02 & 5.608e-05 & 7.063e-03 \\%M512sig0_k6_wav1b
20 &  3 & 2.3e-02 & 1.3 & 4.5  & 3.808e-02 & 1.102e-05 & 1.162e-03 \\%M512sig0_k6_wav3_2p3em2
60 & 10 & 1.0e-02 & 1.4 & 16   & 3.934e-02 & 2.620e-06 & 1.648e-04 \\%M512sig0_k6_wav10_1p0em2
\label{Tsummary512}\end{tabular}}\end{table}

To compare with earlier work, we show in \Fig{EEGW_vs_EEKM} the positions
of our runs in a $\Omega_{\rm GW}^{\rm sat}$ versus $\Omega_{\rm M}^{\max}$
diagram.
For orientation, we also show the data points from the Ref.~\cite{Pol:2019yex}.
We see that the new data points are well above the older once.
This is mainly a consequence of using now smaller values of $\kf$
(2--60, compared to 600 in Ref.~\cite{Pol:2019yex}).

%\begin{itemize}\itemsep-0.2mm
%\item
%\end{itemize}
Hypothesis:
\begin{equation}
\EEGW=\left(q\EEM/\kf\right)^2
\end{equation}
Determine $q$; we find $q\propto\kf^{1/4}$; see \Fig{pres}.
We see that the data can be well described by a fit of the form
$\EEGW=1.2\times10^{-1}\,(\EEM/\kf)^{3/2}$.
Only the data point for $\kf=60$ is slightly above the line.
This could be an artefact of using a finite Reynolds number
in our simulations.

The GW energy in the present simulations is high by comparison with
earlier work of Ref.~\cite{Pol:2019yex}; see \Fig{EEGW_vs_EEKM}.

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp512}
\end{center}\caption[]{
Black, blue, and red are for $\kf=2$, 6, and 20, respectively.
}\label{pcomp512}\end{figure}

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

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pres}
\end{center}\caption[]{
$\EEGW=5.2\times10^{-4}\,(\EEM/\kf)^{1/2}$ (black dotted line), and
$\EEGW=1.2\times10^{-1}\,(\EEM/\kf)^{3/2}$ (red line).
}\label{pres}\end{figure}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{EEGW_vs_EEKM}
\end{center}\caption[]{
Positions of the new runs in a diagram showing
$\Omega_{\rm GW}^{\rm sat}$ versus $\Omega_{\rm M}^{\max}$,
For orientation the old data points of the Ref.~\cite{Pol:2019yex}
are shown as gray symbols.
The red open symbols are for the runs with $512^3$ meshpoints.
}\label{EEGW_vs_EEKM}\end{figure}

In \Tab{Tsummary512}, we have quoted the values of $\EEGW^{\rm sat}$ and
$\hrms^{\rm sat}$ obtained at the end of the simulation $t=t_{\rm end}$.
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 (t_{\rm end}/a_0)^4$;
see Refs.~\cite{Pol:2019yex,Pol:2018pao} for details.
For the aforementioned fiducial parameters of $g_*=100$ and
$T_*=100\GeV$, this factor is $1.64\times10^{-5}$.
The largest value of $\EEGW^{\rm sat}$ quoted in
\Tab{Tsummary_ramp} is $2.4\times10^{-7}$ and corresponds
therefore to $h_0^2\OmGW=4\times10^{-12}$.
Likewise, the values of $\hrms^{\rm sat}$ in \Tab{Tsummary_ramp}
have to be multiplied by $a_0^{-1}=8.0\times10^{-16}$ to
obtain the observable $h_c$ at the present time.
Again, the largest value of $\hrms^{\rm sat}=5.1\times10^{-6}$
corresponds therefore to the observable $h_c=4\times10^{-21}$.

\section{Conclusions}

\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]].

%inflation
\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}, no.10, 103527 (2020)
doi:10.1103/PhysRevD.102.103527
[arXiv:2010.03537 [astro-ph.CO]].

\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{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{Cai:2020qpu}
Y.~Cai and Y.~S.~Piao,
``Intermittent NEC violations during inflation and primordial gravitational waves,''
[arXiv:2012.11304 [gr-qc]].

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

\bibitem{Buchmuller:2020lbh}
W.~B\"uchmuller, 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]].

%phase transitions
\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]].

%QCD 
\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{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{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
  %[nucl-th/0305030].

\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
  %[hep-lat/0611014].

\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{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).
%  [astro-ph/9602031].

\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\ygafd{2020}{114}{130}
{The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence}

\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
  
\bibitem{PC}
The {\sc Pencil Code}, \url{https://github.com/pencil-code},
DOI:10.5281/zenodo.2315093

\end{thebibliography}

\appendix
\section{Earlier low-res data}

\begin{figure}[h!]\begin{center}
\includegraphics[width=\columnwidth]{pcomp}
\end{center}\caption[]{
Black, blue, and red are for $\kf=2$, 6, 20, and 60, respectively.
All runs with $\eta=2\times10^{-3}$ and $64^3$ meshpoints.
}\label{pcomp}\end{figure}

\begin{table}[htb]\caption{
All runs with $\eta=2\times10^{-3}$ and $64^3$ meshpoints.
}\vspace{12pt}\centerline{\begin{tabular}{lccccccc}
$k_1$ & $f_0$ & $\EEM^{\max}$ & $\EEGW^{\rm sat}$ \\%run
\hline
 0.3& 2.0e-01 & 3.985e-02 & 2.762e-04\\%M64sig0_k6_wav03_2p0em1
  1 & 6.0e-02 & 3.874e-02 & 4.189e-05\\%M64sig0_k6_wav1d
  3 & 3.2e-02 & 4.038e-02 & 2.313e-05\\%M64sig0_k6_wav3_3p2em2
 10 & 2.3e-02 & 3.831e-02 & 1.375e-05\\%M64sig0_k6_wav10_2p34em2a
\label{Ttimescale}\end{tabular}}\end{table}
\section{Conclusions}


\vfill\bigskip\noindent\tiny\begin{verbatim}
$Header: /var/cvs/brandenb/tex/emma/k1GW/notes.tex,v 1.8 2021/01/24 05:18:52 brandenb Exp $
\end{verbatim}

\end{document}
