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

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

\newcommand{\neff}{N_{\rm eff}}
\newcommand{\dneff}{\Delta N_{\rm eff}}
\newcommand{\neffv}{N_{\rm eff}^{(\nu)}}
\newcommand{\nv}{N_{\nu}}
\newcommand{\dnv}{\Delta N_{\nu}}

\newcommand{\etaHe}{\eta_{\rm He}}
\newcommand{\etab}{\eta_{10}}
\newcommand{\Yp}{{\rm Y}_{\rm P}}

\newcommand{\B}{\mathbf{B}}
\usepackage{nicefrac}
\usepackage[hidelinks]{hyperref}

\begin{document}
%TK Title tentative

\title{Relic Gravitational Waves from Chiral Plasma Instability \\in the Standard Cosmological Model}

\date{\today}
\preprint{NORDITA-2023-0XX}



%TK author list is alphabetic for now

\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{School 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{Tina~Kahniashvili}
\email{tinatin@andrew.cmu.edu}
\affiliation{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA}
\affiliation{School of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia}
\affiliation{Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia}

\author{Andrew J. Long}
%\email{andrewjlong@rice.edu}
\affiliation{Department of Physics and Astronomy, Rice University, 6100 Main St., Houston, TX 77005, USA}

\author{Guotong Sun}
%\email{guotongs@andrew.cmu.edu}
\affiliation{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA}


\begin{abstract}
To explain the observed matter-antimatter asymmetry it is natural to assume the parity symmetry breaking in the primordial plasma.  
Preferred handedness in the early universe can be manifested though the chiral magnetic effect (CME). 
CME in highly conducting primordial plasma leads to development of a chiral plasma instability (CPI) and magnetohydrodynamic turbulence that generates seed (hyper)magnetic fields and sources the stochastic polarized gravitational wave background. 
In this letter we revisit the CME in the early universe, varying the total normalized chemical potential that is responsible for chiral asymmetry. 
We estimate the resulting magnetic field strength and we compute the stochastic GW background spectral characteristic, including the amplitude and polarization of the signal. 
Finally we discuss the GW detection prospects. 
\end{abstract}

\maketitle

One of the unresolved puzzles in modern cosmology is the observed matter-antimatter asymmetry (matter-antimatter problem). One of the natural possibilities to explain the matter excess over the antimatter today is to assume the baryon asymmetry has been originated in the early universe. Many models of baryogenesis implies the baryon asymmetry, and the parity (mirror) symmetry breaking that is one of the conditions required to produce matter abundance. One of possible manifestation of the parity symmetry breaking is chiral asymmetry. Despite the anomalous non-conservation of a chiral current (a so called chiral anomaly), is prohibited in a classical description, charge-parity (CP) non-conservation (e.g. violation) might be an indication of such a phenomenon in the early universe. Actualy chiral anomaly has been proposed as an explanation of the matter asymmetry.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\paragraph{\textbf{Motivation -- we know that the primordial plasma needed to carry a baryon asymmetry in order to explain the observed excess of matter over antimatter -- we don't know if there may have been other particle asymmetries as well -- such asymmetries do arise in many models of baryogenesis -- it's natural to investigate the observable implications of such an asymmetry -- testing for `by-products' of baryogenesis}} 
%%%%%%%%%%%%%%%%%%%%%%%%%%
\cite{Mavromatos:2018map}. 
One of the consequences of chiral anomaly is the {\it chiral magnetic effect} (CME) - the generation of the electric current along the {\it external} magnetic field \cite{Vilenkin:1980fu}. To account for this quantum effect it is important to describe adequately the early-universe plasma with charged 
fermions (for examples, the quark-gluon  \cite{Kharzeev:2013ffa} or the electroweak plasma \cite{Pavlovic:2016mxq}). 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\paragraph{\textbf{Context 1 -- in a QED plasma the presence of a chiral asymmetry (excess of right-chiral over left-chiral particles) is known to modify the electromagnetic equations of motion, introducing a new term in the electromagnetic current, known as the chiral magnetic effect -- upon extending the equations of magnetohydrodynamics by the CME term, one finds that the presence of the chiral instability opens a tachyonic instability (in long-wavelength modes) and leads to the development of a helical magnetic field -- while this is an interesting scenario to explore the origin of primordail magnetism in its own right, our focus is on the gravitational wave radiation sourced by the magnetic field's inhomogeneous energy density}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Upon extending the equations of magnetohydrodynamics (MHD) by the CME term, one finds that the presence of the chiral instability opens a tachyonic instability (in long-wavelength modes) and leads to the development of a helical magnetic field -- while this is an interesting scenario to explore the origin of primordial magnetism in its own right, our focus is to determine the CME cosmological consequences. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In fact, the conservation of total chirality of the system \cite{Bekenstein:1987} (fermion chirality plus any other parity violating quantity) may lead to different cosmological implications: in particular, one possible explanation for the observed large scale correlated magnetic fields in voids (see Ref.~\cite{Neronov:2010gir} for a pioneering work) is to assume the primordial magnetogenesis \cite{Vachaspati:2020blt} leading to immersing MHD turbulence in the early universe (that is a consequence of coupling between the (hyper)magnetic field and fluid  motions) with non-zero helicity (i.e. chirality). One example of magnetically dynamo-driven (through chiral magnetic instability (CPI), \cite{Akamatsu:2013pjd}) early-universe helical MHD turbulence \cite{Brandenburg:2017rcb} is that due to an asymmetry in the number density of left- and right-handed fermions \cite{Boyarsky:2011uy,Boyarsky:2012ex,Boyarsky:2015faa}. 
In this magnetogenesis scenario the conservation of total chirality implies that while magnetic helicity is immersing while the initial fermions chirality is disappearing \cite{Joyce:1997uy}. 


In what follows we re-address the CME effect in the early universe assuming that the primordial plasma possessed a nonzero initial chiral asymmetry (associated with fermions)  and vanishing magnetic field (i.e. magnetic  helicity). Under assumptions imposed by standard model (SM) of particle physics and cosmology \cite{Kamada:2016cnb,Domcke:2022uue}, we study the growth of a helical magnetic field via the CPI, and we calculate the resulting spectrum of gravitational wave background extending previous studies  of Refs. ~\cite{Brandenburg:2021aln}. 
%TK: more references?
One appealing aspect of our approach is that SM physics operates, such that properties of the primordial plasma are calculable as it is done in Refs. 
%TK references 
Our particular interest is to determine 
whatever CPI limited by the conventional SM parameters  may produce the viable strength of the magnetic fields today and the prospects to detect 
the corresponding gravitational radiation signal. 

We assume the cosmological parameters conventionally used in SM, and we introduce the conformal time $t$ and $t_\star$ -- the moment of the chiral (hyper)magnetic field generation. Along with the time-dependent number of relativistic degrees of freedom $g_\star(T)$ (equal to $106.75$ for the temperature higher than $T>100$ GeV while it decreases quickly at the late time 
with $T<100$ TeV \cite{Husdal:2016haj}) -- the only free parameter is the strength of the (quasi-conserved) chiral asymmetry, parametrised through the hypercharge-weighted chiral chemical potential $\mu_{Y,5}$ and normalized to the temperature $T$.\footnote{
Note, $\mu_{Y,5}$ is the the sum of chemical potentials for right-handed particles minus the sum for left-handed particles, but with different weighting factors,  which is directly related to the initial number of left ($n_{\rm L}$) and right handed ($n_{\rm R}$) fermions.}

Obviously being immersed through CPI, the strength of the field ($B \equiv |{\bf B}_Y|$ and the generation energy-scale (e.g. temperature) will be determined by the total initial chirality $\mu_{Y,5}$  under assumption that no magnetic field was presented for $t \leq t_\star$.  Note, we neglect the slight dependence on $g_\star$.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\paragraph{\textbf{Goal -- in this work we assume that the primordial plasma possessed a nonzero initial chiral asymmetry (and vanishing magnetic helicity), we study the growth of a helical magnetic field via the chiral plasma instability, and we calculate the resulting spectrum of gravitational wave radiation -- see also refs[...] for related recent work with different assumed initial conditions -- one appealing aspect of our approach is that only Standard Model physics operates, such that properties of the plasma are calculable (and in fact, already calculated) -- the only free parameter is the size of the (quasi-conserved) chiral asymmetry, parametrized through the hypercharge-weighted chiral chemical potential $\mu_{Y,5}$}} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%s 
Throughout the paper, we use natural units with  $\hbar=c=k_B=1$. 
We also set the permeability of free space to unity and 
expressing the electromagnetic quantities in Lorentz-Heaviside units such that the magnetic energy density is $B^2/2$. 
Unless otherwise specified, we refer to the comoving (rescaled) quantities and conformal time in order to account for the universe's expansion, and formulate all relevant equations in the 
convenient form \cite{Brandenburg:1996fc}. Notably the magnetic field 
${\bf B} \rightarrow {\bf B}a^2$, $T \rightarrow T a$, $n \rightarrow n a^3$ where $a\equiv a(t)$ is the scale factor. 

To determine the universe's expansion rate we account for the entropy conservation,
$g^{(S)} (aT)^3 ={\rm const}$,  
with $g^{(S)} \approx g_\star$ approximately equal to the entropic degrees of freedom for temperatures above 0.1 MeV. 
Recalling that the Hubble parameter, 
$H(t)=a^{-2}(t)da/dt$ 
satisfies the first Friedmann equation,  
$H^2 = 8\pi G \rho_{\rm rad}/3$ during the radiation-dominated (RD) epoch, and with $\rho_{\rm rad}=(\pi^2/30) g_\star T^4$ representing
the energy density associated with radiation (including all relativistic components), the rescalling $a_\star/a_)$ factor is given through: 
\begin{equation}
    \frac{a_\star}{a_0} = 8 \times 10^{-19} \Big( \frac{106.75}{g_\star}\Big)^{1/3}\Big( \frac{100 {\rm TeV}}{T_\star}\Big)    \label{scale factor}
\end{equation}

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

When considering the CME (e.g. CPI), the spin-flip reaction should be properly accounted for \cite{Kamada:2016eeb,Dvornikov:2020qwx}. 
Each of the
constituent comoving chemical potentials $\mu_i$ evolves in response to the various particle physics reactions
that involve particles of species $i$. In particular, inverse-decay reactions involving the Higgs tend to erase these
chemical potentials, and correspondingly, drive $\mu_{Y,5}$ to zero. 
We model the system with the following kinetic equation: 
$ {d \mu_{Y,5}(t)}/ {dt} = - \Gamma_\mathrm{f}(t) \, \mu_{Y,5}(t)$ , 
where 
$\Gamma_{\rm f}(t)$ is the comoving charge-erasure rate \cite{Bodeker:2019ajh}. This rate corresponds to reactions involving the Higgs that create or destroy right-chiral electrons (or their {CP}-conjugate antimatter partner).  
The strength of these interactions is parameterized by the Yukawa coupling $y_e = \sqrt{2} m_e / v \simeq 3 \times 10^{-6}$ 
%TK what is v?
%EC: higgs vev? (maybe andrew can confirm, doesn't seem to be defined in notes)
and as expected, it is $T(t)$-temperature dependent, that implies \cite{Bodeker:2019ajh}: 
$ \Gamma_{\rm f}(t) \sim (1.3 \times 10^{-2}) y_e^2 T(t) \simeq ( 9.6 \times 10^{-19} {\rm eV}) ( {g_\star}/{106.75})^{-1} $.
Here we normalize the today value of the scale factor to be 1. 
%EC: you mean a_0 or a_\star = 1?
%$g_\star$ is normalized to be 106.75 as assumed at the electroweak energy scale in the SM of particle physics. 
%TK a_0=1
Prior to the electroweak epoch, $\Gamma_{\rm f}(t)$ 
is approximately constant,
%Tk reference
and $\mu_{Y,5} \propto e^{-\Gamma_{\rm f} t}$.  
At early times, the charge-erasure reaction is out of equilibrium, thus $\mu_{Y,5}(t)$ is approximately conserved: 
$e^{- \Gamma_{\rm f} t} \approx 1$ for $\Gamma_{\rm f} t \ll 1$.  
However, once the age of the universe exceeds the interaction time scale, an exponential expression of $\mu_{Y,5}(t)$ will occur.  
We define the fiducial charge erasure conformal time $t_\mathrm{f}$ such that 
$\Gamma_{\rm f}(t_\mathrm{f}) t_{\rm f} = 1$, 
which corresponds to   
$ (T/a)_{\rm f}  \simeq ( 8 \times 10^4 \GeV ) ( {g_\ast}/{106.75})^{-1/2}$. 
%TK Reference 
Here we use 
 the conventional range of 
$g_\star =106.75$  
at the electroweak epoch, and assume that the physical temperature of the primordial plasma at the time when charge erasure occurs to be around $(T/a)_{\rm f} \sim 100$ TeV; i.e. it is not possible to generate a primordial magnetic field from the hypercharge-weighted chiral asymmetry via the CPI after the plasma has cooled below a temperature of $\sim 100$ TeV when using the SM conventional parameters range. 
%
On the other hand, the CPI can be operative at higher temperatures, $T > 100$ TeV, and provide
a source of magnetogenesis \cite{Kamada:2022nyt}, 
as discussed below. 
Moreover, the chiral charge erasure dynamics may be modified under different assumptions about the cosmological expansion
history and particle content of the plasma that is beyond the scope of this letter.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The total chiral asymmetry $\mu_{Y,5}$ value determines the wavenumber $k_{\rm CPI}=2\alpha_{Y} |\mu_{Y,5}| /\pi$ (with the hypercharge fine structure constant $\alpha_Y \sim 0.01$), at which the CPI occurs. 
%TK reference
For the small-scale region, $k>k_{\rm CPI}$, the (hyper)magnetic fields decay exponentially, while for the large-scale regime, $0<k \leq k_{\rm CPI}$, one of the modes is exponentially growing and the chemical potential increases as $\propto e^{kk_{\rm CPI} t/\sigma_Y}$ with $\sigma_{Y}$ -- the comoving hypermagnetic conductivity. The fastest growing mode is characterized by $k=k_{\rm CPI}/2$, and the growth time scale is $\eta_{\rm CPI} = 2 \sigma_Y / k_{\rm CPI}^2 = \pi^2 \sigma_Y / 2 \alpha_Y^2 |\mu_{Y,5}|^2$.  
%TK reference
This means that even if the chemical potential is generated very early, the growth of the (hyper)magnetic field will be ineffective until the age of the universe becomes larger than $t_{\rm CPI}$. 
%TK reference
In a RD universe, this corresponds to the time when the plasma physical temperature 
$(T/a)_{\rm CPI} $ is around 
$(1 \times 10^5 \GeV) 
 ( {\alpha_Y}/{0.01} )^{2} 
 ( {g_\ast}/{106.75})^{-1/2} \times 
 ( {\sigma_Y}/{100 T})^{-1} 
 ([{|\mu_{Y,5}|/T}]/{10^{-3}})^{2} 
$
(that corresponds to $a_{\rm CPI} 
	\simeq (6 \times 10^{-20}) 
	    ( {\alpha_Y}/{0.01})^{-2} 
	( {g_\ast}/{106.75})^{-1/2} \times 
 ( {\sigma_Y}/{100 T}) 
	( [{|\mu_{Y,5}|/T}]/{10^{-3}})^{-2} $).  
	
%GS
% from entropy conservation I got
%$a_{\rm CPI} 
%	\simeq (7.8 \times 10^{-19} a_0) 
%	    ( {\alpha_Y}/{0.01})^{-2} 
%	( {g_\ast}/{106.75})^{1/6} \times %
%	( {\sigma_Y}/{100 T}) 
%	( [{|\mu_{Y,5}|/T}]/{10^{-3}})^{-2} $
% maybe I am crazy
%TK Gutong's expression - must be checked
%AL on Mar 23 -- Guotong, in my newer notes it's aCPI ~ 5 * 10^-19 ... so I think we're in basically agreement ... if you want we can try to sort out the 5 versus 7.8


As we can see the main parameter that governs the CPI dynamics is $|\mu_{Y,5}|/T$. The fiducial value of $|\mu_{Y,5}|/T = 10^{-3}$ was chosen to give $(T/a)_{\rm CPI} = 100$ TeV.  
 %TK reference
%
To ensure self-consistency of our calculations, it's necessary that the processes that would lead to erasure of the hypercharge-weighted chiral chemical potential are out of equilibrium at the time when the CPI develops.  
This condition is expressed as the inequality 
%$\eta_\cpi < \eta_\mathrm{f}$ or equivalently 
$(T/a)_{\rm f} < (T/a)_{\rm CPI}$.  
The system equilibrates when the chemical potential ($\mu_{Y,5}$) is depleted by an $O(1)$ factor, and the helicity density grows by a commensurate factor. This allows us to estimate 
the comoving hypermagnetic helicity density at the time of the chiral plasma instability  
\begin{eqnarray}
	{\mathcal H}_{{\rm M, \rm CPI}} &\simeq &
 \bigl(4 \times 10^{-40}
% {\rm G^2} {\rm Mpc} \, a_0^3 \bigr) 
{\rm G^2} {\rm Mpc}  \bigr) 
 %TK a_0=1
\\
\nonumber
& \times &	
 \biggl( \frac{\alpha_Y}{0.01} \biggr)^{-1} 
 \biggl( \frac{g_\ast}{106.75} \biggr)^{-1} 
 \biggl( \frac{\mu_{Y,5}/T}{10^{-3}} \biggr)   
\label{H-CPI}
\end{eqnarray}
%TK commented eqs with a_0
%\begin{eqnarray}	{\mathcal H}_{{\rm M, \rm CPI}} &\simeq & \bigl(4 \times 10^{-40}  {\rm G^2} {\rm Mpc} \, a_0^3 \bigr) \\ \nonumber & \times &	\biggl( \frac{\alpha_Y}{0.01} \biggr)^{-1} 	\biggl( \frac{g_\ast}{106.75} \biggr)^{-1} \biggl( \frac{\mu_{Y,5}/T}{10^{-3}} \biggr)    \label{H-CPI} \end{eqnarray}
with ${\rm G} \equiv {\rm Gauss}$. 
The corresponding comoving coherence length $\xi_{B, {\rm CPI}}$ and comoving field strength 
$B_{\rm CPI}$ at the time of the CPI are given by 
$\xi_{{B, \rm CPI}} \approx k^{-1}_{\rm CPI}$ and $B_{\rm CPI} \approx \sqrt{{\mathcal H}_{{M, \rm CPI}} / \xi_{{M, \rm CPI}}}$, where the second relation assumes that the field is maximally helical.  
In the presented description the influence of the immersed magnetic field is not accounted for.
%TK maybe present a full expression for equation?
However the presence of primordial (hyper)magnetic field prior the CPI may 
overcome the spin-flip reaction preserving the conservation of total chirality
at temperatures lower than $T_{\rm f}$ (when the characteristic time scale $\Gamma_{\rm f}^{-1}(t)$ is less than the Hubble time scale $H^{-1}(t)$), thus leading to the exponential decrease of $\mu_{Y,5}$). 
Correspondingly the effectiveness of the (hyper)magnetic field generation and its growth depends on the interplay between the chiral asymmetry, temperature, and the spin-flipping modeling \cite{Kamada:2016cnb}. Moreover, in the case when the total chirality is zero (i.e. the fermion chirality is balanced by the helicity of the magnetic field), the evolution of the (hyper)magnetic field follows different patterns, see \cite{Brandenburg:2023aco}. 
Summarizing the dynamics of the CPI induced (hyper)magnetic field might be split into two stages (regime I and II) \cite{Brandenburg:2017rcb}: during the first stage there is an exponential growth of the field that is followed by the inverse cascade of the saturated field. 

After being generated the (hyper)magnetic field, $\B$, along with plasma kinetic motions, experiences MHD decay -- scaling of the (hyper)magnetic energy density ${\mathcal E}_{\rm M}(t) \propto t^{p}$ (and its correlation length $\xi_{\rm M}(t) \propto t^q$) is determined by the initial conditions as well as presence of (hyper)magnetic and kinetic helicity. During the free decay (hyper)magnetic helicity (that is sometimes referred as Hopf invariant 
\cite{Arnol'd:1992} 
in mathematics and corresponds to the linking number in topology \cite{Berger:1983abc}) is a conserved quantity: thus, when the fermion-associated chirality is fully depleted, the magnetic helicity becomes being a conserved quantity 
${\mathcal H}_{\rm M}= \langle \A \cdot \B \rangle ={\rm const}$ (here 
$\A$ is the vector potential, thus 
$\B= {\rm curl} \A$; 
$\langle...\rangle$ denotes ensemble averaging). Moreover, the invariance of (hyper)magnetic helicity ${\mathcal H}_{\rm M}$ leads to 
the inverse cascade phenomena \cite{Frisch:1975} and explains the partial helicity growth (i.e. redistribution of (hyper)magnetic field lines) at large scales until the maximal helical configuration is reached, if the time period allows,  \cite{Tevzadze:2012kk}. In both the fully and partially helical cases the Cauchy–Bunyakovsky–Schwarz inequality 
${\mathcal H}_{\rm M} \leq  2 \xi_{\rm M}(t) {\mathcal E}_{\rm M}(t)$, 
with $\xi_{\rm M}(t)$ the (hyper)magnetic field correlation length, 
is satisfied. The equality holds in the case of fully helical (hyper)magnetic fields ($p+q=0$,
%EC: I have used a different p-q convention in my figure. just making a note here to either change in text or in figure
%TK Axel, please make the choice
see Ref.~\cite{Brandenburg:2016odr} for more details).\footnote{In terms of spectral energy density we have $E_{\rm M}(k)$ (${\mathcal E}_{\rm M} = \int E_{\rm M}(k) dk $) and spectral helicity density $H_{\rm M}(k)$ ((${\mathcal H}_{\rm M} = \int H_{\rm M}(k) dk $)), and,   in addition $\xi_{\rm M} = k_{\rm M}^{-1} =  \int k^{-1} E_{\rm M}(k)/{{\mathcal E}_{\rm M}}$, so correspondingly,  $|H_{\rm M}(k)|\leq 2 k^{-1} E_{\rm M}(k) $. }
Based on the
Cauchy–Bunyakovsky–Schwarz inequality, we characterize fractional helicity through the parameter
$\sigma (t) = {\mathcal H}_{\rm M}/(2\xi_{\rm M} {\mathcal E}_{\rm M}) \leq 1$.\footnote{The situation is 
somewhat 
different if the seed %magnetic 
field is non-helical ($\sigma=0$): the MHD decay still occurs but with a different scaling and with a slower inverse cascading of energy at large scales \cite{Brandenburg:2014mwa}. }

%\paragraph{\textbf{don't know if we need to mention BBN/CMB here -- maybe later when we're talking about observables}}
Big bang nucleosynthesis (BBN) offers various opportunities to test physical processes in the early universe. That includes bounds on additional relativistic degrees of freedom conventionally defined as $\Delta N_{\rm eff} \leq 0.122$ (see Ref.~\cite{Yeh:2022heq} and references therein) 
that can be expressed in terms of additional energy density of radiation -- no more than $\sim 3\%$ can be added to the radiation field (photons) at the moment of BBN in order to satisfy limits from BBN and be consistent with the cosmic microwave background (CMB) {\it Planck} data \cite{Fields:2019pfx}. In the terms of the (hyper)magnetic field, this reads 
a value around $10^{-6}$ G. Obviously the CPI generated magnetic field strength  (XX) is well through these limits. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\paragraph{\textbf{Context 2 -- The preceding discussion is partially thwarted by the the fact that the chiral asymmetry is not conserved, but rather violated by reactions involving the fermion's Dirac mass (or Yukawa interaction) -- nevertheless, at a sufficiently early time in the cosmic history, these chirality-changing reactions would not yet have come into thermal equilibrium -- note that $H \propto T^2$ whereas $\Gamma \propto T$ so $\Gamma \ll H$ at high $T$ -- A careful analysis of the chirality-violating interaction rates reveals that chirality is approximately conserved (in Standard Model processes) in the early universe at temperatures $T \gtrsim 80\,\mathrm{TeV}$ -- at such high temperature, the electroweak symmetry is unbroken, the weak gauge bosons do not yet develop mass through the Higgs field's vacuum expectation value, and it is better to talk about the hypercharge symmetry rather than electromagnetic symmetry and charge -- at plasma temperature $T \gtrsim 80\,\mathrm{TeV}$ the chiral asymmetry is approximately conserved, whereas at lower temperature it is driven to zero exponentially quickly by chirality-violating Standard Model reactions}} asdfasdfasf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\paragraph{\textbf{Setup -- give chiral MHD equations ... mention that everything is comoving ... mention that $\hbar = c = 1$ and Heaviside-Lorentz ... give expressions for all the plasma properties (viscosity, conductivity, spin-flip rate, etc) ... point out which term corresponds to CME ... emphasize that $\mu_{Y,5}/T$ is the only free parameter}}

%\paragraph{\textbf{CPI -- explain how CME term leads to growth of helical magnetic field ... estimate the time when this growth occurs ... estimate the scale of the unstable modes ... estimate the amplitude of the B-field that's generated (assuming maximal helical) ... cite necessary references ... claim that these analytical estimates have been backed up by numerical simulation (right?)}}

%\paragraph{\textbf{Large chiral asymmetry -- emphasize that the CPI must set in before chiral charge erasure occurs at $80\,\mathrm{TeV}$ ... this imposes a lower-limit on $\mu_{Y,5}/T$ which is around $10^{-3}$ ... if the chemical potential is smaller, the asymmetry is erased by scattering and decay processes in the plasma before the magnetic field is amplified, and the resultant gravitational wave amplitude is effectively zero ... we focus on models with a larger chiral asymmetry $\mu_{Y,5}/T \gtrsim 10^{-3}$ ... for references, the baryon asymmetry is $\mu_\textsf{B}/T \sim 10^{-8}$, so we're talking about a large chiral asymmetry}}

%\paragraph{\textbf{GW -- having estimate the generation of primordial magnetic field, we can estimate the corresponding gravitational wave signal -- B-field enters the RHS of Einstein's equation and its anisotropic stress sources GWs ... if all of the GW generation occurs within $O(1)$ Hubble time of the onset of CPI then the frequency of the GWs is set by the CPI scale, which gives $f_\mathrm{CPI} \sim XXX$ ... the amplitude of the GWs depends on how long the source remains active ... if it remains active for $O(1)$ Hubble times, then the GW relic energy today is expected to be $\Omega_\mathrm{gw} \sim XXX$ -- These GWs are very high-frequency (for references LIGO probes gravitational waves around $f \sim 100\,\mathrm{Hz}$) and incredibly weak -- it is not surprising that the GWs are very high-frequency, since they were generated very early in the cosmic history (prior to $T \sim 100\,\mathrm{TeV}$) when the Hubble scale was small ... if the source remains actives for more than $O(1)$ Hubble times, due to the persistent turbulent fluid motions, then the GW amplitude may be larger than the value estimated here -- but we don't expect it to be many orders of magnitude larger}}

%\paragraph{\textbf{Numerical validation -- To validate these estimate, we have performed a numerical simulation of the CPI using the Pencil code ... }}


%Interestingly, the light element abundances (with the bounds on $\neff$) impose limits on the lepton asymmetry in the universe \cite{Barger:2003zg,Abazajian:2004aj,Simha:2008mt,Gelmini:2020ekg}. 





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Since the magnetic field energy density, neglecting the decay, ${\mathcal E}_M$ scales as radiation ($1/a^4$), the magnetic field contributes to the relativistic energy density and thus can be constrained through the BBN bound $\Delta N_{\rm eff}$ implying the effective value of the comoving magnetic field at BBN  $ [B_{\rm eff}^{(\rm BBN)}] = 6\times10^{-7}$ G.  Recently accounting for this BBN bound on $N_{\rm eff}$ and the different decay rate we have determined the upper bounds of the causal  magnetic fields in the RD epoch \cite{Kahniashvili:2021gym},  (see their Fig.~1). Obviously the resulting limits depend on the decay time-scale, and higher temperature magnetogenesis will result in a longer decay period, more effective decrease of the magnetic energy, and more effective growth of the correlation length. 




A separate question is to investigate the gravitational wave signal originated from the CME/CPI efefcts. First, note thaty the gravitational waves will be circularly polarized reflecting the chiral nature of the source \cite{Brandenburg:2021aln}.  Second, as discussed above the CPI implies to the (hyper)magnetic field exponential growth in the relatively short time-scale. As soon as the magnetic field is saturated, the MHD turbulent source leads to the gravitational wave signal generation (through its anisotropic stress projected to be traceless and transverse) with the typical frequency twice the peak frequency (that correspinds to the CPI length of the magnetic source, i.e the comoving linear frequency 
$f_{\rm GW} \equiv \omega/(2\pi) \simeq 2\times k_{\rm CPI}/(2\pi) = 1/\lambda_{\rm CPI}$,
for more details see Ref.~\cite{Brandenburg:2021aln}. 
Using the expression for $\lambda_{\rm CPI}$, we get the typical frequency : 
\begin{equation}
f_{\rm CPI} 	
\simeq 
\big( 10^5 \Hz \big ) 
\bigg ( \frac{\alpha_Y}{0.01} \biggr)
\biggl( \frac{g_\star}{106.75} \biggr)^{-1/3} 
\biggl( \frac{\mu_{Y,5}/T}{10^{-3}} \biggr). %\label{eq:n\u_cpi}
\end{equation}
So this is a $\sim 10 \kHz$ signal.  
For reference, the LIGO interferometer's sensitivity is maximal around $100 \Hz$.  
LIGO's sensitivity is worsening at higher frequencies due to quantum shot noise, but there is still some sensitivity at $10 \kHz$.  
Moreover, there is growing effort to develop technologies that might be able to probe higher-frequency gravitational waves 
\cite{Aggarwal:2020olq}. 
Note, the 
space-based interferometers like LISA are targeting substantially lower frequencies $\sim \mathrm{mHz}$, and not applicable for the CPI operated at high temperatures ($\gg 100 $ GeV).   As far as the signal strength, its smalness (see below) is determined by the low efficiency of the gravitational wave production as well as the weakness of the CPI source itself. 
Since both the gravitational waves and radiation energy densities (as well as the "frozen-in" magnetic field energy density) scales the same with the expansion of the universe, ($\propto a^{-4}$), the ratio between the energy densities ${\mathcal E}_{\rm GW}/\rho_{\rm rad}$ remains constant (the possible damping induced by free-streaming neutrinos operated at frequencies less than $10^{-9}$ Hz is not important in this frequency range, \cite{Weinberg:2003ur}). 
To faciliate the comparison with other sources of relic gravitational waves, we express the gravitational energy density (per logarithmic
unit of frequency) today $\rho_{\rm GW}(w)$ normalized to the critical density ($\rho_{\rm cr}= 8\pi G/(3H^2_0)$ with $H_0$ Hubble parameter today), 
$\Omega_{\rm GW}(f)=(1/\rho_{\rm cr})d{\mathcal E}_{\rm GW}/d{\rm ln}f$ \cite{MaggioreV2}.
There are three ranges of the frequencies as discussed: low frequencies tail, $f<f_H$ -- below 
the Hubble frequency  today, 
\begin{equation}
f_H = \bigl( 1.6 \times 10^{-2}\,\mathrm{Hz} \bigr) \biggl( \frac{T_\star}{100\,\mathrm{TeV}} \biggr) \biggl( \frac{g_\star}{106.75} \biggr)^{1/6}
\nonumber
\end{equation}
Here we normalize temperature to the lowest possible temperature for CPI, $100\,\mathrm{TeV}$; the intermediate frequencies, $f_H \leq f \leq f_{\rm CPI}$, and the high frequencies tail, $f >f_H$. Outside the horizon ($f < f_H$) the causality requirement leads to 
the $\Omega_{GW}(f) \propto f^3$ scaling \cite{MaggioreV2}. 
This can be seen also through the causality and dimension arguments: outside the horizon the gravitational waves cannot be affected by the physical processes inside the horizon, the two point correlation function for the gravitational wave strength $\langle h_{ij} ({\bf x} t_0) h_{lm}({\bf x+r}, t_0  )\rangle \rightarrow 0$ vanishes for the space separation larger than the Hubble horizon, $r=|{\bf r}| \leq H_)^{-1}$ leading to the $f^3$ for the gravitational wave energy density (per logarithmic frequency)  energy, $\rho(f) \propto f^3$ for $f<f_H$. This scaling is universal for any causally generated gravitational waves (with the corresponding $f_H$). It can be also derived through air-acoustic approximation (a direct analogy with the sound waves generation by turbulence \cite{Lighthill:1952,Proudman:1952}) used in analytical computation of Ref.~\cite{Gogoberidze:2007an}. Contrary, the -
intermediate and high frequency spectral shapes are determined by source characteristics as discussed below, and differs depending on the nature what process immerse the gravitational radiation signal. In particular, 
dimensional analysis of the CPI process, see Ref.~\cite{RoperPol:2019wvy}, at high frequencies $f \gg f_\star$ (where $f_\star$ is a comoving linear frequency that corresponds to the twice so called "inertial" wavenumber)  the energy density spectra of the source and the gravitational background are shifted by a factor of $k^{-2}$ implying for Kolmogorov spectrum (${-2/3}$ exponent) $\omega_{\rm GW} (f) h^2 \propto f^{-8/3} $ (with $h$ Hubbke parameter  today normalized to 100km$/$sec$/$Mpc. 
To validate these estimates, we have performed a numerical simulation of the CPI using the Pencil code. 
%TK Maybe setup for sims?
%TK reference for Pencil code
In the intermediate range, $f_H <f < f_\star$ the proper scalling of the spectrum was obtained for the first time by some of authors, \cite{RoperPol:2018sap,RoperPol:2019wvy} for the waves sourced by turbulence and magnetic fields present at the electroweak energy scale and shown $\Omega_{\rm GW}(f) h^2 \propto f^{1}$ \cite{Brandenburg:2021aln}.  For the CME$/$CPI generated gravitational waves, the scalling shown the substantially different  spectrum 
$\Omega_{\rm GW}(f) h^2 \propto f^{-2}$.  Because of the the smallness of the CPI length scale it is impossible to resolve the spectrum, 
so we made extrapolation of $f^{-2}$ till the Hubble frequency, $f_H$.
%Emma, please add about extrapolation. Maybe the figure as well? or ask Axel?

%Guotong, please add your figure. 
%TK I will discuss the dependence of efficiency on viscosity, etc...

Gravitational waves efficiency...

%TK dependence on parameters
\vspace{2mm}
{\bf 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{2mm}
{\bf Acknowledgements}----We thank Kohei Kamada, Sayan Mandal, and Jonathan Stepp for useful discussions and comments. Support through the NASA ATP award 80NSSC22K0825, 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. E.C. acknowledges the Pake Fellowship, and G.S. acknowledges 
support from the Undergraduate Research Office in the
form of a Summer Undergraduate Research Fellowships at Carnegie Mellon University.

%TK GWs from CME
%TK https://arxiv.org/pdf/2101.08178.pdf

%TK https://arxiv.org/pdf/1809.04501.pdf
%TK N_eff from GWs

%\bibliographystyle{JHEP}
\bibliographystyle{h-physrev5}
\bibliography{Chirality-BBN}



 



\end{document}

%EC: placing figure here for now
\begin{figure*}[t]
    \centering
    %\includegraphics[scale=0.55]{fig/Bev-wAck.eps}
    \includegraphics[scale=0.8]{fig/field-evolution-scaled.eps}
    \caption{Turbulent evolution of the comoving magnetic field strength and correlation length normalized by initial temperature and relativistic degrees of freedom. The different combinations of the scaling exponents $p$ and $q$ correspond to different helical properties. Factors of $g$ have been ignored. Normalized to EW scale $T_{\rm EW} = 100$GeV. Gray shaded region corresponds to bounds from Fermi \cite{Fermi-LAT:2018jdy} and are on a scale of log$_{10}(B/{\rm Gauss})$ and log$_{10}(\xi/{\rm Mpc})$ and are independent of $T_*$, $g_*$, and the scaling exponents.}
    \label{fig:fieldevolution}
\end{figure*}

\begin{figure*}[t]
    \centering
    \includegraphics[scale=0.8]{fig/pspecm_wkf6_600_scaled.eps}
    \caption{Old figure rescaled.}
    \label{fig:pspecm_wkf6_600_scaled}
\end{figure*}

Our major goal is to estimate the GW signal keeping the possibility of varying the model parameters, $\mu_{Y, 5}$, $T_\star$ and $g_\star$. Besides these two parameters we should keep in mind that the strength of the GW source might be also parametrized. 
%If the CME limits the total magnetic helicity, the strength of the magnetic field is determined through additional parameter that is defined through the ratio of the magnetic field correlation length $\xi_{\rm M}$ and the Hubble length scale at the moment of CME, $\gamma \equiv \xi_{\rm M}/H_\star^{-1}$ (if not specified explicitly we refer to comoving quantities). The BBN and CMB data impose additional constraints on the total magnetic energy density \cite{Kahniashvili:2021gym}. 

On the other hand, we should to account that the peak frequency will lie in much higher frequency range: in particular, the GW spectrum plateau extend till the Hubble frequency, i.e. 
%TK I put EWPT scales - Emma recalculate/check
%EC: I agree with the frequency expression
\begin{equation}
f_{H_\star} = 1.65 \times 10^{-5} {\rm Hz} \Big( \frac{T_\star}{100{\rm GeV}}\Big ) \Big( \frac{g_\star}{100}\Big )^{\frac{1}{6}},
\end{equation}
while the peak frequency will be $f_\star = 2 f_{H_\star} \gamma^{-1}$. 
%TK Emma, please recalculate 
In the standard scenarios of the electroweak phase transitions it is assumed that $\gamma^{-1} \simeq 100$, however in the case of the CME $\gamma^{-1}$ might be substantially higher. Accounting 
$T_\star \approx 100 {\rm TeV}$ we might reach the detection prospects by LIGO/VIRGO. 
%EC: even 10 TeV would reach LIGO (I mean lowest frequency 10Hz) if we are assuming gamma^-1=100
%TK Axel, what will be the \Omega_{\rm GW}

In fact there is a possibility to consider CME at higher energies ($T \gg T_{\rm CPI}$ 
\cite{Domcke:2022uue}. In this case the stronger GW signal can be obtained with a shifted range of the detection prospects. 


including photons and relativistic fermions, 
and assuming that the initial chirality sourced by 
the magnetic field via the CME MHD instability,
the maximal magnetic helicity 
%TK measured today - we probably do not need quote today because it is constant 
is \cite{Brandenburg:2017rcb}, 
\begin{equation}
B^2 \xi_{\rm M}=%\frac{1}{4\alpha_{\rm em}} \frac{g_0}{g_\star}n_{\gamma 0} N_{\rm f}
%TK Axel, I am confused about the g dependence
%TK B^2 \xi_M = const
\label{chirality}
\end{equation}
that is substantially smaller then ${\mathcal H}_{\rm M, CPI}$. 
%TK where  $g_0=3.36$ is the number of relativistic number of freedom today, and $n_{\gamma 0}= 2\zeta(3)/\pi^2 T_0^3$ (with $\zeta(3)=1.202$ Riemann Zeta function) is the number of photons today ($411{\rm cm}^{-3}$) with today-temperature, $T_0 \approx 2.73{\rm K}$. 
%however the question of baryon asymmetry is beyond of the current work scopes and is deferred to future studies. 
%TK more here

{It is convenient to parametrize the MHD turbulence strength through the corresponding comoving Alfv\'en velocity $v_A$, 
%$v_A^2 = 1/{1+(\rho+P)/{(2\mathcal{E}_{\rm M}})}$
$v_A^2 = 1/(1+(\rho+P)/{2\mathcal{E}_{\rm M}})$ 
\cite{Gedalin:1993} (with all quantities are rescalled to be comoving) that in the non-relativistic limit leads to 
the well-known expression
$v_A^2 = B_{\rm eff}^2/{4\pi (\rho+P)}={2{\mathcal{E}_M}/(\rho+P)}$. To avoid spoil the 
expansion during the RD epoch, we assume that $\chi_{\rm M} \equiv {\mathcal E}_{\rm M} / \rho_{\rm rad} \leq 1$, and we consider $\chi_{\rm M}$ as a free parameter of our model, $v_A^2 = 3\chi_{\rm M}/(3\chi_{\rm M}+2)$ (equal $1$ if $\chi_{\rm M} \rightarrow \infty$ and $3\chi_{\rm M}/2$ if $\chi_{\rm M} \ll 1$).
}
%GS: I am confused about the above limits. Didn't we assume that \chi_{\rm M} is always less than one?
%TK In general chi might be larger than 1
We also emphasize that accounting the decay makes $\chi_{\rm M}$ time dependent, in particular, $\chi_{\rm M} \propto t^{p/2}$, 
, and thus starting with the maximal allowed $\chi_{\rm M,0}=1$, 
at the electroweak scale, $T_{\rm EW}$ (and with $g_{\rm EW}$), the $\chi_{\rm M} = (T_{\rm EW}/T_\star)^{-p} (g_{\rm EW}/g_\star)^{-p/3}$ and the correlation will reach $\xi_{\rm M} (t_{\rm EW}) = \gamma H_\star^{-1} (T_{\rm EW}/T_\star)^{-q} (g_{\rm EW}/g_\star)^{-q/3}$. For fully helical case, $p+q=0$ ($p=-2/3$) and ${\mathcal E}_{\rm M}\xi_{\rm M} = {\rm const} = \mu_0$. There are important consequences: i) if the CME is preferable at high temperatures with higher values of $N_{\rm f}$ (i.e. $g_\star$), the magnetic field strength that satisfies the BBN and CMB bounds will be weaker at the electroweak scale by the factor of $(T_{\rm EW}/T_\star)^{-1/3} (g_{\rm EW}/g_\star)^{-1/9}$; the dependence on $g_\star$ is somehow weak, however using the CME scale at the $100$TeV 
the magnetic field amplitude will be 10 times weaker. ii) The {\it maximal} correlation length ($\gamma=1$ that is unrealistically high because CME operated at micro-physics scales)
%TK reference
will be also smaller (compared to that when CME occurs at the electroweak energy scale) because of the Hubble horizon $H_\star^{-1}$ scales as $Tg^{1/6}$ while the maximal growth rate is proportional to $T^{-2/3} g^{-2/9}$. The overall dependence will be $T^{1/3}g^{-1/18}$. The extremely week dependence on $g$ allows us to discard it, being focused only on the temperature dependence. If the CME at the electroweak energy does not insures production of viable primordial magnetic fields that satisfy the lower bounds imposed by Fermi data
%TK reference
the high-energy CME with higher numbers of chiral fermions has even more difficulties to reach out necessary values: This discussion is summarized on Figure 1.
%TK reference


The GW spectrum gets saturated in relatively short time scales (determined by a largest eddy turnover-time scale) due to the source decay effect, and because of the low efficiency of the magnetic field conversion into GWs, the gravitational radiation  does not contribute significantly to the total energy density of radiation. 
%TK reference



  
 

