\documentclass[prl,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\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}[3]{, Geophys. Astrophys. Fluid Dynam. {\bf #2}, #3 (#1).}
\newcommand{\pgafd}[2]{, Geophys. Astrophys. Fluid Dynam., in press, arXiv:#2 (#1).}
\newcommand{\sgafd}[2]{, Geophys. Astrophys. Fluid Dynam., submitted, arXiv:#2 (#1).}
\newcommand{\yzfa}[3]{, Zeitschr. f. Astrophys. {\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{Circular Polarization of Gravitational Waves from Early-Universe Helical Turbulence}

\date{\today}
\preprint{NORDITA-2020-102}

\author{Tina~Kahniashvili}
\email{tinatin@andrew.cmu.edu}
%\affiliation{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, PA 15213, USA}
\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, 3-5 Cholokashvili St., 
%0194 Tbilisi, Georgia}
\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia}
%\affiliation{Abastumani Astrophysical Observatory, 
%M. Kostava st. 47/57, 
%Tbilisi, GE-0179, Georgia}
\affiliation{Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia}
%\affiliation{Department of Physics, Laurentian University, Ramsey Lake Road, Sudbury, ON P3E 2C, Canada}
\affiliation{Department of Physics, Laurentian University, Sudbury, ON P3E 2C, Canada}

\author{Axel~Brandenburg}
%\email{brandenb@nordita.org}
%\affiliation{Nordita, KTH Royal Institute of Technology and Stockholm University, 
%Roslagstullsbacken 23, 
%10691 Stockholm, Sweden}
\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, 3-5 Cholokashvili St.,
%0194 Tbilisi, Georgia}
\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, 
%5000 Forbes Ave, 
%Pittsburgh, PA 15213, USA}
\affiliation{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA}

\author{Giga Gogoberidze}
%\email{Grigol.Gogoberidze@iliauni.edu.ge}
%\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili St., 
%0194 Tbilisi, Georgia}
\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia}


\author{Sayan~Mandal}
%\email{sayanm@andrew.cmu.edu}
%\affiliation{Physics and Astronomy Department, Stony Brook University, Stony Brook, New York 11794, USA}
%\affiliation{Physics and Astronomy Department, Stony Brook University, 100 Nicolls Road, Stony Brook, New York 11794, USA}
\affiliation{Physics and Astronomy Department, Stony Brook University, Stony Brook, New York 11794, USA}
%\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili St., 
%0194 Tbilisi, Georgia}
\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia}

\author{Alberto~Roper~Pol}
%\email{alberto.roperpol@colorado.edu}
\affiliation{Universit\'e de Paris, CNRS, Astrophysique et Cosmologie, Paris, F-75013, France}
%\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili St.,
%0194 Tbilisi, Georgia}
\affiliation{Faculty of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia}

\begin{abstract}
We perform direct numerical simulations to compute the net circular
polarization of gravitational waves from helical (chiral) turbulent 
sources in the early universe for a variety of initial conditions, 
including driven (stationary) and decaying turbulence.
We investigate the resulting gravitational wave signal assuming different 
turbulence geneses such as magnetically or kinetically driven cases.
Under realistic physical conditions in the early universe we compute
numerically for the first time the total (integrated over all
wavenumbers) polarization degree of the gravitational waves. 
We find that the spectral polarization degree
strongly depends on the initial conditions. 
The peak of the spectral polarization degree occurs 
at twice the typical wavenumber of the source, as expected,
and for fully helical decaying turbulence, it reaches
its maximum (100\%) {\it only} at the peak.
We determine the temporal evolution of the turbulent sources as well as
the resulting gravitational waves, showing that the dominant contribution 
to their spectral energy density happens shortly after the source activation. 
Only through an artificially prolonged decay of the turbulence can further
increase of the gravitational wave amplitude be achieved. 
We estimate the detection prospects for the net polarization arguing that
its detection contains {\it clean} information (including the generation 
mechanisms, time and strength) about the sources of possible parity
violations in the early universe.
\end{abstract}
\pacs{98.70.Vc, 98.80.-k}

\maketitle

%\section{Introduction}
A remarkable possible source of stochastic gravitational waves (GWs) 
from the early universe is turbulence in the primordial plasma 
induced either from cosmological first-order (electroweak or QCD) phase transitions
\cite{Witten:1984rs,Hogan:1986qda,Kamionkowski:1993fg},
or from the primordial magnetic fields that are coupled to the cosmological plasma \cite{Brandenburg:1996fc,Christensson:2000sp,Kahniashvili:2010gp,Brandenburg:2017rnt}.
The GW signal is potentially detectable 
by the Laser Interferometer Space Antenna (LISA) in the case of
strong enough turbulent motions present at the electroweak scale
(assuming that the total energy in turbulence is up to 1--10\% 
of the total thermal energy at the moment of generation)
\cite{Kosowsky:2001xp,Caprini:2006jb,Gogoberidze:2007an,Kahniashvili:2008pe,Caprini:2009yp}. 
Interestingly, if the GWs are produced during the QCD phase transitions,
they are potentially detectable through Pulsar Timing Arrays (PTAs)
and this has been proposed to explain the reported NANOGrav signal
\cite{Neronov:2020qrl}.
Since GWs propagate {\it almost} freely from the moment of generation 
until today\footnote{We discard the GW damping due to neutrino free 
streaming \cite{Durrer:1997ta,Weinberg:2003ur} or from anisotropic
stresses \cite{Deryagin:1986qq}.}, the detection of GWs sourced by 
primordial turbulence will open a new window to understand physical 
processes in the very early stages of the evolution of the universe (at the 
femtoseconds timescale); see Ref.~\cite{Caprini:2019egz}
and references therein.
Moreover, several theoretical extensions of the standard model (SM)
of particle physics and cosmology (which is insufficient to explain the
matter-anti-matter asymmetry in the universe) imply
parity symmetry violation at the electroweak energy scale being
possibly manifested through helical (chiral) turbulent motions and/or
magnetic fields \cite{Long:2013tha,Dorsch:2016nrg}.
As expected, such parity-violating turbulent
sources will produce circularly polarized GWs
\cite{Kahniashvili:2005qi,Kahniashvili:2008er,Kisslinger:2015hua,
Alexander:2018fjp,Anand:2018mgf,Niksa:2018ofa,Ellis:2020uid}
analogously to the GWs produced via Chern-Simons coupling
\cite{Lyth:2005jf},
and, if detected, the GW polarization can be a {\it
clean} measure of the deviations from the SM and will provide us with a
better understanding of the nature of parity symmetry and its violation.

The detection of circular polarization of the stochastic GW 
background is a challenging task \cite{Seto:2006hf,Smith:2016jqs}, and
the planar interferometers cannot measure net polarization in the case
of isotropic backgrounds \cite{Seto:2007tn,Seto:2008sr}.
However, the dipolar anisotropy induced by our proper motion
with respect to the cosmic reference frame makes it possible to
measure the net circular polarization of the stochastic GW background
\cite{Seto:2006dz,Crowder:2012ik}, and recently it has been shown that 
the net polarization of GWs could be detected with a signal-to-noise ratio 
of order one by LISA if the strength of the signal achieves $h_0^2
\Omega_{\rm GW}\sim 10^{-11}$ (with $\Omega_{\rm GW}$ the fraction
between the GW energy density ${\mathcal E}_{\rm GW}$ and the critical
density today ${\mathcal E}_{\rm cr}=3H_0^2/(8\pi G)$, with
$H_0=100 h_0{\rm km} {\rm s}^{-1} {\rm Mpc}^{-1}$ the Hubble parameter today and
$G$ is the gravitational constant) \cite{Domcke:2019zls}. 
These findings make it extremely important to properly compute all 
characteristics (such as the amplitude, the spectral shape, and the
polarization degree) of the GW signal from primordial helical (chiral)
sources.

When computing the GW signal from early-universe turbulent sources,
previous studies (with the exception of Ref.~\cite{Pol:2019yex}) assumed
stationary hydrodynamic turbulence and a turbulence duration set by a
fraction of the Hubble time at the moment of generation ($H_\star^{-1}$),
making it possible to use the simplified GW equation with a discarded
term that describes the expansion of the universe 
($\sim H=a^{-1}{da/dt_{\rm ph}}$, where $t_{\rm ph}$ denotes the physical
time, and $a$ is the scale factor). 
Moreover, the magnetic field and primordial plasma coupling, and 
correspondingly, the turbulence decay have been neglected, making it
impossible to study the GW source dynamics and the temporal dependence
of the GW amplitude and spectral characteristics beyond the dilution due
to the universe expansion, i.e., ${\mathcal E}_{\rm GW} \propto a^{-4}$.

%\section{Model}
Recently we have performed for the first time direct numerical simulations
of magnetohydrodynamic (MHD) turbulence in the early universe accounting
for the expansion of the universe using the {\sc Pencil Code} \cite{PC},
and numerically computed the resulting stochastic GW background and
relic magnetic fields \cite{Pol:2019yex}. 
To scale out expansion effects, we use appropriately scaled comoving 
variables and conformal time.
The full set of MHD equations is then similar to the usual MHD equations
\citep{Brandenburg:1996fc}, and the GW equation is written for the scaled
strains and the comoving total (magnetic and kinetic) traceless-transverse 
stress tensor $T_{ij}^{\rm TT}$.
These simulations allowed us to conclude that the proper inclusion of
coupling effects results in the extension of the GW spectrum at lower frequencies 
due to the power transfer at large scales,
until the causal horizon determined by the comoving Hubble frequency,
$f_H=1.65 \cdot 10^{-5}{\rm Hz} \, (g_\star/100)^{1/6} (T_\star/100{\rm GeV})$
with $g_\star$ and $T_\star$ being the relativistic degrees of freedom
and the temperature at the moment of the GW source activation.
In terms of the GW energy density per logarithmic frequency interval
as a fraction of critical density today, $h_0^2\Omega_{\rm GW}(f)$, the
spectrum in the frequency range $f_{\rm H} \leq f \leq f_{\rm S}$ (where 
the typical frequency of the source is 
$f_{\rm S}=2 N f_{\rm H}$\footnote{The factor ``2'' is due to the quadratic nature of
the turbulent source.}, and $N$ determines the number of turbulent 
eddies per linear Hubble length scale)\footnote{The parameter 
$N$ is determined by the physical stirring wave number of the source, 
$k_0^{\rm phys}$, via $H_\star/N =2\pi/k_0^{\rm phys}$, and to the
comoving peak wave number through $k_0=k_0^{\rm phys} (a_\star/a_0)$.
Here $a_\star$ and $a_0$ are scale factors corresponding to the 
moment of generation and today, and $a_\star/a_0 = 8.0 \times 10^{-16}
(100/g_\star)^{1/3}(100{\rm GeV}/T_\star)$.} is $\propto f$ as
opposed to $f^{3}$ for the causal low-frequency ($f<f_H$) tail obtained analytically
\cite{Kahniashvili:2008pe}, 
while in terms of the {\it comoving} dimensionless strain
amplitude $h_c$, conventionally written as
$h_c(f) = 1.263 \cdot 10^{-18} (f/1{\rm \, Hz})^{-1} \left [h_0 \Omega_{\rm 
GW}(f)\right]^{1/2}$ \cite{Romano:2016dpx},
we observe the scaling $h_c \propto f^{-1/2}$
in the frequencies range $f_H < f < f_S$.
At this point we distinguish three parts of
$\Omega_{\rm GW}(f)$: the low-frequency region
below causal horizon $f<f_H$, 
the intermediate region $f_H < f < f_S$, and 
the high-frequency region $f>f_S$. 
However, due to computational limitations, we are unable
to reproduce the entire spectrum in our numerical simulations. 
At high frequencies, $f>f_{\rm S}$, numerical simulations agree
well with the analytical estimate; see Ref.~\cite{Pol:2019yex} for 
more details: in particular, for Kolmogorov turbulence with a slope $-5/3$,
the high frequency tail $h_c(f)$ ($\Omega_{\rm GW}(f)$) scales as 
$\propto f^{-7/3}$ ($\propto f^{-8/3}$). 
Earlier work \cite{Niksa:2018ofa} showed that the high-frequency
scaling of the GW spectrum depends strongly on the assumptions about
the turbulence and the modeling of the time-decorrelation function.
Since the spectrum of the stochastic GW background 
determines the detectability of GWs in a given experiment, realistic
turbulence simulations are essential for establishing the sensitivity of 
upcoming GW experiments to early-universe physics \cite{Romano:2016dpx}.

\begin{figure*}[t!]\begin{center}
\includegraphics[width=.9\textwidth]{pe_comp}
\end{center}\caption[]{\small
Evolution of magnetic energy (left) and growth of GW energy density
(right) for simulations where the driving is turned off at $t=1.1$ 
(black dotted line), or the strength of the driving is reduced linearly
in time over the duration $\tau=0.2$ (green),  $0.5$ (blue),
$1$ (red), or $2$ (black). 
Time is in units of the Hubble time at the moment of source activation.
}\label{pe_comp}\end{figure*}

\begin{figure*}[t!]\begin{center} 
\includegraphics[width=\textwidth]{pcomp_M512k2_both_fixed}
\end{center}\caption[]{\small 
Evolution of (a) $\OmK$, (b) $\epsK$, and (c) $\OmGW$ for kinetically
driven cases with $\sigma=0$ (black), $0.5$ (blue), and $1$ (red),
and of (d) $\OmM$, (e) $\epsM$, and (f) $\OmGW$ for magnetically
driven cases with $\sigma=0$ (black), $0.3$ (blue), and $1$ (red).
}\label{pcomp_M512k2}\end{figure*}

Interestingly, the spectral form of the GW spectrum at low frequencies
is independent of the initial conditions and the turbulence model while
the amplitude of the signal strongly depends on the model chosen.
The universal form of the GW spectrum does not allow us to discriminate
between helical and nonhelical sources and thus limits our ability to
determine the parity violation in the early universe. 
This leads us to the present study with its main focus on GW circular
polarization estimates and the question whether the detection of
polarization can help in the identification of distinct properties of
the source.

As it was shown in Ref.~\cite{Pol:2019yex}, the GW spectrum becomes
stationary shortly after the driving of the source ends (i.e.,
when the free decay stage of the source starts),
while the energy density of the source is still present. 
To demonstrate this, we drive magnetic fields with an electromotive
force consisting of plane waves that are delta-correlated in time.
In Fig.~\ref{pe_comp}, we show the temporal 
evolution of the source and the growth of the GW energy density 
for the driven ($1<t<1.1$) and decaying stages ($t>1.1$), where 
the driving decreases linearly for a time $\tau=0.1$--$2$,
although $\tau>0.5$ may be unrealistic.
During the statistically stationary stage, the GW energy density growth rate is 
proportional to the duration of turbulence, as was estimated 
through analytical modeling of 
Ref.~\cite{Gogoberidze:2007an}.
In reality, the driving stage is short compared to the 
Hubble time-scale, and consists of the few largest eddy turnover times.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pgrah_comp}
\end{center}
\caption[]{ 
Degree of circular polarization for (a) kinetically
and (b) magnetically forced cases with $\sigma=0$ (black)
$0.1$ (blue), $0.3$ (green), $0.5$ (orange), and $1$ (red). 
Approximate error bars based on the temporal fluctuations and
statistical spread for different random seeds of the forcing
are shown as solid black lines for
$\sigma=0$ and as dotted lines otherwise.
}\label{pgrah_comp}\end{figure}


Moreover, the GW generation can be split into three cases
(see Fig.~\ref{pcomp_M512k2} for kinetically and magnetically driven turbulence):
non-helical, partially helical (the fractional helicity is described 
by the ratio $\sigma<1$ of actual to
maximally allowed helicity), and fully helical. 
Surprisingly, the kinetically driven turbulence is more efficient in producing
GW energy; see Figs.~\ref{pcomp_M512k2}(c) and (f).
However in this case the presence of kinetic helicity does 
not affect the source amplitude -- contrary to the magnetically 
driven case where the amplitude 
of the source increases substantially with increasing fractional helicity;
see Figs.~\ref{pcomp_M512k2}(a) and (d). 
We also present the kinetic and magnetic energy dissipation rates,
$\epsK$ and $\epsM$, respectively.
The dissipation rates remain almost unchanged during the
stationary stage as we can expect. In addition, 
we see that they are almost unaffected by the presence of helicity.
One may have expected a correlation between $\epsK$ (or $\epsM$)
and $\Omega_{\rm GW}(f)$, but in the magnetically driven case, 
larger values of $\sigma$ produce even slightly less dissipation at
early times.
Nevertheless, $\Omega_{\rm GW}(f)$ clearly increases with $\sigma$.

To estimate the polarization degree, we follow the standard procedure 
described in Refs.~\cite{Romano:2016dpx,SM1}. 
Using the usual circular polarization basis tensors,
we decompose the Fourier transform of the GW strains 
into two states -- right- ($h_+$) and left-handed ($h_-$) 
circularly polarized GWs, and the 
polarization degree is defined  ${\mathcal P}(k)={\mathcal H}(k)/H(k)$, 
where $H({k})=\langle h^{\star}_+({\mathbf k}) h_{+}({\mathbf k'})
-h^{\star}_-({\mathbf k}) h_{-}({\mathbf k'}) \rangle$
and ${\mathcal H}({k}) = \langle h^{\star}_+({\mathbf k}) h_{+}({\mathbf k'}) +
h^{\star}_-({\mathbf k}) h_{-}({\mathbf k'}) \rangle$ characterize 
the GW amplitude and polarization, respectively.

The GW polarization degree depends, as expected,
on the fractional helicity of the source.
On the other hand, the evolution of helical sources is determined by the
initial fractional helicity: the coupling between the partially helical magnetic
field and the plasma motions leads to a reconfiguration of the magnetic
field at large scales through free decay, resulting in the growth of
the fractional helicity due to the increase of the correlation length
and the corresponding decrease of the magnetic energy until the fully
helical stage is developed and inverse cascading starts
\cite{Tevzadze:2012kk}.
However, for weakly helical sources, a substantial time-period is needed
to reach a fully helical configuration.
On the other hand, our simulations show that the dominant contribution
to the GW signal occurs shortly after the source has reached its maximum:
the subsequent decay of the magnetic field causes a
decline of the turbulent driving of GWs.
This decline is further enhanced by the expansion of the universe,
although this effect is small if the decay time of the turbulence is 
short compared with the Hubble time; see \Fig{pe_comp}. 
We see that, even with a substantially extended decay phase of the
turbulence of twice the Hubble time, the final GW production is 
enhanced by only a factor of about four.
Thus, the GW polarization will retain information about the
{\it initial} fractional helicity of the source.

As we have highlighted above, previous works
\cite{Kahniashvili:2005qi,Kisslinger:2015hua} considered stationary
turbulence with two different models for helical turbulence realization:
(i) Kolmogorov-like helical turbulence with two different spectral 
slopes, ${-5/3}$ and $-11/3$ for the spectral energy 
and helicity densities, $E_{\rm M}(k)$ and $H_{\rm M}(k)$, respectively,
and (ii) helicity-transferring
turbulence (if helicity transfer and small-scale helicity dissipation
dominates \cite{Kraichanan:1973}) with the spectral indices 
$-7/3$ and $-10/3$ for $E_{\rm M}(k)$ and $H_{\rm M}(k)$ 
\cite{Moiseev:1996}\footnote{Note that Ref.~\cite{Kahniashvili:2005qi}
used the different convention for the spectral indices referring to the
power spectra for the symmetric ($P_S(k) \propto E_{\rm M}(k)/k^2 $)
and the antisymmetric parts ($P_A(k) \propto H_{\rm M}(k)/k$).}.
The former case seems most suitable for describing the usual nonhelical
turbulence experiencing forward cascading \cite{Borue:1997}.
The difference in using these spectral shapes is determined by the 
effect of helicity on the energy dissipation length.
For highly helical turbulence, the helicity dissipation length 
is larger, so the two region description \cite{Kahniashvili:2005qi}
might be justified using Kolmogorov-like turbulence at large
wavenumbers and approximating the low wave number tail with helicity
transfer turbulence \cite{Ditlevsen:2001}. 
Following this description, Ref.~\cite{Kahniashvili:2008er} and later also 
Ref.~\cite{Ellis:2020uid} discussed two stages --
first the fractionally helical one and later the fully helical one with 
inverse transfer to describe more precisely the GW generation by
helical MHD turbulence.
However, these estimates suffer due to (i) the assumption of stationary
turbulence; (ii) neglecting the decay and correspondingly temporal
dynamical effects (the GW spectrum becomes stationary shortly after
source activation).

We show the polarization degree spectra in \Fig{pgrah_comp} for 
continuous pumping of kinetic or magnetic energy and 
helicity at intermediate scales.
We see a substantial difference between the kinetic and magnetic 
initial sources at the low frequency tail due to the inverse cascade 
for the magnetic sources that is absent in the kinetically driven case: 
more precisely the energy density spectrum is unchanged after the decay 
stage starts, while the transfer of magnetic helicity to the large scales
results in the increasing of the polarization degree. 
These results confirm that the polarization degree 
is scale dependent: $\propto k^{-0.5}$ 
at large wavenumbers, which is shallower than the $k^{-1}$
expected for Kolmogorov-like helical turbulence with different
spectral indices for the magnetic spectral energy density ($n_S=-5/3$)
and the spectral helicity density ($n_H=-11/3$). 
In our simulations, the actual indices are a bit smaller \cite{SM1},
which also explains the shallower slope in the polarization degree.
The departure from the theoretical predictions is 
due to the assumption of a scale-independent 
time-decorrelation function 
for magnetic energy and helicity densities. 
Interestingly, the spectral shape of the polarization degree is 
independent of the actual indices for energy and helicity, but depends
on the difference between them 
\cite{Kahniashvili:2005qi}. 
Obviously, in a realistic case, the proper consideration of
time-decorrelation and its dependence on wave numbers is required.
In fact, even in the simplified description, different forms of the
time-decorrelation function for different models of turbulence (including
both compressible and incompressible cases) and its scale-dependence
leads to different scaling of the GW spectrum at high frequencies
\cite{Niksa:2018ofa}.

%\section{Conclusions}
In this Letter we present the first numerical simulations of the
circular polarization degree of GWs generated through parity
violating turbulent sources in the early universe.
We present our results for GWs generated at electroweak energy scales,
but the formalism is not limited to the specific moment
of GW generation, 
and can be adjusted to primordial turbulence sources 
at any time after inflation and before recombination epochs.
We have confirmed that the GW signal reaches its maximal strength
faster then the turbulence decays. 
We have also shown that the slope of the low frequency tail 
of the GW spectrum
is independent of the nature of the turbulent source (i.e., the 
nature of initial driving, the presence of helicity, etc). 
This restricts the
discrimination between helical and non-helical sources 
as well between kinetic and magnetic drivings
if the signal will be detected.
On the other hand, the polarization spectrum not only 
retains information about 
the initial characteristic frequency (as well as 
the GW signal does) and the strength of parity violation of the source,
but also manifests the dependence on the driving mechanisms.
We have shown that the previously used assumption of stationary turbulence
does not predict the GW polarization spectrum for realistic turbulence
(the scaling of the spectrum at low and high frequencies, the peak
position, etc), and thus might result in inadequate estimates of the
detection prospects.
Fortunately, numerical simulations have now become an affordable tool
to address these and other questions of relic GW generation 
and give a more complete picture for the detection prospects by 
LISA (as for electroweak phase transitions) and by PTAs (as for 
QCD phase transitions), and/or any future planned missions, 
including atomic interferometry \cite{Dimopoulos:2007cj}. 


%\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.~\cite{DATA}.


\acknowledgements
%\vspace{3mm}
We thank Arthur Kosowsky and Andrii Neronov 
for useful discussions. 
Support through the Swedish Research Council, grant 2019-04234,
and Shota Rustaveli GNSF  (grants FR/18-1462 and FR/19-8306)
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.

\vspace{-5mm}
%references
\begin{thebibliography}{}
\vspace{-5mm}
\bibitem{Kamionkowski:1993fg}
  M.~Kamionkowski, A.~Kosowsky and M.~S.~Turner,
 % ``Gravitational radiation from first order phase transitions,''
  Phys.\ Rev.\ D {\bf 49}, 2837 (1994).
  %[astro-ph/9310044].

\bibitem{Witten:1984rs}
  E.~Witten,
%  ``Cosmic Separation of Phases,''
  Phys.\ Rev.\ D {\bf 30}, 272 (1984).

\bibitem{Hogan:1986qda}
  C.~J.~Hogan,
%  ``Gravitational radiation from cosmological phase transitions,''
  %Mon.\ Not.\ Roy.\ Astron.\ Soc.\  
  MNRAS 
  {\bf 218}, 629 (1986).

\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{Christensson:2000sp}
  M.~Christensson, M.~Hindmarsh and A.~Brandenburg,
  %``Inverse cascade in decaying 3-D magnetohydrodynamic turbulence,''
  Phys.\ Rev.\ E {\bf 64}, 056405 (2001).
%  [astro-ph/0011321].

\bibitem{Kahniashvili:2010gp}
  T.~Kahniashvili, A.~Brandenburg, A.~G.~Tevzadze and B.~Ratra,
  %``Numerical simulations of the decay of primordial magnetic turbulence,''
  Phys.\ Rev.\ D {\bf 81}, 123002 (2010).
%  [arXiv:1004.3084 [astro-ph.CO]].

\bibitem{Brandenburg:2017rnt}
  A.~Brandenburg, {\it et al.}
  %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{Kosowsky:2001xp}
  A.~Kosowsky, A.~Mack and T.~Kahniashvili,
  %``Gravitational radiation from cosmological turbulence,''
  Phys.\ Rev.\ D {\bf 66}, 024030 (2002).
%  [astro-ph/0111483].


\bibitem{Caprini:2006jb}
  C.~Caprini and R.~Durrer,
  %``Gravitational waves from stochastic relativistic sources: Primordial turbulence and magnetic fields,''
  Phys.\ Rev.\ D {\bf 74}, 063521 (2006).
%  [astro-ph/0603476].


\bibitem{Gogoberidze:2007an}
  G.~Gogoberidze, T.~Kahniashvili and A.~Kosowsky,
  %``The Spectrum of Gravitational Radiation from Primordial Turbulence,''
  Phys.\ Rev.\ D {\bf 76}, 083002 (2007).
%  [arXiv:0705.1733 [astro-ph]].

\bibitem{Kahniashvili:2008pe}
T.~Kahniashvili, {\it et al. }
%T.~Kahniashvili, L.~Campanelli, G.~Gogoberidze, Y.~Maravin and B.~Ratra,
%``Gravitational Radiation from Primordial Helical Inverse Cascade MHD Turbulence,''
Phys. Rev. D \textbf{78}, 123006 (2008).
%doi:10.1103/PhysRevD.78.123006
%[arXiv:0809.1899 [astro-ph]].

\bibitem{Caprini:2009yp}
  C.~Caprini, R.~Durrer and G.~Servant,
  %``The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition,''
  JCAP {\bf 0912}, 024 (2009).
%  doi:10.1088/1475-7516/2009/12/024
%  [arXiv:0909.0622 [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]].

%AR: Added NANOGrav reference
%\bibitem{NANOGrav:2020}
%Z.~Arzoumanian, {\it et al. }
%[arXiv:2009.14174 [astro-ph]]

\bibitem{Durrer:1997ta}
R.~Durrer and T.~Kahniashvili,
%``CMB anisotropies caused by gravitational waves: A Parameter study,''
Helv. Phys. Acta \textbf{71}, 445 (1998).
%[arXiv:astro-ph/9702226 [astro-ph]].

\bibitem{Weinberg:2003ur}
S.~Weinberg,
%``Damping of tensor modes in cosmology,''
Phys. Rev. D \textbf{69}, 023503 (2004).
%doi:10.1103/PhysRevD.69.023503
%[arXiv:astro-ph/0306304 [astro-ph]].


\bibitem{Deryagin:1986qq}
D.~V.~Deryagin, D.~Y.~Grigoriev, V.~A.~Rubakov and M.~V.~Sazhin,
%``Possible Anisotropic Phases in the Early Universe and Gravitational Wave Background,''
Mod. Phys. Lett. A \textbf{1}, 593 (1986).
%doi:10.1142/S0217732386000750

\bibitem{Caprini:2019egz}
C.~Caprini, {\it et al.}
%M.~Chala, G.~C.~Dorsch, M.~Hindmarsh, S.~J.~Huber, T.~Konstandin, J.~Kozaczuk, G.~Nardini, J.~M.~No, K.~Rummukainen, P.~Schwaller, G.~Servant, A.~Tranberg and D.~J.~Weir,
%``Detecting gravitational waves from cosmological phase transitions with LISA: an update,''
JCAP \textbf{03}, 024 (2020).
%doi:10.1088/1475-7516/2020/03/024
%[arXiv:1910.13125 [astro-ph.CO]].

\bibitem{Long:2013tha}
    A.~J.~Long, E.~Sabancilar and T.~Vachaspati,
%    ``Leptogenesis and Primordial Magnetic Fields,''
    JCAP {\bf 1402}, 036 (2014).
    %[arXiv:1309.2315 [astro-ph.CO]].

\bibitem{Dorsch:2016nrg}
  G.~C.~Dorsch, S.~J.~Huber, T.~Konstandin and J.~M.~No,
%  ``A Second Higgs Doublet in the Early Universe: Baryogenesis and Gravitational Waves,''
  JCAP {\bf 1705}, 052 (2017).
  %doi:10.1088/1475-7516/2017/05/052
  %[arXiv:1611.05874 [hep-ph]]

\bibitem{Kahniashvili:2005qi}
T.~Kahniashvili, G.~Gogoberidze and B.~Ratra,
%``Polarized cosmological gravitational waves from primordial helical turbulence,''
Phys. Rev. Lett. \textbf{95}, 151301 (2005).
%doi:10.1103/PhysRevLett.95.151301
%[arXiv:astro-ph/0505628 [astro-ph]].

\bibitem{Kahniashvili:2008er}
  T.~Kahniashvili, G.~Gogoberidze and B.~Ratra,
  %``Gravitational Radiation from Primordial Helical MHD Turbulence,''
  Phys.\ Rev.\ Lett.\  {\bf 100}, 231301 (2008).
%  doi:10.1103/PhysRevLett.100.231301
%  [arXiv:0802.3524 [astro-ph]].

\bibitem{Kisslinger:2015hua}
L.~Kisslinger and T.~Kahniashvili,
%``Polarized Gravitational Waves from Cosmological Phase Transitions,''
Phys. Rev. D \textbf{92}, 043006 (2015).
%doi:10.1103/PhysRevD.92.043006
%[arXiv:1505.03680 [astro-ph.CO]].

\bibitem{Alexander:2018fjp}
  S.~Alexander, E.~McDonough and D.~N.~Spergel,
%  ``Chiral Gravitational Waves and Baryon Superfluid Dark Matter,''
  JCAP {\bf 1805}, 003 (2018).
  %doi:10.1088/1475-7516/2018/05/003
  %[arXiv:1801.07255 [hep-th]].

\bibitem{Anand:2018mgf}
  S.~Anand, J.~R.~Bhatt and A.~K.~Pandey,
%  ``Chiral plasma instability and primordial Gravitational wave,''
Eur.\ Phys.\ J.\ C {\bf 79},  119 (2019).
%  doi:10.1140/epjc/s10052-019-6619-5
%  [arXiv:1801.00650 [astro-ph.CO]].

\bibitem{Niksa:2018ofa}
P.~Niksa, M.~Schlederer and G.~Sigl,
%``Gravitational Waves produced by Compressible MHD Turbulence from Cosmological Phase Transitions,''
Class. Quant. Grav. \textbf{35},  144001 (2018).
%doi:10.1088/1361-6382/aac89c
%[arXiv:1803.02271 [astro-ph.CO]].

\bibitem{Ellis:2020uid}
J.~Ellis, {\it et al.}
%J.~Ellis, M.~Fairbairn, M.~Lewicki, V.~Vaskonen and A.~Wickens,
%``Detecting circular polarisation in the stochastic gravitational-wave background from a first-order cosmological phase transition,''
%[arXiv:2005.05278 [astro-ph.CO]].
JCAP \textbf{10}, 032 (2020).

\bibitem{Lyth:2005jf}
D.~H.~Lyth, C.~Quimbay and Y.~Rodriguez,
%``Leptogenesis and tensor polarisation from a gravitational Chern-Simons term,''
JHEP \textbf{03}, 016 (2005).
%doi:10.1088/1126-6708/2005/03/016
%[arXiv:hep-th/0501153 [hep-th]].

\bibitem{Seto:2006hf}
N.~Seto,
%``Prospects for direct detection of circular polarization of gravitational-wave background,''
Phys. Rev. Lett. \textbf{97}, 151101 (2006).
%doi:10.1103/PhysRevLett.97.151101
%[arXiv:astro-ph/0609504 [astro-ph]].

\bibitem{Smith:2016jqs}
  T.~L.~Smith and R.~Caldwell,
  %``Sensitivity to a Frequency-Dependent Circular Polarization in an Isotropic Stochastic Gravitational Wave Background,''
  Phys.\ Rev.\ D {\bf 95}, 044036 (2017).
  %doi:10.1103/PhysRevD.95.044036
  %[arXiv:1609.05901 [gr-qc]].

\bibitem{Seto:2007tn}
N.~Seto and A.~Taruya,
%``Measuring a Parity Violation Signature in the Early Universe via Ground-based Laser Interferometers,''
Phys. Rev. Lett. \textbf{99}, 121101 (2007).
%doi:10.1103/PhysRevLett.99.121101
%[arXiv:0707.0535 [astro-ph]].

\bibitem{Seto:2008sr}
  N.~Seto and A.~Taruya,
  %``Polarization analysis of gravitational-wave backgrounds from the correlation signals of ground-based interferometers: Measuring a circular-polarization mode,''
  Phys.\ Rev.\ D {\bf 77}, 103001 (2008).
  %doi:10.1103/PhysRevD.77.103001
  %[arXiv:0801.4185 [astro-ph]].

\bibitem{Seto:2006dz}
N.~Seto,
%``Quest for circular polarization of gravitational wave background and orbits of laser interferometers in space,''
Phys. Rev. D \textbf{75}, 061302 (2007).
%doi:10.1103/PhysRevD.75.061302
%[arXiv:astro-ph/0609633 [astro-ph]].

\bibitem{Crowder:2012ik}
S.G. Crowder, {\it et al.} 
%S.~G.~Crowder, R.~Namba, V.~Mandic, S.~Mukohyama and M.~Peloso,
%``Measurement of Parity Violation in the Early Universe using Gravitational-wave Detectors,''
Phys. Lett. B \textbf{726}, 66 (2013).
%doi:10.1016/j.physletb.2013.08.077
%[arXiv:1212.4165 [astro-ph.CO]].

\bibitem{Domcke:2019zls}
V.~Domcke, {\it et al.}
%J.~Garcia-Bellido, M.~Peloso, M.~Pieroni, A.~Ricciardone, L.~Sorbo and G.~Tasinato,
%``Measuring the net circular polarization of the stochastic gravitational wave background with interferometers,''
JCAP \textbf{05}, 028 (2020).
%doi:10.1088/1475-7516/2020/05/028
%[arXiv:1910.08052 [astro-ph.CO]].


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

\bibitem{PC}
The {\sc Pencil Code}, \url{https://github.com/pencil-code},
DOI:10.5281/zenodo.2315093

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

\bibitem{SM1}
See the Supplemental Material for the spectra of the symmetric and antisymmetric
parts of the GW amplitude and polarization, $H({k})$ and ${\mathcal H}({k})$,
respectively, as well as the underlying energy and helicity spectra.

\bibitem{Tevzadze:2012kk}
A.~G.~Tevzadze, L.~Kisslinger, A.~Brandenburg and T.~Kahniashvili,
%``Magnetic Fields from QCD Phase Transitions,''
Astrophys. J. \textbf{759}, 54 (2012). 
%doi:10.1088/0004-637X/759/1/54
%[arXiv:1207.0751 [astro-ph.CO]].

\bibitem{Kraichanan:1973}
R. H. Kraichnan, J. Fluid Mech., 59, 745 (1973).

\bibitem{Moiseev:1996}
 S. S. Moiseev and O.  Chkhetiani, 
%Sov. JETP, {\bf 83}, 192, (1996)
 JETP {\bf 83}, 192, (1996).
 
\bibitem{Borue:1997}
V. Borue and S. A. Orszag, 
%"Spectra in helical three-dimensional homogeneous isotropic turbulence"
Phys. Rev. E {\bf 55}, 7005 (1997). 
 \bibitem{Ditlevsen:2001}
 P. D. Ditlevsen and P. Giuliani, Phys. Rev. E {\bf 63}, 036304
(2001).

\bibitem{Dimopoulos:2007cj}
S.~Dimopoulos, P.~W.~Graham, J.~M.~Hogan, M.~A.~Kasevich and S.~Rajendran,
%``Gravitational Wave Detection with Atom Interferometry,''
Phys. Lett. B \textbf{678}, 37 (2009). 
%doi:10.1016/j.physletb.2009.06.011
%[arXiv:0712.1250 [gr-qc]].

\bibitem{DATA}
T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal, and A. Roper~Pol,
Datasets for ``Circular Polarization of Gravitational Waves from Early-Universe Helical Turbulence''
(v2020.11.07), \url{https://doi.org/10.5281/zenodo.4256906}.

\end{thebibliography}

\newpage 
\onecolumngrid
\appendix

\clearpage
\section*{\LARGE Supplemental Material}

\subsection{1 Polarization Degree}
To compute the polarization degree we follow the standard procedure:
choosing the coordinate system
so that unit vector ${\hat {\mathbf e}}_3$ points in the GW propagation
direction (${\bf \hat k} ={\bf k}/k$ with $k=|{\bf k}|$ is the unit
vector) and using the usual circular polarization basis tensors
$e^{\pm}_{ij} = -({\bf e}_1 \pm i{\bf e}_2)_i \times ({\bf e}_1 \pm i {\bf e}_2)_j/\sqrt{2}$,
we decompose the Fourier transform of the GW strains 
$h_{ij}({\bf k}) = \int d^3x \, e^{i{\bf k}\cdot {\bf x}} h_{ij}({\bf x})$
into two states -- right- ($h_+$) and left-handed ($h_-$) 
circularly polarized GWs $h_{ij}=h_+e^+_{ij} + h_- e^-_{ij}$.
The GW circular polarization degree is given by
\cite{Kahniashvili:2005qi}\footnote{As an alternative we can use the decomposition using the {\it linear} polarization basis tensors $e^+_{ij}({\bf \hat k})=e_i^1 e_j^1 -e_i^2 e_j^2$ and $e^\times_{ij}({\bf \hat k})=e_i^1 e_j^2 + e_i^2 e_j^1$ as  $h_{ij}({\bf k}) = h_+({\bf k}) e_{ij}^+({\bf \hat k}) + h_\times({\bf k}) e_{ij}^\times({\bf \hat k})$ where $h_+$ and $h_\times$ are gauge independent components corresponding to two polarization modes. In this case the 
polarization degree will be equal to 
${\mathcal P}(k) =   {\langle h^{\star}_\times ({\mathbf k})
h_{+}({\mathbf k'}) -
 h^{\star}_+({\mathbf k}) h_\times({\mathbf k'}) \rangle}/
{\langle h^\star_+({\mathbf k}) h_+({\mathbf k'}) +
 h^{\star}_\times({\mathbf k}) h_\times({\mathbf k'}) \rangle}$
} 
\begin{equation}
{\mathcal P}(k) =  \frac {\langle h^{\star}_+({\mathbf k})
h_{+}({\mathbf k'}) -
 h^{\star}_-({\mathbf k}) h_{-}({\mathbf k'}) \rangle}
{\langle h^{\star}_+({\mathbf k}) h_{+}({\mathbf k'}) +
 h^{\star}_-({\mathbf k}) h_{-}({\mathbf k'}) \rangle}
=\frac{{\mathcal H}(k)}{H(k)},~ \label{degree}
\end{equation}
where $H ({k})$ and ${\mathcal H}({k})$ characterize 
the GW amplitude and polarization (chirality) defined through the 
Gaussian-distributed GWs wave number-space two-point 
function (for simplicity of notations we omit the time-dependence):
$ 
\langle h^{\star}_{ij}({\mathbf k}) h_{lm} ({\mathbf k'})\rangle
=(2\pi)^3 \delta^{(3)}({\bf k}-{\bf k'}) \left[ {\mathcal
M}_{ijlm} H(k) + i{\mathcal A}_{ijlm} {\mathcal H} (k) \right],
$
with $4 {\mathcal M}_{ijlm} ({\mathbf{\hat k}}) \equiv
P_{il}P_{jm}+P_{im}P_{jl}-P_{ij}P_{lm}$, and
 $8 {\mathcal A}_{ijlm}({\mathbf{\hat k}})
\equiv {\hat {\bf k}}_q (P_{jm} \epsilon_{ilq} + P_{il}
\epsilon_{jmq} + P_{im} \epsilon_{jlq} + P_{jl} \epsilon_{imq})$ are tensors,
$P_{ij}({\bf \hat k}) \equiv \delta_{ij} - {\hat k}_i{\hat k}_j$
is the projector operator, $\delta_{ij}$ and $\epsilon_{ijl}$ are the
Kronecker delta and the fully antisymmetric tensor, respectively.

\subsection{2 Energy and helicity spectra}

In \Fig{pgrah_spec_comp} we show the numerator and denominator
of the degree of polarization, $H(k)$ and ${\cal H}(k)$, respectively.
The underlying kinetic and magnetic energy and helicity spectra are
shown in \Fig{phel_spec_comp}.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=.85\textwidth]{pgrah_spec_comp}
\end{center}\caption[]{
$H(k)$ (left) and $|{\cal H}(k)|$ (right)
for (a,b) kinetically and (c,d) magnetically forced cases
with $\sigma=0$ (black) $0.1$ (blue), $0.3$ (green),
$0.5$ (orange), and $1$ (red).
}\label{pgrah_spec_comp}\end{figure*}

\begin{figure}[h!]\begin{center}
\includegraphics[width=.85\columnwidth]{phel_spec_comp}
\end{center}\caption[]{\small 
(a) $2E_K(k)$ and (b) $|H_K(k)/k|$ 
for kinetically forced cases, and
(c) $2E_M(k)$, and (d) $|kH_M(k)|$
for magnetically forced cases
with $\sigma=0$ (black) $0.1$ (blue), $0.3$ (green),
$0.5$ (orange), and $1$ (red).
}\label{phel_spec_comp}\end{figure}

\FFig{pgrah_spec_comp} shows that the GW polarization is larger at
larger length scales.
To see whether this could be related to inverse cascading in the
magnetic case, we now show in \Fig{phel_spec_comp} the corresponding
energy and helicity spectra.
They show clear inverse cascading of the magnetic
energy and helicity spectra in the magnetically driven case.
Inverse cascading is clearly absent in the kinetic energy and helicity
spectra in the kinetically driven case.

\end{document}

