% $Id: paper.tex,v 1.59 2021/01/15 06:19:17 roper Exp $
\documentclass[a4paper,11pt]{article}
%AB: ApJ style
%\documentclass[preprint2,numberedappendix]{emulateapj-rtx4}
%\pdfoutput=1
\usepackage{jcappub}
%AB: outcommented
%\usepackage[T1]{fontenc} % if needed

\usepackage{amssymb}
\usepackage{graphicx,bm,color}
\graphicspath{{./fig/}{./png/}}

%|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
%             Customized Commands
%|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
%  mathematical abbreviations
%  =========================
%
% 12-feb-2003: wd introduced \BoldVec macro to get size of vectors in
%              exponents right
%
\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}$}}{}}{}
% 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}$} {}}
%
\newcommand{\GGGG}{{\bf G} {}}
%
%  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\EEK{{\cal E}_{\rm K}}
\def\EEKM{{\cal E}_{\rm K/M}}
\def\EEM{{\cal E}_{\rm M}}
\def\EEGW{{\cal E}_{\rm GW}}
\def\OmGW{{\Omega}_{\rm GW}}
\def\XiGW{{\Xi}_{\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\Sh{S_h}
\def\Shd{S_{\dot{h}}}
\def\Ah{A_h}
\def\Ahd{A_{\dot{h}}}
\def\cs{c_{\rm s}}
\def\xiM{\xi_{\rm M}}
\def\xiK{\xi_{\rm K}}
\def\kf{k_*}
\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}}
%
\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{\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}}
\newcommand{\GeV}{\,{\rm GeV}}
\newcommand{\MeV}{\,{\rm MeV}}
%
\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}}
\def\tmax{t_{\rm max}}
%
%  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{\sprl}[2]{, Phys.\ Rev.\ Lett.\, submitted, arXiv:#2 (#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{\sprd}[2]{, Phys. Rev. D, 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{Polarization of gravitational waves from helical MHD turbulent sources}

%\date{Received 5 July 2016; published 3 February 2017,~ $ $Revision: 1.59 $ $}
%\preprint{NORDITA-2016-82 --- Phys. Rev. Lett. 118, 055102 (2017)}
\date{\today,~ $ $Revision: 1.59 $ $}
%\date{\!$ \, $Revision: 1.59 $ $\!}
%\preprint{NORDITA-2018-xx}

\author[a,b,1]{Alberto~Roper~Pol}
\emailAdd{roperpol@apc.in2p3.fr; $ $Revision: 1.59 $ $}
\affiliation[a]{Universit\'e de Paris, CNRS, Astrophysique et Cosmologie, Paris, F-75013, France}
\affiliation[b]{Faculty of Natural Sciences and Medicine, Ilia State University,
3-5 Cholokashvili Ave, Tbilisi, GE-0194, Georgia}

\author[c,b]{Sayan~Mandal}
%\email{sayan.mandal@stonybrook.edu}
\affiliation[c]{Physics and Astronomy Department, Stony Brook University, Stony Brook,
New York 11794, USA}

\author[d,e,f,b]{Axel Brandenburg}
%\email{brandenb@nordita.org}
%\affiliation{Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, CO 80303, USA}
\affiliation[d]{Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden}
\affiliation[e]{Department of Astronomy, AlbaNova University Center, Stockholm University, 10691 Stockholm, Sweden}
\affiliation[f]{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, PA 15213, USA}

\author[f,b,g,h]{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[g]{Abastumani Astrophysical Observatory, Tbilisi, GE-0179, Georgia}
\affiliation[h]{Department of Physics, Laurentian University, Ramsey Lake Road, Sudbury, ON P3E 2C, Canada}

%\author[i]{Arthur~Kosowsky}
%\email{kosowsky@pitt.edu}
%\affiliation[i]{Department of Physics and Astronomy, University of Pittsburgh, and
%Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT PACC), Pittsburgh PA 15260}

\abstract{
%\begin{abstract}
We use direct numerical simulations of decaying primordial hydromagnetic
turbulence with helicity to compute the resulting production of
gravitational waves (GWs) and their degree of polarization.
%The magnetic field is either given initially, or it is produced by
%AR: added text to include velocity fields too
The turbulence sourcing field, which can be magnetic or kinetic,
is either given initially, or it is produced by
%AR
driving it with an electromotive force applied during a short initial
time interval.
In both types of simulations, we find a clear dependence of the degree
%of circular polarization of the resulting GWs on the fractional magnetic
of circular polarization of the resulting GWs on the fractional magnetic/kinetic
helicity.
Furthermore, the spectral peak is higher in the helical case, which
facilitates its observational detection with the planned Laser
Interferometer Space Antenna (LISA).
In both cases, the low frequency tail is shallower than what was expected
based on earlier analytical work.
For frequencies above the peak, the two types of simulations show opposite
trends in that the spectral energy is higher (lower) in the helical case
than in the nonhelical one if an initial spectrum is given (driven).
The degree of polarization approaches zero at frequencies below the peak,
but stays finite at higher frequencies if the initial magnetic field is
given to be fully helical all wavenumbers, and it approaches zero if the
initial field is driven by a fully helical monochromatic forcing function.
}
%\pacs{98.70.Vc, 98.80.-k}
%\end{abstract}

%\keywords{
%magnetic fields --- MHD --- turbulence --- polarization
%}

\begin{document}
\maketitle
\flushbottom

\section{Introduction}

Primordial magnetic fields produced during phase transitions can drive stochastic
%background gravitational waves (GWs) in the milliHertz range and could in
%AR: made more general than mHz which only considers EWPT and k_* = 100 
background gravitational waves (GWs) in a broad range of frequencies,
from nHz to Hz, depending on the energy scale of the phase transition
($T \sim 100 \GeV$ for the electroweak, and $T \sim 100
\MeV$ for the QCD phase transition), and the characteristic scale of the
magnetic field $\kf$, which is taken to be a fraction of the Hubble scale
at the time of magnetic field generation.
These signals could in
%AR
%future be observable with the planned Laser Interferometer Space Antenna
%(LISA).
%AR: added PTA
future be observable with gravitational wave detectors, e.g., pulsar timing
arrays (PTA), or the planned Laser Interferometer Space Antenna
(LISA).
%AR: more can be commented on this for different detectors
These magnetic fields could have survived until the present
time \cite{BKMRPTV17} and are thus a leading candidate for explaining the
lower limits on the strength of intergalactic magnetic fields inferred
from the non-observation of cascade photos from TeV photons of blazars
interacting with the cosmic background light \cite{NV10}.
Such fields are expected to be turbulent and partially helical.

Helical magnetic fields produce circularly polarized GWs \cite{KGR05}.
However, the degree of polarization of GWs has been a matter of
uncertainty owing to the approximations made in the analytical
calculations available to date.
The work of Ref.~\cite{KK15} suggested that the circular polarization of
GWs would be about 80\% for a maximally helical magnetic field, by
assuming Kolmogorov-type turbulence, with spectral index of $-11/3$ for
the magnetic power spectrum and $-14/3$ for the helicity power spectrum,
following Ref.~\cite{K65}.
The peak is observed at twice the magnetic spectral peak.
For smaller values of the fractional magnetic helicity, the maximum
degree of polarization also diminishes.
Ref.~\cite{KGR05} showed that the maximum circular polarization depends
%on the spectral index and can actually reach nearly 100\%, but only for
%a Moiseev-type spectrum, with equal spectral indices for the magnetic
%field and the magnetic helicity, following Ref.~\cite{MC96}.
%AR: slightyl changed
on the relation between the magnetic field and magnetic helicity power spectra,
and it can actually reach nearly 100\% in the case
when the spectral indices for both are equal, which we call a Moiseev-type spectrum,
following Ref.~\cite{MC96}.
On the other hand, this type of spectrum generates a degree of polarization
that does not decay to zero for large wave numbers, but it stays constant
and is approximately equal to the fractional helicity.
%AR: added
More recently, the two types of spectra mentioned above have been studied
in Ref.~\cite{EFLVW20} in the context of detectability with LISA.
We show that the circular degree of polarization computed from direct
numerical simulations does not follow any of the previous analytical models,
even though it retains some similarities, depending on the how the magnetic
field is generated/driven at the initial time.

%In two recent papers \cite{RPBKKM18,RPMBKK19}, we have described the
%AR: modified to 'the authors'
In two recent papers \cite{RPBKKM18,RPMBKK19}, the authors have described the
implementation of a GW solver into the {\sc Pencil Code} \cite{PC} and
have presented numerical computations for primordial hydrodynamic and
%hydromagnetic turbulence from phase transitions.
%AR: added
hydromagnetic turbulence from phase transitions,
and their subsequent generation of a stochastic gravitational wave
background, also computed numerically.
%AR
Some of those solutions were fully helical, but helicity was not a focus
of their studies.
%AR: added reference to PRL submitted paper
In Ref.~\cite{KBGMRP20} the degree of circular polarization for
kinetically and magnetically forced turbulence was presented.
We extend the study in the present work by including the effect of
fractional helicity to the magnetic field power spectrum, and the variation
of the polarization for different assumptions of magnetic field production.
%AR

We present here two types of simulations, similar to the two types of
hydromagnetic simulations presented in Ref.~\cite{RPMBKK19}: one where
a primordial magnetic field was assumed given as initial condition,
and one where a magnetic field was generated by an electromotive force
$\emf(\xx,t)$ that depends on time $t$ and position $\xx$.
This driving was then applied during a short time interval,
and then switched off, i.e, set to zero, $\emf=\bm{0}$, such that
turbulence decays for later times.
By using suitably scaled variables and conformal time, the governing
equations describing the evolution of GWs and turbulent magnetic fields
in an expanding universe in the radiation era can be brought into a
form that is best suited for numerical simulations \cite{RPBKKM18}.

We begin by summarizing our approach and the equations solved.
We then present the magnetic and GW energy spectra, the degree of
polarization in relation to the magnetic helicity, and compare with the
sensitivity limits for LISA.
We conclude by highlighting unique properties of the GW spectra and
their polarization properties.

Throughout this work, electromagnetic quantities are expressed in
Lorentz--Heaviside units where the vacuum permeability is unity.
Einstein index notation is used so summation is assumed over repeated
indices.
Latin indices $i$ and $j$ refer to spatial coordinates 1 to 3.
The linear polarization modes are indicated by $\lambda = +, \times$.

\section{The model}

We use the linear polarization modes $+, \times$ to describe
the two gauge independent components of the tensor-mode perturbations,
$\tilde{h}_{ij}(\kk) = \tilde{h}_+(\kk) e_{ij}^+(\kk) +
\tilde{h}_\times(\kk) e_{ij}^\times (\kk)$,
where the tilde indicates that this decomposition is performed
in Fourier space.
The linear polarization basis tensors are
\begin{equation}
e_{ij}^+ (\kk) = e_i^1 e_j^1 - e_i^2 e_j^2,
\quad
e_{ij}^\times (\kk) = e_i^1 e_j^2 + e_i^2 e_j^1,
\end{equation}
where $\ee^1$ and $\ee^2$ form a basis with the units vectors $\kk/|\kk|$
\cite{RPBKKM18}.
We solve the non-dimensional
GW equation in the radiation era for the scaled strains
$\tilde{h}_\lambda$, where $\lambda = +, \times$,
using conformal time,
normalized to unity at the initial time of magnetic energy generation $t_\ast$,
and comoving wave vector, normalized by the corresponding scale
%$c t_\ast$
%AR: changed to
$1/(c t_\ast)$
\cite{RPBKKM18,RPMBKK19}
\begin{equation}
\left( \partial^2_t + \kk^2\right)
\tilde{h}_\lambda (\kk, t)\, =\, {6 \over t}
\tilde{T}_\lambda^{\rm TT} (\kk, t),
\label{GW}
\end{equation}
where $\tilde{T}_\lambda^{\rm TT}$ is the comoving stress tensor, projected
into the traceless and transverse (TT) gauge, described by
the linear polarization modes, and normalized by the
energy density at $t_\ast$.
%The scale factor evolves linearly with conformal time, i.e.,
%$a = t$, using the convention $a_\ast = 1$.
%This leads to $t_\ast = H_\ast^{-1}$, where
%$H = a^{-2} \partial_t a$ is the Hubble rate.
%AR: omitted details
The scaled strains are tensor-mode perturbations over
the Friedmann-Lema\^itre-Robertson-Walker metric tensor.
%Hence, the metric tensor is $a^2 \delta_{ij} + a h_{ij}$
%AR: modified
Hence, the metric tensor is such that the line element is
$\dd s^2 = a^2(-\dd t^2 + \delta_{ij}
+ h_{ij}/a)$
%AR
\cite{RPBKKM18,RPMBKK19}.
%AR: added
During radiation-dominated epoch, the equation of state is
$p = \fourthird \rho$, which leads to a linear evolution of
the scale factor $a$ with $t$.
%AR

The stress is composed of magnetic and kinetic contributions.
In physical space it is
\begin{equation}
T_{ij} (\xx) =\frac{4}{3}\frac{\rho u_i u_j}{1 - \uu^2} - B_i B_j +
\left(\frac{\rho}{3} + \frac{\BB^2}{2}\right)
\delta_{ij},
\label{Tij}
\end{equation}
where $\uu$ is the plasma velocity,
$\rho$ is the energy density, and
$\BB$ is the magnetic field.
%We have used an ultrarelativistic equation of state, so the pressure
%is $p = \rho/3$ for the plasma such that the total enthalpy is
%$w = p + \rho = \fourthird \rho$.
%AR: added above now
The total enthalpy is $w = p + \rho = \fourthird \rho$.
Since $T_{ij}$ refers to comoving and normalized stress tensor, the
MHD fields $\rho, \uu, \BB$ are accordingly normalized and
comoving.
%AR

%We refer to magnetic fields normalized by $\sqrt{\EErad^\ast}$
%and comoving with $a^2$.
%The energy density is normalized by $\EErad^\ast/c^2$,
%and the pressure (and the enthalpy) are normalized by $\EErad^\ast$.
%The energy density, the pressure and the enthalpy are comoving with $a^4$.
%This is consistent with normalizing the stress by the radiation
%energy density, $\EErad^\ast = 3H_\ast^2 c^2/(8 \pi G)$,
%at $t_\ast$, where $c$ is the speed of light and $G$ is
%Newton's constant.
%The stress is also comoving with $a^4$.
%The velocity is normalized with $c$.
%AR: this can be reduced and simplified by avoiding giving so many
%AR: details since we refer to our previous papers, see now text
%AR: above
%
The non-dimensional and comoving MHD equations for an ultrarelativistic
gas in a flat expanding universe in the radiation-dominated
era after the electroweak phase transition are given by
\cite{RPBKKM18,RPMBKK19,BEO96,BKMRPTV17,DN13}
\begin{eqnarray}
{\partial\ln\rho\over\partial t}
&=&-\frac{4}{3}\left(\nab\cdot\uu+\uu\cdot\nab\ln\rho\right)
+{1\over\rho}\left[\uu\cdot(\JJ\times\BB)+\eta \JJ^2\right],
%\label{dlnrhodt}\\
\nonumber \\
{\partial\uu\over\partial t}&=&-\uu\cdot\nab\uu
+{\uu\over3}\left(\nab\cdot\uu+\uu\cdot\nab\ln\rho\right)
%-{\uu\over\rho}\left[\uu\cdot(\JJ\times\BB)+\eta \JJ^2\right]\nonumber \\
%&&-{1\over4}\nab\ln\rho+{3\over4\rho}\JJ\times\BB+{2\over\rho}\nab\cdot\left(\rho\nu\SSSS\right),
+{2\over\rho}\nab\cdot\left(\rho\nu\SSSS\right) \nonumber \\
&&-{1\over4}\nab\ln\rho -{\uu\over\rho}\left[\uu\cdot(\JJ\times\BB)+\eta\JJ^2\right] +{3\over4\rho}\JJ\times\BB,
%\label{dudt} \\
\nonumber \\
{\partial\BB\over\partial t}&=&\nab\times(\uu\times\BB-\eta\JJ+\emf),\quad
\JJ=\nab\times\BB,
%\label{dAdt}
\nonumber
\end{eqnarray}
where ${\sf S}_{ij}=\half(u_{i,j}+u_{j,i})-\onethird\delta_{ij}\nab\cdot\uu$
are the components of the rate-of-strain tensor with commas denoting
partial derivatives, $\JJ$ is the current density,
$\nu$ is the kinematic viscosity, and $\eta$ is the
magnetic diffusivity.
%Similarly to the other fields, the current density
%is comoving with $a^3$ and
%normalized by $H_\ast \sqrt{\EErad^\ast}/c$.
%The kinematic viscosity and the magnetic diffusivity
%are comoving with $a^{-1}$ and normalized by $c^2/H_\ast$.
%AR: not necessary

The electromotive force, $\emf$, is used
to model the magnetic field generation with fractional magnetic helicity.
%AR: added some text
In some of our runs we will omit this term, and we will start with a fully
developed turbulent spectrum of the magnetic field as our initial condition.
On the other hand, we also consider magnetic fields generated dynamically
using $\emf$.
%AR
%Specifically, we write
%AR: modified
Specifically, in the latter case, we write
\begin{equation}
\emf(\xx,t)={\rm Re}\{{\cal N}\tilde{\ff}(\kk(t))\exp[i\kk\cdot\xx+i\varphi]\},
\end{equation}
where the wave vector $\kk(t)$ and the phase $\varphi(t)$ change randomly
from one time step to the next.
This forcing function is therefore white noise in time and consists of plane
waves with average wavenumber $k_*$ such that $|\kk|$ lies in an interval
$k_*-\delta k/2\leq|\kk|<k_*+\delta k/2$ of width $\delta k$.
%Here, ${\cal N}$ is a normalization factor and $\delta t$ is the time step.
%AR: no need to define \delta t
Here, ${\cal N}$ is a normalization factor.
%The Fourier amplitudes are
%AR: added 'of the forcing'
The Fourier amplitudes of the forcing are
\begin{equation}
\tilde{f}_i=\left(\delta_{ij}-i\sigma\epsilon_{ijl} k_l/k\right)
\tilde{f}_j^{(0)}/\sqrt{1+\sigma^2},\;
\end{equation}
where $\tilde{\ff}^{\rm(0)}({\kk})=(\kk\times\ee)/[\kk^2-(\kk\cdot\ee)^2]^{1/2}$
is a nonhelical forcing function.
Here, $\ee$ is an arbitrary unit vector that is not aligned with $\kk$.
Note that $|\ff|^2=1$.
The values $\sigma=0$ and $\sigma=1$ correspond to nonhelical and maximally helical cases,
%but we also consider intermediate values.
%AR: added
but we also consider intermediate values, which correspond to
fractional helicities.
%The forcing is only enabled during the time interval $1\leq t\leq 1.1$.
%AR: slightly changed
The forcing is only enabled during an arbitrarily short time interval
$1\leq t\leq \tmax$.
We chose $\tmax = 1.1$.
%The resulting fractional magnetic helicity is giving approximately by
%$2\sigma/(1+\sigma^2)$.
%AR: added definition for clarity
The resulting fractional magnetic helicity $\kf \bra{\AAA \cdot \BB}/
\bra{\BB^2}$ is given approximately by
$2\sigma/(1+\sigma^2)$, being $\AAA$ the magnetic vector potential,
and angle brackets denoting ensemble average over stochastic realizations.
%AR

When the magnetic field is given as initial condition, we generate
a random three-dimensional vector field in Fourier space,
\begin{equation}
B_i({\kk})=B_0\left(P_{ij}(\kk)-\ii\sigma\epsilon_{ijl}
k_l/k\right) g_j({\kk})\, g_0(k),
\label{Bikk}
\end{equation}
where $g_j(\kk)$ is the Fourier transform of a $\delta$-correlated
%vector field in three dimensions with Gaussian fluctuations, $k_*$
%AR: added info about g_j
vector field in three dimensions with Gaussian fluctuations,
i.e., $g_i(\xx) g_j(\xx') = \delta_{ij} \delta^3(\xx - \xx')$, $k_*$
is now identified with the initial wavenumber of the energy-carrying
%eddies, and $S(k)$ determines the spectral shape with slope $\alpha/2$
%AB: S(k) -> g_0(k)
%eddies, and $S(k)$ determines the spectral shape with slope $\alpha/2$
%AR: changed to g_0, and modified such that details are given below now
%in the low wave number limit and Kolmogorov-type spectral slope  $-5/3$
%in the high wave number range \cite{BKMRPTV17}.
eddies, and $g_0(k)$ determines the spectral shape \cite{BKMRPTV17}.
\begin{equation}
g_0(k)={k_*^{-3/2} (k/k_*)^{\alpha/2-1}
\over[1+(k/k_*)^{2(\alpha+5/3)}]^{1/4}}.
\label{Sfunction}
\end{equation}
%AR: added text and details moved here
The power spectrum of magnetic fields is computed from $\half \bra{\BB^2}$.
Hence, it is proportional to $k^2 g_0^2 (k)$ in Fourier space, which
correspond to a spectral index $\alpha$ in the low wave number limit
(subinertial range) and Kolmogorov-type spectral slope $-5/3$ in the high
wave number range.
%
Here, the $k_*^{-3/2}$ prefactor ensures that the resulting magnetic
energy is independent of the value of $k_*$.

%AR: added definition of spectra for magnetic field
The equal time correlation function of the magnetic field,
assuming statistical homogeneity and isotropy, and a
Gaussian-distribution in space, is
\begin{equation}
\bra{B_i^*(\kk, t) B_j(\kk', t)} = (2\pi)^6 \delta^3
(\kk - \kk') \biggl[P_{ij} \frac{\EM (k, t)}{4 \pi k^2} +
i \epsilon_{ijl} k_l \frac{\HM (k, t)}{8 \pi k^2} \biggr],
\label{equaltime}
\end{equation}
where $P_{ij} = \delta_{ij} - k_i k_j/k^2$, and
$\EM (k, t)$, $\HM (k, t)$ are the magnetic and helical
power spectra, respectively.
%
%The GW and magnetic energies are given by
%$\EEGW=\bra{\dot{h}_{+,\rm phys}^2+\dot{h}_{\times,\rm phys}^2}\,c^2/32\pi G$
%AB: is the factor ok? How should we handle this?
%AR: added details on magnetic field above, and details on GW below
%and $\EEM=\bra{\BB^2}/2$.
%We work here with magnetic and GW energy spectra per linear wave number
%interval, and in the case of GW spectra, we also use GW energy spectra
%per logarithmic wave number interval, $\Omega_{\rm GW}(k)=k\EEGW(k)$.
%AR: only magnetic since text is modified now
We work here with magnetic energy spectra per linear wave number
interval.

%AR: added text
The GW energy density (in physical units) is
\begin{equation}
\EEGW(t) = \frac{c^2}{32 \pi G} \bra{\dot{h}_{ij}^{\rm phys}
\dot{h}_{ij}^{\rm phys}},	
\end{equation}
where $h_{ij}^{\rm phys} = h_{ij}/a$ are the physical strains,
and a dot represents derivative with respect to physical time.
In terms of the normalized and comoving units used in \Eq{GW},
the ratio of comoving GW energy density to critical energy
density $\OmGW = \EEGW/\EEcrit^0$ is
\begin{equation}
a^4 \OmGW(t) = \frac{1}{12} \biggl( \frac{H_*}{H_0} \biggr)^2
\biggl\langle \frac{\partial h_{ij}}{\partial t}
\frac{\partial h_{ij}}{\partial t} + \frac{1}{t^2} h_{ij} h_{ij} -
\frac{2}{t} h_{ij} \frac{\partial h_{ij}}{\partial t}\biggr\rangle,
\end{equation}
where $\EEcrit^0 = 3 H_0^2 c^2/(8 \pi G)$, with
$H_0 = 100\,h_0 \km \s^{-1} \Mpc^{-1} \approx 3.241
\times 10^{-18}\,h_0 \s^{-1}$ being the Hubble rate
at the present time, and $h_0$ takes into account the
uncertainties in its exact value.
To simplify the notation, we will use prime to denote
derivative with respect to $t$.

The equal time correlation function for the strains and
their derivatives (again assuming the statistical fields to
be homogeneous, isotropic and Gaussian) can be expressed as
\begin{align}
\bra{h'_{ij} (\kk, t) h'_{ij} (\kk', t) } = & \, \frac{1}{4}
\bigl[{\cal M}_{ijlm} S_{h'} (k, t) +
i {\cal A}_{ijlm} A_{h'} (k, t) \bigr], \nonumber \\
\bra{h_{ij} (\kk, t) h_{ij} (\kk', t) } = & \, \frac{1}{4}
\bigl[{\cal M}_{ijlm} S_{h} (k, t) +
i {\cal A}_{ijlm} A_{h} (k, t) \bigr], \nonumber \\
\bra{h_{ij} (\kk, t) h'_{ij} (\kk', t) } =
\bra{h'_{ij} (\kk, t) h_{ij} (\kk', t) } = & \, \frac{1}{4}
\bigl[{\cal M}_{ijlm} S_{\rm mix} (k, t) +
i {\cal A}_{ijlm} A_{\rm mix} (k, t) \bigr],
\end{align}
%AR: to be continued

\section{Results}

We have computed solutions for a range of values of $\sigma$, using
both given and driven initial fields.
In \Fig{pspecm_hel} we compare energy and helicity spectra of the
magnetic field and the resulting GWs for $\sigma=\pm1$ and $\pm0.01$.
The spectra are averaged over the time interval $1.1\leq t\leq 1.2$,
which is just after the maximum magnetic energy has been reached and
the GW energy begins to fluctuate around an approximately statistically
steady state.

We see that in all cases, the causal $k^4$ magnetic energy spectrum
\cite{DC03} is established.
This implies a flat spectrum for the GW energy \cite{RPMBKK19,BB19}.
For $k>k_*$, the slope of the magnetic energy spectrum is a little steeper
This is because of the finite time driving during the rather short time
interval, $1\leq t\leq 1.1$.
The $k^{-5/3}$ Kolmogorov-type spectrum would emerge if the driving
was continued over a longer time interval.
than the expected $k^{-5/3}$ spectrum.
%AB: Mention MY75 for helical Kolmogorov spectrum
%For $k>k_*$, the slope of the magnetic energy spectrum follows the 
%the helical Kolmogorov $k^{-8/3}$ spectrum \cite{MY75} for $\sigma=\pm1$,
%while for $\sigma=\pm0.01$ we have the usual Kolmogorov $k^{-5/3}$ spectrum.
%Our result for the fully helical case is closer to the $k^{-8/3}$ spectrum
%\cite{MY75} than to the $k^{-7/3}$ spectrum of Ref.~\cite{MC96}.
%AB.
%AB: taken out again

In \Fig{pspecm_hel}, the magnetic helicity spectra are multiplied by
$k/2$, which allows us to see that the spectral realizability condition
\cite{Mof78},
\EQ
|k\HM(k,t)/2|\leq\EM(k,t),
\label{RealSpec}
\EN
is nearly saturated at $k=k_*$, so the magnetic energy is here nearly
fully helical for $\sigma=\pm1$; see the upper two panels of \Fig{pspecm_hel}.
For $\sigma=\pm0.01$ (lower two panels) there is a clear separation between
$\EM(k,t)$ and $k\HM(k,t)/2$ and the sign of the magnetic helicity tends
to fluctuate noticeably.
We also see a systematic sign flip at higher $k$ both in $\HM$ and
in $H_{\rm GW}$.
Such sign flips are typical of decaying helical turbulence are a
consequence of the fact that the fractional magnetic helicity is there
already extremely small.

\begin{figure*}[t!]\begin{center}
%\includegraphics[width=.49\textwidth]{pspecm_F1152d2_sig1_t11_M2c_double}
%\includegraphics[width=.49\textwidth]{pspecm_F1152sigm1a}
%AB: new figures
\includegraphics[width=.49\textwidth]{pspecm_M1152e_exp6k4}
\includegraphics[width=.49\textwidth]{pspecm_M1152e_exp6k4_sig01}
\includegraphics[width=.49\textwidth]{pspecm_F1152sig001a}
\includegraphics[width=.49\textwidth]{pspecm_F1152sigm001a}
\end{center}\caption[]{
Magnetic and GW energy and helicity spectra.
}\label{pspecm_hel}\end{figure*}

%AB: The polarization is actually defined in terms of Ah and Sh, not Ahdot and Shdot.
%AB: I've now plotted both, and we should adapt the text below to define both.
%In \Fig{ppol_comp} we plot the GW polarization
%AB: plural
In \Fig{ppol_comp} we plot the GW polarizations
\begin{equation}
%\mathcal{P} (k) \,=\, \frac{\XiGW (k)} {\OmGW (k)}\, = \,\frac{\Ahd (k)}{\Shd (k)}\,.
%AB: can we not just omit A and S?
\mathcal{P} (k) \,=\, \frac{\XiGW (k)} {\OmGW (k)},
\end{equation}
where
\begin{equation}
%AB: combine
%\quad\mbox{where}\quad
%\XiGW (k) \,=\, \frac{k \Ahd (k)}{6 H_*^2}, \hskip 8mm \text{such that \ }\hskip 4mm
%\int_{-\infty}^{\infty} \XiGW(k) \, \dd \ln k\, =\, \XiGW\,, \hskip 8mm
%AB: modified, and used 2/6 = 1/3.
\XiGW (k) \,=\, \frac{k}{3 H_*^2}
 \int_{4\pi} {\rm Im} \left( \dot{\tilde{h}}_{+}^{\rm phys}
\dot{\tilde{h}}_{\times}^{{\rm phys}, \ast} \right)
\, k^2 \, \dd  \Omega_k\,,
\label{XiGW_app}
\end{equation}
%where
%AB: changed
%\begin{align}
%A_{\dot{h}} (k, t)\, &\, = \,\int_{\varOmega_D}\left(
%\dot{\tilde{h}}_{+}^{\rm phys}
%\dot{\tilde{h}}_{\times}^{{\rm phys}, \ast} -
%\dot{\tilde{h}}_{+}^{{\rm phys},\ast}
%\dot{\tilde{h}}_{\times}^{\rm phys} \right) \, k^{D - 1} \, \dd \varOmega_k
%\nonumber \\
%&\, =\, 2 \int_{\varOmega_D} {\rm Im} \left( \dot{\tilde{h}}_{+}^{\rm phys}
%\, =\, 2 \int_{\varOmega_D} {\rm Im} \left( \dot{\tilde{h}}_{+}^{\rm phys}
%\dot{\tilde{h}}_{\times}^{{\rm phys}, \ast} \right)
%\, k^{D - 1} \, \dd  \varOmega_k\,,
%\label{Adoth_app}
%\end{align}
%and $\varXi_{\rm GW}(t)$ is the total normalized helical energy density.
%AB: simpler
with $\int\Xi_{\rm GW}(k)\,d\ln k=\Xi_{\rm GW}$ being the total
normalized helical energy density.
%AB.
We also define the fractional magnetic helicity
${\cal P}_{\rm M}(k)=k\HM(k)/2\EM(k)$.

We see that $\mathcal{P}_{\rm GW}(k)$ reaches $\pm1$ at
$k\approx2k_*\approx1200$ when $\sigma=\pm1$.
This is more than what was found in Ref.~\cite{KK15}.
This discrepancy can probably be explained by the use of the approximations
made in the analytic calculations.
Toward smaller $k$, there are systematic fluctuations around
${\cal P}_{\rm GW}(k)=0$, especially for $k=200$.

\begin{figure*}[t!]\begin{center}
%\includegraphics[width=.70\textwidth]{ppol_comp}
\includegraphics[width=\textwidth]{ppol_comp2}
%AB: commented out for now
%\includegraphics[width=\textwidth]{pol_analytical_comp}
\end{center}\caption[]{
(a) given initial field with different fractional helicity.
(b) driven  initial field with different fractional helicity.
+ analytical formula.
%The black lines show ${\cal P}_{\rm GW}(k)$ for various values of $\sigma$.
%The red line gives the result for the runs similar to Run~ini1 of Ref.~\cite{RPMBKK19},
%for $\sigma=1$, 0.5, and 0.1.
}\label{ppol_comp}\end{figure*}

We also plot the total polarization integrated over all wave numbers, i.e.,
\begin{equation}
\tilde{\mathcal{P}}_{\rm GW}=\frac{\Xi_{\rm GW}}{\Omega_{\rm GW}},\quad
\tilde{\mathcal{P}}_{\rm M}=\frac{\HH_{\rm M}}{2\xiM\EEM}.
\end{equation}
The result is shown in \Fig{pppol1}.
We see that the numerical data can well be fitted with
the curve
\begin{equation}
%\tilde{\mathcal{P}}_{\rm GW}\approx0.6\tan \tilde{\mathcal{P}}_{\rm M}.
%AB: change 0.6 to general coefficient
\tilde{\mathcal{P}}_{\rm GW}\approx\tilde{\mathcal{P}}_0
\tan \tilde{\mathcal{P}}_{\rm M},
\end{equation}
where $\tilde{\mathcal{P}}_0\approx0.6$ for the driven case and 
$\tilde{\mathcal{P}}_0\approx0.76$ for the case with an initial
magnetic field.
This relation between $\tilde{\mathcal{P}}_{\rm GW}$ and
$\tilde{\mathcal{P}}_{\rm M}$ is unexpected and the agreement with the
actual data points is remarkable.
On the other hand, for a Beltrami field, which is monochromatic
and defined only for a single value of $k$, we have an exact
analytical solution for both $\tilde{\mathcal{P}}_{\rm GW}$ and
$\tilde{\mathcal{P}}_{\rm M}$; see Appendix~A.
The resulting dependence of $\tilde{\mathcal{P}}_{\rm GW}$ on
$\tilde{\mathcal{P}}_{\rm M}$ has a rather different shape; compare
with the dotted line for the Beltrami field.
This is surprising, but it may be related to the smaller fractional
polarization away from the peak near $k_*$; see \Fig{ppol_comp}.
This feature becomes more pronounced for smaller values of $\sigma$.
On the other hand, close to the maximally helical fields,
$\tilde{\mathcal{P}}_{\rm GW}$ shows a much more dramatic increase
with $\tilde{\mathcal{P}}_{\rm M}$.

\begin{figure}[t!]\begin{center}
\includegraphics[width=.60\textwidth]{pppol1}
%\includegraphics[width=\columnwidth]{pppol1}
\end{center}\caption[]{
GW polarization versus magnetic polarization.
}\label{pppol1}\end{figure}

\section{Possibility of detection with LISA}

To assess the observational prospects of detecting GWs for different
values of $\sigma$, we now plot $\Omega_{\rm GW}$ for $\sigma=1$,
$0.7$, $0.5$, $0.3$, and $0.1$; see \Fig{pspecm_scl5_hel_sig}.
The result for $\sigma=1$ was already shown in the top left
panel of \Fig{pspecm_hel}.
The spectra are similar, but we now also see that for smaller values
of $\sigma$, the jump in $\Omega_{\rm GW}(f)$ near the peak is
less pronounced, so for larger frequencies, i.e., to the right of the peak,
$\Omega_{\rm GW}(f)$ increases (decreases) for smaller (larger)
values of $\sigma$.
%This corresponds to the inertial range of turbulence,
%where the slope for the nonhelical or weakly helical case
%is clearly shallower (compatible with $f^{11/3}$) than
%for the fully helical case (compatible with $f^{14/3}$).
%AB: omit

For smaller frequencies, to the left of the peak, which corresponds to the
subinertial range, we have the aforementioned shallow spectrum with
$\Omega_{\rm GW}(f)\propto f$, which is approximately independent of
the value of $\sigma$.
Again, positive (negative) values are indicated by a red (blue) symbols.

\begin{figure*}[h!]\begin{center}
\includegraphics[width=.49\textwidth]{pspecm_scl5_hel_sig}
\includegraphics[width=.49\textwidth]{pspecm_scl5_hel_sig_Me}
\end{center}\caption[]{
Diagnostic diagram of $h_0^2\Omega_{\rm GW}(f)$ for
$\sigma=1$ (red), 0.7 (orange), 0.5 (green), 0.3 (blue), 0.1 (black)
for the driven case (left) and the case with initial condition (right).
The LISA sensitivity curve is shown as a green dash-dotted line and slopes
proportional to $f^{-4}$ and $f^{-5}$, as well as $f^{-8/3}$
are shown for orientation.
}\label{pspecm_scl5_hel_sig}\end{figure*}

\begin{figure*}[h!]\begin{center}
\includegraphics[width=\textwidth]{pspecm_scl5_hel_comp3}
\end{center}\caption[]{
Diagnostic diagram of $h_0^2\Omega_{\rm GW}(f)$ and $h_0^2\Xi_{\rm GW}(f)$ 
for given initial magnetic fields with (a) $\sigma=\pm1$ and (b) $\pm0.1$, and
initially driven magnetic fields with (c) $\sigma=\pm1$ and (d) $\pm0.1$.
}\label{pspecm_scl5_hel_comp3}\end{figure*}

%\begin{figure*}[h!]\begin{center}
%\includegraphics[width=\textwidth]{pspecm_scl5_hel_comp2_Me}
%\end{center}\caption[]{
%Similar to \Fig{pspecm_scl5_hel_comp2}, but for the case with
%initial magnetic energy spectra.
%}\label{pspecm_scl5_hel_comp2_Me}\end{figure*}

Finally, we plot in \Fig{pspecm_scl5_hel_comp3} the energy spectra together with
the helicity spectra $\Sigma_{\rm GW}$.
In this logarithmic representation, the fractional GW polarization
appears still rather noticeable for $\sigma=\pm0.1$.
We see that for fractional magnetic helicity, the GW power decreases and
the values of $\Xi_{\rm GW}(f)$ are now close to the LISA sensitivity limit,
$\Xi_{\rm GW}^{\rm LISA}(f)$.
On the other hand, the polarization signal is still rather clean in the
sense that the amount of statistical fluctuations in the sign of
$\Xi_{\rm GW}^{\rm LISA}(f)$ is weak.

\section{Conclusions}

Our numerical simulations have confirmed that there is a direct correspondence
between the sign of magnetic helicity and the sign of circular polarization
of the GWs produced from the resulting magnetic stress.
Our work has also revealed a new functional form of the dependence of
circular polarization of GWs on the fractional magnetic helicity.
The dependence is shallower for weakly helical magnetic fields
and steeper for nearly fully helical fields.
The precise form of this curve can be important if one wants to infer the
magnetic helicity from circular polarization measurements of GWs.
Cosmological magnetic fields may well be close to fully helical,
because the magnetic helicity is a conserved quantity while the
magnetic energy decays, so the ratio always increases until it
reaches nearly 100\%.

The subinertial range of the GW spectrum is hardly affected by the
presence or absence of magnetic helicity. However, the peak of
the GW spectrum is higher in the case with magnetic helicity.
This is another reason why magnetic helicity is beneficial for
enhancing the detection prospects of GWs, provided LISA is indeed
able to measure the GW energy at the peak frequency.
On the other hand, at larger wave numbers, corresponding to the inertial
range of the underlying turbulence, the trend is the other way around,
i.e., the spectral GW energy is lower in the presence of helicity and
the slope steeper.

\appendix

\section{Analytical model for polarization}

%AB: do you still plan to add something else here?

%AB: I don't think we need this now
%\begin{figure*}[h!]\begin{center}
%\includegraphics[width=.6\textwidth]{pol_analytical}
%\end{center}\caption[]{
%Comparison of polarization for analytical model (red) with $n_S = n_A =-13/3$
%with $h = 0.5$ and 1,
%and $n_S = -11/3$, $n_A = -14/3$ with $h = 1$, with numerical results (black).
%}\label{pol_analytical}\end{figure*}
%\begin{figure*}[h!]\begin{center}
%\includegraphics[width=\textwidth]{pol_analytical_comp_numerical}
%\end{center}\caption[]{
%}\label{pol_analytical}\end{figure*}

\section{Beltrami field}

Simple magnetic field with partial helicity
\begin{equation}
\BB(\xx, t) = B_0 \Theta(t) \left( \begin{array}{c}
0 \\ \sigma \sin k_0 x \\ \cos k_0 x
\end{array} \right),
\end{equation}
where $\sigma \in [0, 1]$ is a parameter to modify the helicity
of the magnetic field, $k_0 \in {\cal R}$ is a characteristic
wavenumber, and $\Omega_{\rm M} = B_0^2/2$ is the amplitude.
The Heaviside function $\Theta(t)$ is included to indicate that
this field is zero for $t \leq 0$.

In Fourier space this is
\begin{equation}
\tilde{\BB}(\kk, t) = \frac{1}{2} B_0 \Theta(t) \left( \begin{array}{c}
0 \\ i \sigma[ \delta^3(\kk + \kk_0) - \delta^3(\kk - \kk_0)] \\
 \delta^3(\kk + \kk_0) + \delta^3(\kk - \kk_0)
\end{array}  \right),
\end{equation}
where $\kk_0 = (k_0, 0, 0)$, and $\delta^3(\kk)$ is the Dirac's
delta function in three dimensions.
We have used the convention
\begin{equation}
\tilde{\BB}(\kk, t) = \int \BB(\xx, t) e^{-i \kk \cdot \xx} \, d \xx, \quad
\BB(\xx, t) = \frac{1}{(2 \pi)^3} \int \tilde{\BB} (\kk, t)
e^{-i \kk \cdot \xx} \, d \kk.
\end{equation}

The magnetic spectrum is
\begin{equation}
E_{\rm M} (k, t) = \frac{1}{2} \int_{\Omega_1} \tilde{\BB}(\kk, t)
\cdot \tilde{\BB}^*(\kk ,t)\,
d \Omega_1 = \frac{\Omega_{\rm M}}{2} \Theta(t) (1 + \sigma^2) \delta(k - |k_0|),
\end{equation}
where $\Omega_1 = 2$ is the one-dimensional solid angle,
and $\delta(k)$ is the one-dimensional Dirac's delta function.

The total magnetic energy density is
\begin{equation}
\EEM (t) = \int_0^{\infty} E_{\rm M} (k, t) \, d k = 
\frac{\Omega_{\rm M}}{2} \Theta(t) (1 + \sigma^2).
\end{equation}

The magnetic helicity spectrum is
\begin{equation}
H_{\rm M} (k, t)= \int_{\Omega_1} \tilde{\BB}(\kk, t) \cdot (i \kk
\times \tilde{\BB}(\kk, t))^* \, d
\Omega_1 = 2 \Omega_{\rm M} \Theta(t) |k_0| \sigma \delta(k - |k_0|).
\end{equation}

The total magnetic helicity density is
\begin{equation}
\HHM (t) = \int_0^{\infty} H_{\rm M} (k, t) \, d k = 
2 \Omega_{\rm M} \Theta(t) |k_0| \sigma.
\end{equation}

The magnetic spectral peak is
\begin{equation}
k_{\rm M} = \frac{ \EEM (t)}
{\int_0^{\infty} k^{-1} E_{\rm M} (k, t)\, d k} = |k_0|.
\end{equation}

The fractional magnetic helicity is
\begin{equation}
\varepsilon_h = \frac{\bra{\BB(\xx, t) \cdot (\nabla \times
\BB(\xx, t)}}{k_{\rm M} \bra{\BB^2 (\xx, t}} = 
\frac{\HHM(t)}{2 k_{\rm M} \EEM(t)} = 
\frac{2 \sigma}{1 + \sigma^2},
\end{equation}
which reduces to $+1, 0,$ and $-1$,
for $\sigma = +1, 0,$ and $-1$, respectively.

The stress tensor for magnetic fields is
\begin{equation}
T_{ij} (\xx, t) = - B_i B_j + \frac{1}{2} \delta_{ij} \BB^2,
\end{equation}
where
\begin{equation}
\BB^2 (\xx, t) = \Omega_{\rm M} \Theta(t)\left[(1 + \sigma^2) +
\cos 2k_0 x (1 - \sigma^2)\right],
\end{equation}
and
\begin{equation}
-B_i B_j (\xx, t) = - \Omega_{\rm M} \Theta(t) \left( \begin{array}{c c c}
0 & 0 & 0 \\
0 & \sigma^2 (1 - \cos 2 k_0 x) & \sigma \sin 2 k_0 x \\
0 & \sigma \sin 2k_0 x & 1 + \cos 2 k_0 x
\end{array} \right).
\end{equation}

Combining both we obtain the stress tensor
\begin{equation}
T_{ij} (\xx, t) = \frac{1}{2} \Omega_{\rm M} \Theta(t)\left(
\begin{array}{c c c} T_{11} (\xx) & 0 & 0 \\
0 & T_{22} (\xx) & T_{12} (\xx) \\ 0 & T_{12} (\xx)
& -T_{22} (\xx) \end{array} \right),
\label{Tij}
\end{equation}
where
\begin{gather}
%\begin{equation}
T_{11} (\xx) = (1 + \sigma^2) + (1 - \sigma^2) \cos 2k_0 x, \cr
T_{22} (\xx) = (1 - \sigma^2) + (1 + \sigma^2) \cos 2k_0 x, \cr
T_{12} (\xx) = 2 \sigma \sin 2k_0 x.
\end{gather}
%\end{equation}

Since GWs are produced by linear perturbations over the metric tensor,
and the stress tensor components are also perturbations over background
fields, constant modes ($\kk = \nullvector$) do not produce GWs.
For this reason, we can neglect the terms that are homogeneous in space.
This can be done by computing $\tilde{T}_{ij}^{\kk} (\kk, t)
= \tilde{T}_{ij} (\kk, t) (1 - \delta^3 (\kk))$.
The stress obtained is the same as in \Eq{Tij}, replacing $T_{11} (\xx)$
and $T_{22} (\xx)$ by  
%\begin{gather}
\begin{equation}
T_{11}^{\kk} (\xx) = (1 - \sigma^2) \cos 2k_0 x, \\
T_{22}^{\kk} (\xx) = (1 + \sigma^2) \cos 2k_0 x. \\
%\end{gather}
\end{equation}

In Fourier space, the resulting stress tensor is
\begin{equation}
\tilde{T}^{\kk}_{ij} (\kk, t) = \frac{1}{2} \Omega_{\rm M} \Theta(t)
\left( \begin{array}{c c c}
\tilde{T}_{11}^{\kk} (\kk) & 0 & 0 \\
0 & \tilde{T}_{22}^{\kk} (\kk) & \tilde{T}_{12} (\kk) \\
0 & \tilde{T}_{12} (\kk)       & - \tilde{T}_{22}^{\kk} (\kk)
\end{array} \right),
\end{equation}
where
%\begin{equation}
\begin{gather}
\tilde{T}_{11}^{\kk} (\kk) = \frac{1}{2} (1 - \sigma^2)
[\delta^3 (\kk + 2 \kk_0) + \delta^3 (\kk - 2 \kk_0)], \cr
\tilde{T}_{22}^{\kk} (\kk) = \frac{1}{2} (1 + \sigma^2)
[\delta^3 (\kk + 2 \kk_0) + \delta^3 (\kk - 2 \kk_0)], \cr
\tilde{T}_{12} (\kk) = i \sigma  
[\delta^3 (\kk + 2 \kk_0) - \delta^3 (\kk - 2 \kk_0)].
%\end{equation}
\end{gather}
%
%If we neglect the terms proportional to $\delta_{ij}$
%
%\begin{equation}
%T_{ij}^{\rm TT} (\xx, t) = - \Omega_{\rm M} \Theta(t) \left(\begin{array}{c c c}
%0 & 0 & 0 \\
%0 & - \frac{1}{2} (1 + \sigma^2) \cos 2 k_0 x & 2 \sigma \sin 2k_0 x \\
%0 & 2\sigma \sin 2k_0 x & \frac{1}{2} (1 + \sigma^2) \cos 2k_0 x
%\end{array} \right).
%\end{equation}

The transverse and traceless projection of the stress tensor is
\begin{equation}
\tilde{T}_{ij}^{\rm TT} (\kk, t) = (P_{il} P_{jm} - \frac{1}{2}
P_{ij} P_{lm}) \tilde{T}_{lm}^{\kk} (\kk, t),
\end{equation}
where $P_{ij} = \delta_{ij} - \hat{k}_i \hat{k}_j$ is the
projection operator being $\hat{\kk} = \kk/\sqrt{\kk^2}$.

The only wave vectors that leads to non-zero stress in
Fourier space are $\kk = 2 \kk_0 = (\pm 2 k_0, 0, 0)$.
The projection operator for these wave vectors is
$P_{ij} (\pm 2\kk_0) = \delta_{ij} - \delta_{i1} \delta_{j1}$.

%\acknowledgements
\vspace{3mm}

%r e f
\begin{thebibliography}{}

\bibitem{NV10}
A. Neronov and I. Vovk\ysci{2010}{328}{73}
%\bibitem[Neronov \& Vovk(2010)]{NV10}
%Neronov, A., \& Vovk, I.\ysci{2010}{328}{73}

\bibitem{MY75}
A. S. Monin and A. M. Yaglom. {\it Statistical Fluid Mechanics Mechanics of Turbulence}, vol. 2. MIT press, Cambridge, 1975.

\bibitem{RPBKKM18}
A.~Roper Pol, A.~Brandenburg, T.~Kahniashvili, A.~Kosowsky, and S.~Mandal\ygafd{2020}{114}{130}
{The timestep constraint in solving the gravitational wave equations sourced by hydromagnetic turbulence}

\bibitem{RPMBKK19}
A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili, and A. Kosowsky,
\yprd{2020}{102}{083512}

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

\bibitem{KBGMRP20}
T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal, and A. Roper Pol,
\sprl{2020}{2011.05556}

%\bibitem{Sub16}
%K. Subramanian\yrpp{2016}{79}{076901}
%
%\bibitem{NV10}
%A. Neronov and I. Vovk\ysci{2010}{328}{73}
%
%\bibitem{DKOT11}
%K. Dolag, M. Kachelriess, S. Ostapchenko, and R. Tom\`as\yapjl{2011}{727}{L4}
%
%\bibitem{Der86}
%D.~V.~Deryagin, D.~Y.~Grigoriev, V.~A.~Rubakov and M.~V.~Sazhin\yjour{1986}{Mod. Phys. Lett. A}{1}{593}
%
%\bibitem{DGN02}
%A. D. Dolgov, D. Grasso and A. Nicolis, Phys. Rev. D 66, 103505 (2002).
%
%\bibitem{Gogo07}
%G. Gogoberidze, T. Kahniashvili, and A. Kosowsky\yprd{2007}{76}{083002}
%
%\bibitem{Gri74}
%L. P. Grishchuk\yjour{1974}{Sov. Phys. JETP}{40}{409}
%
%\bibitem[Kamionkowski et al.(1994)]{Kamion94}
%Kamionkowski, M., Kosowsky, A., \& Turner, M. S.\yprd{1994}{49}{2837}
%
%\bibitem[Kosowsky et al.(2002)]{KMK02}
%Kosowsky, A., Mack, A., \& Kahniashvili, T.\yprd{2002}{66}{024030}
%
%\bibitem{Caprini16}
%C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller, G. Servant, and D. J. Weir\yjour{2016}{04}{001}
%
%\bibitem{STY18}
%S. Saga, H. Tashiro, \& S. Yokoyama\yprd{2018}{98}{083518}
%
\bibitem{BEO96}
A. Brandenburg, K. Enqvist, and P. Olesen\yprd{1996}{54}{1291}

\bibitem{DN13}
R. Durrer and A. Neronov\yanar{2013}{21}{62}

\bibitem{DC03}
R. Durrer and C. Caprini\yjour{2003}{J.\ Cosmol.\ Astropart.\ Phys.}{0311}{010}

\bibitem{Mof78}
H. K. Moffatt\ybook{1978}
{Magnetic Field Generation in Electrically Conducting Fluids}{Cambridge: Cambridge Univ.\ Press}

\bibitem{KK15}
L. Kisslinger and T. Kahniashvili\yprd{2015}{92}{043006}
{Polarized gravitational waves from cosmological phase transitions}

\bibitem{K65}
R. H. Kraichnan\ypf{1965}{8}
{Inertial-range spectrum of hydromagnetic turbulence}

\bibitem{EFLVW20}
J. Ellis, M. Fairbairn, M. Lewicki, V. Vaskonen, and
A. Wickens\yjcap{2020}{10}{032}
%
%\bibitem{CDK04}
%C. Caprini, R. Durrer, and T. Kahniashvili\yprd{2004}{69}{063006}
%
%\bibitem{Mof61}
%H. K. Moffatt\yjfm{1961}{11}{625}
%
%\bibitem{BS05b}
%A. Brandenburg, K. Subramanian\yana{2005}{439}{835}
%
%\bibitem{BM99}
%D. Biskamp and W.-C. M\"uller\yprl{1999}{83}{2195}
%
%\bibitem{CHB01}
%M. Christensson, M. Hindmarsh, and A. Brandenburg\ypre{2001}{64}{056405}
%
%\bibitem{WB16}
%J. M. Wagstaff and R. Banerjee\yjcap{2016}{01}{002}
%
\bibitem{BKMRPTV17}
A. Brandenburg, T. Kahniashvili, S. Mandal, A. Roper Pol, A. G. Tevzadze, and T. Vachaspati\yprd{2017}{96}{123528}

\bibitem{BB19}
A. Brandenburg and S. Boldyrev\yapj{2020}{892}{80}
{The turbulent stress spectrum in the inertial and subinertial ranges}

%\bibitem{CCFTV15}
%W. Chen, B. D. Chowdhury, F. Ferrer, H. Tashiro, T. Vachaspati\ymn{2015}{450}{3371}
%
\bibitem{KGR05}
T. Kahniashvili, G. Gogoberidze, and B. Ratra\yprl{2005}{95}{151301}

\bibitem{MC96}
S. S. Moiseev and O. G. Chetiani\yjetp{1996}{83}{192}

%\bibitem{KGR08}
%T. Kahniashvili, G. Gogoberidze, and B. Ratra\yprl{2008}{100}{231301}
%
%\bibitem{KCGMR08}
%T. Kahniashvili, L. Campanelli, G. Gogoberidze, Y. Maravin, and B. Ratra\yprd{2008}{78}{123006}
%
%\bibitem{TWW92}
%M. S. Turner, E. J. Weinberg, and L. M. Widrow\yprd{1992}{46}{2384}
%
%\bibitem{EarlyTimes}
%See Supplemental Material, Sec.~2, for a demonstration of the $k^2$
%subinertial range in $\EGW(k,t)$ at early times.
%
%\bibitem{Comp}
%See Supplemental Material, Sec.~3, for a comparison with the spectrum
%obtained in Ref.~\cite{Gogo07}.
%
%\bibitem{AC14}
%M. R. Adams and N. J. Cornish\yprd{2014}{89}{022001}
%
%\bibitem{MB06}
%A. J. Mee and A. Brandenburg\ymn{2006}{370}{415}
%
%\bibitem{HHRW15}
%M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir\yprd{2015}{92}{123009}
%
%\bibitem{HHRW17}
%M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir\yprd{2017}{96}{103520}
%
%\bibitem{NSS18}
%P. Niksa, M. Schlederer, and G. Sigl\ycqg{2018}{35}{144001}

\end{thebibliography}

\end{document}
