% $Id: paper.tex,v 1.59 2021/01/15 06:19:17 roper Exp $
\documentclass[a4paper,11pt]{article}
\usepackage{jcappub}

\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{\tilff}{\BoldVec{\tilde 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\OmM{{\Omega}_{\rm M}}
\def\OmMmax{{\Omega}_{\rm M}^{\rm max}}
\def\OmGWsat{{\Omega}_{\rm GW}^{\rm sat}}
\def\PPM{{\cal P}_{\rm M}}
\def\sigM{\sigma_{\rm M}}
\def\PPGW{{\cal P}_{\rm GW}}
\def\PPh{{\cal P}_h}
\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\FM{F_{\rm M}}
\def\GM{G_{\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}}}
\def\onetwelveth{{\textstyle{1\over12}}}
%
\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}[4]{, {\em #4}, Science {\bf #2}, #3 (#1).}
\newcommand{\ysph}[3]{, Solar Phys. {\bf #2}, #3 (#1).}
\newcommand{\yjetp}[4]{, {\em #4}, Sov. Phys. JETP {\bf #2}, #3 (#1).}
\newcommand{\yspd}[3]{, Sov. Phys. Dokl. {\bf #2}, #3 (#1).}
\newcommand{\ysov}[3]{, Sov. Astron. {\bf #2}, #3 (#1).}
\newcommand{\ysovl}[3]{, Sov. Astron. Lett. {\bf #2}, #3 (#1).}
\newcommand{\ymn}[3]{, Mon.\ Not.\ R.\ Astron.\ Soc.\ {\bf #2}, #3 (#1).}
\newcommand{\ymhd}[3]{, Magnetohydrohydrodyn. {\bf #2}, #3 (#1).}
\newcommand{\yqjras}[3]{, Quart. J. Roy. Astron. Soc. {\bf #2}, #3 (#1).}
\newcommand{\ynat}[3]{, Nature {\bf #2}, #3 (#1).}
\newcommand{\yjfm}[3]{, J. Fluid Mech. {\bf #2}, #3 (#1).}
\newcommand{\pjfm}[1]{, J. Fluid Mech., in press (#1).}
\newcommand{\sjfm}[1]{, J. Fluid Mech., submitted (#1).}
\newcommand{\ypr}[3]{, Phys.\ Rev.\ {\bf #2}, #3 (#1).}
\newcommand{\yprd}[4]{, {\em #4}, {\em Phys.\ Rev.} {\bf D #2} (#1) #3.}
\newcommand{\pprr}[3]{, {\em #3}, {\em Phys.\ Rev.\ Res.} (#1) arXiv#2.}
\newcommand{\ypre}[3]{, Phys.\ Rev.\ E {\bf #2}, #3 (#1).}
\newcommand{\yprl}[4]{, {\em #4}, 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}[4]{, {\em #4}, {\em Astrophys.\ J.} {\bf #2}, (#1) #3.}
\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}[4]{, {\em #4}, 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]{, {\em #3}, Phys. Fluids {\bf #2} (#1).}
\newcommand{\yphy}[3]{, Physica {\bf #2}, #3 (#1).}
\newcommand{\ygafd}[4]{, {\em #4}, 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}[5]{, {\em #5}, #2 {\bf #3}, #4 (#1).}
\newcommand{\pjour}[3]{, #2, in press (#1).}
\newcommand{\sjour}[3]{, #2, submitted (#1).}
\newcommand{\yprep}[2]{, #2, preprint (#1).}
\newcommand{\pproc}[3]{, (ed. #3), #2 (#1) (to appear).}
\newcommand{\yproc}[4]{, (ed. #4), pp. #2. #3 (#1).}
\newcommand{\ybook}[3]{, {\em #2}. #3 (#1).}

%\begin{document}

\title{Polarization of gravitational waves from helical MHD turbulent sources}


\date{\today,~ $ $Revision: 1.59 $ $}
%\preprint{NORDITA-2018-xx}

\author[a,b,1]{Alberto~Roper~Pol}
\emailAdd{roperpol@apc.in2p3.fr}
\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}
\emailAdd{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}
\emailAdd{brandenb@nordita.org} 
\affiliation[d]{Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfv\'ens v\"ag 12, 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}
\emailAdd{tinatin@andrew.cmu.edu} 
\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}
%\emailAdd{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}

\emailAdd{\today}

\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 turbulent source term, which can be magnetic or kinetic,
is either given initially, or it is produced by 
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 
helicity of the source (magnetic or kinetic).
In both cases, we find a frequency tail below the spectral peak
shallower than the $f^3$ scaling expected at super-horizon scales,
in agreement with similar recent numerical simulations.
This type of spectrum facilitates its observational detection with
the planned Laser Interferometer Space Antenna (LISA).
We show that the magnetic fields driven for a short time produce
GWs more efficiently than the case with initial magnetic fields, leading
to larger values of the spectral peak.
The helicity does not seem to have a huge impact on the spectral peak,
being slightly larger for smaller fractional helicity.
For frequencies above and below the peak, the two types of simulations
show that the spectral energy is lower in the helical case than in the
%nonhelical one, being the differences more noticeable when the magnetic
%AB: moved position of "being"
nonhelical one, the differences being more noticeable when the magnetic
field is 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 at all wave numbers, and it approaches zero if the
initial field is driven by a fully helical monochromatic forcing function.
We also show that in all cases with a driven turbulence spectrum,
there is a sharp drop of GW energy
above the spectral peak
that is in agreement with the previously obtained results.
The drop becomes sharper for larger fractional helicity.
We predict that the spectral peak of the GW signal can be detected
by LISA if the magnetic energy density is $\OmM=0.1$ and
$\OmM=0.05$ for the initial and driven field cases, respectively.
Additionally, the polarization signal around the peak could be
detected by cross-correlating detections of the LISA-Taiji network
%in the driven case with initial $\OmM=0.1$ and magnetic helicity of
%AB: inserted fractional
in the driven case with initial $\OmM=0.1$ and fractional magnetic helicity of
0.5 to 1.
%Second-generation space-base GW detectors, as BBO and DECIGO, would
%AB: added "such"
Second-generation space-base GW detectors, such as BBO and DECIGO, would
%allow the detection of a larger range of the GW spectrum and smaller
%AB: added "for"
allow for the detection of a larger range of the GW spectrum and smaller
amplitudes of the magnetic field.
}
%\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 a stochastic
background of 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 (EW), 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.
%Future space-based GW detectors as the Laser Interferometer Space Antenna
%AB: added such
Future space-based GW detectors such as the Laser Interferometer Space Antenna
%LISA), planned to launch in 2034, as well as TianQin and Taiji, will be sensitive to
%AB: launch -> to be launched 
LISA), planned to be launched in 2034, as well as TianQin and Taiji, will be sensitive to
GWs in frequencies ranging from 10$\mu$Hz to a few Hz, with a peak
sensitivity around 1 mHz (LISA, Taiji) and a few mHz (TianQin).
%In addition, the next generation of space-based GW detectors are planned
%AB: generation is singular, so are -> is
In addition, the next generation of space-based GW detectors is planned
to improve the sensitivity to GW signals, and covering the range from
mHz to 10 Hz, e.g., the Deci-hertz Interferometer Gravitational Wave
%AB: shouldn't Hertz be uppercase?
%AR: I've seen it written like this in their papers (e.g.,https://arxiv.org/pdf/2006.13545.pdf), alhough personally I agree it looks better with Hertz
%mHz to 10 Hz, e.g., the Deci-Hertz Interferometer Gravitational Wave
Background (DECIGO), and the Big Bang Observer (BBO).
%TK reference for websites?
%AR: yes, more details and references to be added
%AB: oops, I just deleted your comment above, but they should not be as comments
%AB: but as references.
In this range of frequencies we expect to observe GW signals generated
at the EW phase transition.
%In the lower regime of frequencies, measuring the time of arrival of a
%AB: this sentence didn't read well; changed of -> by
In the lower regime of frequencies, measuring the time of arrival by a
%pulsar timing array (PTA) can allow to detect GW signals in the range from
%AB: changed can allow -> allows one 
pulsar timing array (PTA), allows one to detect GW signals in the range from
$10^{-9}$--$10^{-7}$ Hz, allowing to detect signals generated during
the QCD phase transition.
%Currently, the International Pulsar Timing Array (IPTA) is a collaboration that
%unifies the effort of the European Pulsar Timing Array (EPTA),
%%TK reference for websites? 
%the North American Nanohertz Observatory for Gravitational Waves (NANOGrav),
%and the Parkes Pulsar Timing Array (PPTA).
%AR: maybe not really needed, since we focus on EWPT for polarization detection
%AB: ok
Even smaller frequencies can be probed by the indirect detection of
B-modes in the CMB polarization, which can be produced by inflationary
produced GW signals.
%AR: to be added more details and references
%TK agree
%AB: We don't have currently any refereces in the first paragraph

There is various kind of evidence for primordial magnetic fields.
They could have survived until the present
time \cite{Brandenburg:2017neh} and are thus a leading candidate for explaining the
lower limits on the strength of extragalactic magnetic fields inferred
from the non-observation of cascade from TeV photons of blazars
interacting with the cosmic background light;
see Ref.~\cite{Neronov:1900zz} for a pioneering work and Ref.~\cite{Vachaspati:2020blt}
for a review and references therein.
%Such fields are expected to be turbulent and partially helical.
%TK some words why are they turbulent. I think helicity should be mentioned later, after sources are identified
%AR: modified, and text added
These fields dynamically interact with the primordial plasma, leading to
magnetohydrodynamic (MHD) turbulence.
Parity violating (or charge-parity violating) processes in the early Universe
can lead to the generation of helical magnetic fields via sphaleron decay
(see Refs.~\cite{Vachaspati:1991nm, Cornwall:1997ms,Vachaspati:2001nb}
for pioneering work, and Ref.~\cite{Vachaspati:2020blt} for a recent review),
and due to inhomogeneities in the Higgs field in low-scale EW hybrid inflation \cite{GarciaBellido:2002aj}.
Helical magnetic fields can also be produced
%during inflation; see Ref.~\cite{Sharma:2019jtb} for recent work.
%AB: added OF21 here too
during inflation; see Refs.~\cite{Sharma:2019jtb,OF21} for recent work.
Such fields
produce circularly polarized GWs \cite{Kahniashvili:2005qi, Kisslinger:2015hua, Kahniashvili:2020jgm}.
In particular, the detection of circularly polarized GWs will shed light
on phenomena of fundamental symmetry breaking in the early universe,
such as parity violation, and potentially can serve as an explanation of the
lepto- and baryogenesis problem.
Additionally, the evolution of the coherence length scale of the magnetic fields
depends on the magnetic helicity, since inverse cascading is enhanced by helicity,
leading to larger scales at the present time for helical fields.
This could explain the origin of the present-time fields observed by Fermi.

The degree of polarization of GWs, which is obviously
determined by the helicity fraction of the source,
has been a matter of uncertainty owing to the approximations made in the analytical
calculations available to date. 
In particular, previous works (see
Refs.~\cite{Kahniashvili:2005qi,Kisslinger:2015hua}) suggested that the circular polarization of
GWs would be at most about 80\% for a maximally helical magnetic field, by
assuming Kolmogorov-type turbulence, with spectral index of $-5/3$ for
the magnetic energy density and $-8/3$ for the helicity%
\footnote{The spectral indices in Refs.~\cite{Kahniashvili:2005qi,Kisslinger:2015hua}
refer to spectra defined with a $1/k^2$ factor with respect to those used
in the present work, where a scale-invariant spectrum is $\propto k^{-1}$.
Hence, the spectral indices $-5/3$ and $-8/3$ correspond to $-11/3$ and
$-14/3$ in their works.},
following phenomenological modeling of Ref.~\cite{Kraichnan:1965zz}.
The GW spectral peak is at twice the magnetic 
field energy density spectral peak\footnote{The GW energy density is
sourced by the stress tensor, computed from the convolution of
the magnetic field in Fourier space.
This produces the GW peak to be twice the magnetic peak.}.
For low values of the fractional magnetic helicity, the maximum
degree of polarization also diminishes.
Reference~\cite{Kahniashvili:2005qi} showed that the maximum circular polarization depends
on the relation between the magnetic energy and magnetic helicity spectra,
and it can actually reach nearly 100\% in the case
when the spectral indices are equal, which we call a Moiseev-Chkhetiani 
type spectrum, following Ref.~\cite{MC96} or {\it helical turbulence}
spectrum.
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. 
More recently, the two types of spectra mentioned above have been studied
in Ref.~\cite{Ellis:2020uid} in the context of detectability with LISA,
by using the dipole response induced by the proper motion of the solar system,
as proposed in Refs.~\cite{Seto:2006hf,Seto:2006dz}, and recently applied to
LISA in Ref.~\cite{Domcke:2019zls}.
The detection can be improved by cross-correlating two space-based GW detectors
as, for example, LISA and Taiji \cite{Orlando:2020oko}, or LISA and TianQin.
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 whether the magnetic
field is present at the initial time, or forced for a short duration.

Recently, in Refs.~\cite{Pol:2018pao,Pol:2019yex},
the authors have described the
implementation of a GW solver into the {\sc Pencil Code} \cite{PC} and
have presented direct numerical simulations for modeling development and dynamics of
primordial hydrodynamic and
hydromagnetic turbulence from phase transitions,
and subsequent generation of a stochastic GW
background, also computed numerically.
Their simulations included fully helical sources, but
the estimation of the GW polarization degree spectra 
as well as the polarization detection prospects
were not the focus of their studies.
We present here two types of simulations, similar to the two types of
hydromagnetic simulations presented in Ref.~\cite{Pol:2019yex}: one where
a primordial magnetic field was assumed given as initial condition,
and one where a magnetic field was generated by an electromotive force
$\FFF(\xx,t)$ that depends on time $t$ and 
position $\xx$.
In particular, regarding the first type,
we study the cases with an initial stochastic
magnetic field with different values of the fractional helicity,
up to the fully helical case.
In Ref.~\cite{Kahniashvili:2020jgm}, the degree of circular polarization for
kinetically and magnetically forced turbulence was presented,
which is similar to our second type of simulations.
We complement this study in the present work by studying the variation
of the polarization degree 
in the different scenarios of the magnetic field generation.
Furthermore, we study cases where turbulence is driven for times
significantly shorter than what was considered in
Ref.~\cite{Kahniashvili:2020jgm}.
The driving is then applied during a short time interval
(around a 10\% of the Hubble time),
and then switched off 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{Pol:2018pao}.
We explore the detectability of the generated GW signal and its
polarization with planned space-based GW detectors.

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}

\subsection{Basic equations}

We perform direct numerical simulations of the MHD turbulence starting
at the time of generation (which can be appropriately scaled to, e.g.,
the EW phase transition).
At every time step of the MHD simulation, we compute the
contributions from velocity and magnetic fields to the stress tensor
$T_{ij}$.
Then, we solve the GW equation to compute the strains $h_{ij}$,
sourced by the traceless and transverse projection of the stress tensor.
The details of the numerical setup and application to the EW scale
is described in Refs.~\cite{Pol:2018pao,Pol:2019yex}.

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\footnote{We use the Fourier convention, such that
$\tilde h (\kk) = \int h(\xx) e^{-i \kk \cdot \xx} \dd^3 \xx$, and
the inverse Fourier transform is $h (\xx) = (2 \pi)^{-3}
\int \tilde h (\kk) e^{i \kk \cdot \xx} \dd^3 \kk$.}.
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,
\label{basis_pol}
\end{equation}
where $\ee^1$ and $\ee^2$ form a basis with the units vectors $\kk/|\kk|$
\cite{Pol:2018pao}.
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$\footnote{We take the scale factor $a$ to be unity at the time
of generation, which allows to simply write $a = t$, and $t_*=H_*^{-1}$,
where $H_*$ is the Hubble rate at the time of generation.
For different convention, one should substitute $t$ by $a$ in
the denominator of the sourcing term of \Eq{GW}.},
and comoving wave vector, normalized by the corresponding scale
$1/(c t_\ast)$
\cite{Pol:2018pao,Pol:2019yex}
\begin{equation}
\left( \partial^2_t + \kk^2\right)
\tilde{h}_\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 scaled strains are tensor mode perturbations over
the Friedmann-Lema\^itre-Robertson-Walker metric tensor.
Hence, the metric tensor is such that the line element is
$\dd s^2 = a^2(-\dd t^2 + \delta_{ij}
+ h_{ij}/a)$ \cite{Pol:2018pao,Pol:2019yex}. 
During the radiation-dominated epoch, the equation of state is
$p = \onethird \rho$, which leads to a linear evolution of
the scale factor $a$ with $t$, and allows to get rid of the
damping term, which should be included otherwise in \Eq{GW}.

The stress is composed of magnetic and kinetic contributions
and computed in physical space as
\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.
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.

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{Pol:2018pao,Pol:2019yex,Brandenburg:1996fc,Brandenburg:2017neh,Durrer:2013pga}
\begin{eqnarray}
{\partial\ln\rho\over\partial t}
&=&-\frac{4}{3}\left(\nab\cdot\uu+\uu\cdot\nab\ln\rho\right)
+{1\over\rho}\left[\uu\cdot(\JJ\times\BB)+\eta \JJ^2\right],
\nonumber \\
{\partial\uu\over\partial t}&=&-\uu\cdot\nab\uu
+{\uu\over3}\left(\nab\cdot\uu+\uu\cdot\nab\ln\rho\right)
+{2\over\rho}\nab\cdot\left(\rho\nu\SSSS\right) \nonumber \\
&&-{1\over4}\nab\ln\rho -{\uu\over\rho}\left[\uu\cdot(\JJ\times\BB)+\eta\JJ^2\right] +{3\over4\rho}\JJ\times\BB,
\nonumber \\
{\partial\BB\over\partial t}&=&\nab\times(\uu\times\BB-\eta\JJ+\FFF),\quad
\JJ=\nab\times\BB,
\nonumber
\end{eqnarray}
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.
The electromotive force $\FFF$ is used to model the generation of magnetic
fields.

\subsection{Magnetic fields present at the initial time}

In the first type of runs, we consider the magnetic field to be present
at the initial time of the simulation, so
we set the electromotive force term $\FFF$ to be zero at all times.
We generate a random three-dimensional vector field in Fourier space,
\begin{equation}
B_i({\kk})=B_0\left(P_{ij}(\kk)-i \sigma\epsilon_{ijl}
k_l/k\right) g_j({\kk})\, g_0(k),
\label{Bikk}
\end{equation}
where $B_0$ is the magnetic field amplitude,
$g_j(\kk)$ is the Fourier transform of a $\delta$-correlated
vector field in three dimensions with Gaussian fluctuations,
i.e., $g_i(\xx) g_j(\xx') = \delta_{ij} \delta^3(\xx - \xx')$,
$\sigma$ is a parameter that allows to control the magnetic helicity,
$P_{ij} = \delta_{ij} - k_i k_j/k^2$ is the projection operator,
and $g_0(k)$ determines the spectral shape \cite{Brandenburg:2017neh}.
\begin{equation}
g_0(k)={k_*^{-3/2} (k/k_*)^{\alpha/2-1}
\over[1+(k/k_*)^{2(\alpha-\beta)}]^{1/4}},
\label{Sfunction}
\end{equation}
where $\kf$ is the spectral peak, which is identified with
the initial wave number of the energy-carrying eddies.
The magnetic energy spectrum is computed from $\half \bra{\BB^2}$%
\footnote{Angle brackets denote ensemble average over stochastic realizations,
which can be approximated as the average of the random field over the
physical domain.},
such that the initial magnetic energy density is $\OmMmax = \half B_0^2$.
Due to the normalization used, this corresponds to a fraction of the
radiation energy density at the time of magnetic field generation, and since
this case corresponds to decaying turbulence, $\OmMmax$ is the maximum
value of the magnetic energy density.
Hence, it is proportional to $k^2 g_0^2 (k)$ in Fourier space, which
corresponds to a spectral index $\alpha=4$ in the low wave number limit
(subinertial range) for a Batchelor spectrum,
and Kolmogorov-type spectral slope $\beta=-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_*$.
The exponents $2$ and $1/4$ in the denominator determine the transition
from one slope to the other around the spectral peak.
%AR: moved here from below, modified, and added footnote
%AB: What we later define as P_M is different.
%AR: agree, this is clarified in the footnote, do you think it is clear enough?
%AB: ok
The helicity's degree of the initial magnetic field, $\PPM =
\kf \bra{\AAA \cdot \BB}/ \bra{\BB^2}$\footnote{The general definition
of $\PPM$ uses a characteristic wave number $k$ computed from the integration
of the helical spectrum over wave numbers, which is, in general, different than the
spectral peak $\kf$.}, is given by
$2\sigma/(1+\sigma^2)$, being $\AAA$ the magnetic vector potential,
such that $\BB = \nabla \times \AAA$.

\subsection{Magnetic fields forced at the initial time}

In the second type of simulations, the electromotive force, $\FFF$, is used
to model the magnetic field generation with fractional magnetic helicity,
which is non-zero for a short amount of time, and its value is given by
\begin{equation}
\FFF(\xx,t)={\rm Re}\{{\cal A}{\tilff}(\kk(t))\exp[i\kk\cdot\xx+i\varphi]\},
\label{FFF_model}
\end{equation}
where the wave vector $\kk(t)$ and the phase $\varphi(t)$ change randomly
from one time step to the next.
This forcing function is therefore white noise in time and consists of plane
%waves with average wavenumber $k_*$ such that $|\kk|$ lies in an interval
%AB: here and elsewhere, isn't JCAP writing wave number as 2 words (unlike ApJ)?
waves with average wave number $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.
%AR: what is the value of N? This should be the amplitude of the forcing, right? If so, maybe using A (for amplitude) or similar might be more appropriate?
%AB: yes, I guess it is fine like so.
Here, ${\cal A}$ is the amplitude of the forcing term.
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},\;
\label{forcing}
\end{equation}
where ${\tilff}^{\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 parameter $\sigma \in [-1,1]$ is related to the helicity fraction of the 
forcing term, with $\sigma=0$ and $\sigma=\pm 1$ corresponding to nonhelical
and maximally helical cases.
The forcing is only enabled during an arbitrarily short time interval
$1\leq t\leq \tmax$, to reproduce the more realistic scenario in which the magnetic
field does not appear abruptly, but it is built up to its
maximum value $\OmMmax$ at $\tmax$, and then it decays.
We chose $\tmax = 1.1$ in the present work,
which corresponds to a 10\% of the Hubble time.
In Ref.~\cite{Kahniashvili:2020jgm}, the authors consider forcing
up to $\tmax = 3$, so the forcing is active for 2 Hubble times, although
the amplitude is considered to decrease linearly.

\subsection{Characterization of stochastic magnetic and strain fields}

The magnetic fields considered, the resulting velocity fields,
and the tensor mode perturbations $h_{ij}$, are all stochastic fields, which
can be described using the several correlation functions.
The equal time autocorrelation function of the magnetic field,
assuming statistical homogeneity and isotropy, and a
Gaussian-distribution in space\footnote{For random fields with Gaussian 
distribution, the autocorrelation two-point function is sufficient to describe
its statistical properties.}, is
\begin{equation}
\bra{\tilde B_i^*(\kk, t) \tilde 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 $\EM (k, t)$, $\HM (k, t)$ are the magnetic and helical
spectra, respectively.
We work here with magnetic energy spectra per linear wave number
interval.

The magnetic field is either given at the initial time (see \Eq{Bikk}), or
driven using the forcing term, $\FFF$, described in \Eqs{FFF_model}{forcing}
for a short time, $1 \leq t \leq \tmax$, being $\tmax=1.1$.
In both cases, we have introduced a parameter $\sigma$ that allows
to control the fractional helicity of the initial magnetic
field or the initial forcing term \cite{Brandenburg:2017neh},
\begin{equation}
\PPM (t) = \frac{\displaystyle\int_0^\infty k \HM (k, t) \, \dd k}
{2\displaystyle\int_0^\infty \EM(k, t) \, \dd k} = \frac{2\sigM (t)}
{1 + \sigM^2 (t)}.
\label{PPM}
\end{equation}
For later times, the values of $\sigM(t)$, and $\PPM(t)$ are given by the
dynamical evolution of the MHD fields.
The magnetic polarization spectrum is directly computed from the helical and
magnetic spectra,
$\PPM(k, t) = \half k \HM(k, t)/\EM(k, t)$.
When we compare magnetic and helicity spectra, the latter are multiplied by
$k/2$, due to the realizability condition
\cite{Mof78},
\EQ
|k\HM(k,t)/2|\leq\EM(k,t),
\label{RealSpec}
\EN
such that the magnetic polarization $\PPM$, defined in \Eq{PPM},
takes values from -1 to 1.
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,
\label{OmGW}
\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$ from now on.
The equal time correlation function for general tensor fields
$\Pi_{ij}^a$, $\Pi_{ij}^b$ (assuming isotropic and homogeneous random fields) is
expressed as \cite{Caprini:2003vc}
\begin{align}
\bra{\tilde \Pi^a_{ij} (\kk, t) \tilde \Pi^b_{ij} (\kk', t) } = & \, \quarter
\biggl[{\cal M}_{ijlm} \frac{S (k, t)}{4\pi k^2} +
i {\cal A}_{ijlm} \frac{A (k, t)}{4\pi k^2} \biggr] (2\pi)^6 \delta^3(\kk - \kk'),
\label{etc}
\end{align}
where
\begin{align}
{\cal M}_{ijlm} = & \, P_{il} P_{jm} + P_{im} P_{jl} -
P_{ij} P_{lm}, \nonumber \\
{\cal A}_{ijlm} = & \, \half \hat{\kk}_q (P_{jm} \varepsilon_{ilq} +
P_{il} \varepsilon_{jmq} + P_{im} \varepsilon_{jlq} +
P_{jl} \varepsilon_{imq}).
\end{align}
To compute the spectral functions of the GW energy density, we
use \Eq{etc} applied to the strains $h_{ij}$ and its time derivatives
$h'_{ij}$: when $\Pi^a_{ij} = \Pi^b_{ij} = h_{ij}$, we define
$S_h, A_h$; when $\Pi^a_{ij} = \Pi^b_{ij} = h'_{ij}$, we define
$S_{h'}, A_{h'}$, and when $\Pi^a_{ij} = h_{ij}, \Pi^b_{ij} =
h'_{ij}$ (or viceversa, note that \Eq{etc} is symmetric in $a, b$),
we define $S_{\rm mix}, A_{\rm mix}$.
The GW energy density, $\OmGW(k, t)$, and the helical/chiral spectra,
$\XiGW (k, t)$, are
\begin{align}
a^4 \OmGW(k, t) = & \, \frac{1}{12} \biggl( \frac{H_*}{H_0} \biggr)^2
k \biggl[S_{h'} (k, t) + \frac{1}{t^2} S_h (k, t) - \frac{2}{t}
S_{\rm mix} (k, t) \biggr] \nonumber \\
= & \, \bigl(H_*/H_0 \bigr)^2 \, k \EGW (k, t),
\label{OmGW_k} \\
a^4 \XiGW(k, t) = & \, \frac{1}{12} \biggl( \frac{H_*}{H_0} \biggr)^2
k \biggl[A_{h'} (k, t) + \frac{1}{t^2} A_h (k, t) - \frac{2}{t}
A_{\rm mix} (k, t) \biggr] \nonumber \\
= & \, \bigl(H_*/H_0 \bigr)^2 \, k \HGW (k, t),
\label{XiGW_k} 
\end{align}
such that $\OmGW(t) = \int \OmGW(k, t) \, \dd \ln k = \int \EGW(k, t) \, \dd k$.
Using the $+,\times$ polarization basis defined in \Eq{basis_pol},
the functions $S$ and $A$ for a generic tensor field $\Pi_{ij}$ can be
expressed as \cite{Caprini:2003vc, Pol:2018pao}
\begin{align}
2\bra{\tilde \Pi_+^a (\kk) \tilde \Pi_+^{b,*}(\kk') +
\tilde \Pi^a_\times(\kk) \tilde \Pi_\times^{b,*} (\kk')} = &
\, (2\pi)^6 \delta^3(\kk - \kk') \frac{S(k)}{4 \pi k^2}, \label{SPiij}\\
2 \bra{\tilde \Pi^a_\times(\kk) \tilde \Pi_+^{b,*} (\kk') -
\tilde \Pi^a_+(\kk) \tilde \Pi_\times^{b,*} (\kk')} = &
\, i (2 \pi)^6 \delta^3(\kk - \kk') \frac{A(k)}{4\pi k^2},
\label{APiij}
\end{align}
which allows to compute the spectral functions in \Eqs{OmGW_k}{XiGW_k}
from $\tilde h_\lambda$, $\tilde h'_\lambda$, solutions to the GW \Eq{GW},
via shell-integration\footnote{The spectral functions $S(k)$ and $A(k)$
defined as in \Eqs{SPiij}{APiij} correspond to shell-integrated functions
of the tensor fields.
For example,
$$S(k,t) = \frac{2}{(2\pi)^6}  \int 4 \pi k^2 (\tilde \Pi_+^2(\kk, t) +
\tilde \Pi_\times^2 (\kk, t)) \dd \Omega_{k},$$
where $\Omega_{k}$ is the solid angle of the shell of size $k$,
such that $\int S(k, t) \dd k = \bra{\Pi_{ij} (\xx,t) \Pi_{ij} (\xx, t)}$.}.
The degree of circular polarization of the GW background is
$\PPh (k, t) = A_h(k, t)/S_h(k, t)$ \cite{Kahniashvili:2005qi}.
We also define the polarization using the GW energy density spectral
functions: $\PPGW (k, t) = \HGW(k, t)/\EGW (k, t)$\footnote{%
In previous works, these two definitions are used interchangeably, assuming that
$\EGW = k^2 S_h$, and $\HGW = k^2 A_h$,
which holds in the absence of sources.
However, turbulent sources can modify the dispersion relation of the strains,
and $\PPh \neq \PPGW$, in general.}.
The total GW polarization can be expressed as the ratio of the integrated
spectra over wave numbers, $\PPGW (t) = \XiGW(t)/\OmGW(t)$.

\section{Numerical results}

We have computed solutions for a range of values of $\sigma$, using
both given and driven initial fields (see
Table~\ref{runs} for a summary of the different runs).
In general, the different GW modes grow up to $\delta t = t - \tmax \sim {\cal O}
(k^{-1})$; see Fig.~3 of Ref.~\cite{Pol:2019yex}, when they start to oscillate,
taking $\tmax = 1$ in the runs with an initial magnetic field.
We have set up all the runs to have a spectral peak $\kf \approx 600$
(see \Tab{runs}), and the smallest wave number of the simulations is
$k_0\approx 100$.
Hence, the GW spectrum stops growing and enters the oscillatory stage
at $\delta t \sim 10^{-2}$, and after that time we average the spectra
over oscillations in time to obtain the saturated GW spectra,
and the value of $\OmGWsat$, which corresponds to the saturated value of
the total GW energy density.

\begin{table}[b!]
    \centering
    {\footnotesize
    \begin{tabular}{l|l|l|c|c|l|l|c}
         Type & \multicolumn{1}{|c|}{$\sigma$} & \multicolumn{1}{|c|}{$\sigM$} &
            $\OmMmax$ & $\OmGWsat$ & \multicolumn{1}{|c|}{$\PPM$} &
            \multicolumn{1}{|c|}{$\PPGW$} & $\kf/(2\pi)$\\
         \hline
         ini  & \ 0.1  & \ 0.0977    & $3.959 \times 10^{-3}$ & $2.530 \times 10^{-11}$
            & \ 0.1936 & \ 0.1837
            & 100 \\  % M1152e_exp6k4_sig01
         ini  & \ 0.3  & \ 0.2995    & $4.264 \times 10^{-3}$ & $3.337 \times 10^{-11}$
            & \ 0.5497   & \ 0.5440
            & 100 \\ % M1152e_exp6k4_sig03
         ini  & \ 0.5  & \ 0.5017    & $4.884 \times 10^{-3}$ & $4.414 \times 10^{-11}$
            & \ 0.8016   & \ 0.7828
            & 100 \\ % M1152e_exp6k4_sig05
         ini  & \ 0.7  & \ 0.7076    & $5.817 \times 10^{-3}$ & $6.539 \times 10^{-11}$
            & \ 0.9430   & \ 0.9080 
            & 100 \\ % M1152e_exp6k4_sig07
         ini  & \ 1   & \ 0.9007    & $7.805 \times 10^{-3}$ & $1.180 \times 10^{-10}$
            & \ 1.0045   & \ 0.9504 
            & 100 \\ % M1152e_exp6k4
         forc & -0.01 & -0.002921   & $1.440 \times 10^{-2}$ & $2.654 \times 10^{-9}$
            & \ 0.005843 & \ 0.02373
            & 100 \\ % F1152sigm001a
         forc & \ 0.01  & \ 0.01023 & $1.441 \times 10^{-2}$ & $2.645 \times 10^{-9}$
            & \ 0.02045  & \ 0.06038
            & 100 \\ % F1152sig001a
         forc & \ 0.3   & \ 0.3084  & $1.440 \times 10^{-2}$ & $3.114 \times 10^{-9}$
            & \ 0.5632   & \ 0.5915
            & 100 \\ % F1152sig03a
         forc & \ 0.5   & \ 0.5215  & $1.932 \times 10^{-2}$ & $3.542 \times 10^{-9}$
            & \ 0.8200   & \ 0.7849
            & 100 \\ % F1152sig05a
         forc & \ 0.7   & \ 0.7521  & $2.106 \times 10^{-2}$ & $3.954 \times 10^{-9}$
            & \ 0.9608   & \ 0.8874
            & 100 \\ % F1152sig07c  0.02125 4.147314037823668e-09
         forc & -1    & -0.9523     & $2.125 \times 10^{-2}$ & $4.147 \times 10^{-9}$
            & -0.9988  & -0.9111
            & 100 % F1152sigm1a
    \end{tabular}
    }
    \caption{Summary of runs}
    \label{runs}
\end{table}

\subsection{Runs with initial magnetic field}

In \Fig{pspecm_hel_initial}, we compare the GW energy and helicity spectra of the
magnetic field and the resulting GWs for different values of $\sigM$, which
refers to the value of $\sigma$ at the initial time of the simulation.
The spectrum of GW energy density fluctuates around an approximately statistically
steady spectrum.
We show the spectra averaged over oscillations, and the shaded region corresponds
to the maxima and minima of the oscillations at every wave number.
The magnetic field spectrum is shown at the initial time of the simulation,
when the magnetic energy density has its maximum, and it decays for
later times.

\begin{figure*}[t!]\begin{center}
%\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}
%AR: modified and added figures
%AB: It would be good for me to be able to check which runs those figures
%AB: refer to. Therefore I keep the 4 lines above for now. Or do you have a list?
%AB: Maybe it is already clear from the table, where the directories are listed.
\includegraphics[width=.49\textwidth]{i_s01EM_EGW.pdf}
\includegraphics[width=.49\textwidth]{i_s03EM_EGW.pdf}
\includegraphics[width=.49\textwidth]{i_s05EM_EGW.pdf}
\includegraphics[width=.49\textwidth]{i_s1EM_EGW.pdf}
\end{center}\caption[]{
Magnetic and GW energy (solid lines) and helicity (dashed lines)
spectra for the runs with an initial magnetic field, for different values of
$\sigM$.
The magnetic spectra are shown at the initial time of the simulation, and
the GW spectra are the saturated ones (averaged over times $t > 1.01$),
while the shaded region shows the variation over time between the maxima and minima.
Positive values of the helicity spectra (magnetic and GW) are shown in red, while
negative values in blue.
The dotted lines are added for comparison with power laws at different wave numbers.
}\label{pspecm_hel_initial}\end{figure*}

%We see that in all cases, the causal $k^4$ magnetic energy spectrum
%\cite{Durrer:2003ja} is established.
%AR: This would refer to the cases with forced magnetic field, so moved below
%This implies a flat spectrum for the GW energy, as shown from numerical
%simulations in Ref.~\cite{Pol:2019yex}, and obtained theoretically for
%a constant-in-time stress spectrum in Ref.~\cite{Brandenburg:2019uzj}.
%AR: slightly modified, and added comment on Gaussianity
%AB: ok
We observe a flat spectrum for the GW energy, as shown from numerical
simulations in Ref.~\cite{Pol:2019yex}, and obtained theoretically for
a constant-in-time stress with a white noise spectrum ($\sim k^2$).
The latter is expected for a Gaussian magnetic field with a $k^4$ Batchelor
spectrum \cite{Brandenburg:2019uzj}.
%AR: Axel, does the result in your paper with Boldyrev depend on helicity?
%AR: I mean the analytical calculation.
%AR: added text
%AB: yes, a little; this is the green lines in that paper.
The small deviations from an exact flat spectrum might be due to the small
number of points in wave number space when we reach the largest scales of the
simulation, as well as to the oscillations over time.

The magnetic helicity spectrum shows the same power law scalings as the spectrum
of the magnetic energy density in the inertial range, which correspond to
Kolmogorov type spectra ($k^{-5/3}$).
In the subinertial range, we observe that as we decrease the partial helicity,
the helical spectrum becomes slightly steeper than the $k^4$ Batchelor spectrum
observed in the magnetic energy density, being both spectra identical in the fully
helical case, as expected.
Note that the magnetic helicity spectra are multiplied by $k/2$ to take into
account the realizability condition.
It is also for large scales (small wave numbers) that we observe a decay of the
GW helical spectrum, $\HGW(k)$, with respect to the flat spectrum observed for
$\EGW(k)$.
As we decrease the fractional helicity of the magnetic field, we observe this
decay to be at approximately twice the smallest wave number of the magnetic
field (expected from the convolution).
This could be due to the steeper slope of the helical magnetic spectrum at low
wave numbers, impacting the helical spectrum of GWs.
Hence, we can expect the helicity spectrum of GWs to decay with respect to the
flat spectrum in the subinertial range, and omit the smallest wave number of the
simulation, $k_0$, since the correct computation of this mode would require the
computation of the magnetic field at $k_0/2$.
We defer further discussion of GW polarization $\PPGW(k) = \HGW(k)/\EGW(k)$ to the next
section.

\subsection{Runs with forced magnetic field at initial times}

In \Figs{pspecm_hel_forc001}{pspecm_hel_forc}, we compare the GW energy and
helicity spectra of the magnetic field and the resulting GWs for the runs in which we 
drive the magnetic field for a short duration $\delta t = \tmax - 1 = 0.1$
(i.e., 10\% of the Hubble time).
In \Fig{pspecm_hel_forc001}, we show runs that are almost non-helical
($\sigma = \pm 0.01$), and in \Fig{pspecm_hel_forc} we show runs with larger
$\sigma$\footnote{%
For the case with forced magnetic fields, $\sigma$ corresponds
to the helicity of the forcing term, see \Eq{forcing}, while $\sigM$ refers to
the helicity of the driven magnetic field at $\tmax$.
The two values are very similar, but not identical in general (see \Tab{runs}).},
up to the fully helical case (which is negative in this case).
We show that magnetic fields with negative helicity drive GW with the same
GW energy density, but negative chiral spectrum.

The magnetic spectra are shown at the time $\tmax$ when the magnetic energy
density reaches its maximum value $\OmMmax$, and when we switch off the forcing
term.
At earlier terms, the magnetic spectrum corresponds to a spike around the forcing
wave number $\kf$, which is distributed to a turbulent spectrum due to the MHD
dynamical evolution of the magnetic field.
We see that in all cases, the causal $k^4$ magnetic energy spectrum
\cite{Durrer:2003ja} is established.
The spectra are averaged over the time interval $1.1\leq t\leq 1.11$,
which is just after the maximum magnetic energy has been reached and
the GW energy begins to fluctuate around an approximately statistically
steady state.
For $k>k_*$, the slope of the magnetic energy spectrum is
steeper than a Kolmogorov-type spectrum ($k^{-5/3}$).
This is because of the finite time driving during the rather short time
interval, $1\leq t\leq 1.1$.
The consequences of the finite forcing time are the appearance of a smoothed
spike around $\kf$ which has not completely disappeared, and the steeper
slopes, especially in the inertial range.
This effect seems to be enhanced with helicity.
For the almost non-helical runs, we observe $\sim k^{-1.9}$ right below the peak,
which then turns to $k^{-8/3}$, and smoothly turns to a steeper slope down to the
diffusive scale.
For the cases with $|\sigma| = 0.3, 0.5, 0.7, 1$, the slopes right below the peak
are $k^{-2.2}$, $k^{-2.67}$, $k^{-3.2}$, $k^{-3.4}$, respectively.
Again, this turns to steeper scales down to the the diffusive scales.
If the driving was continued over a long time interval
(long enough for the magnetic field to be
processed, i.e., $\tmax \sim 1$) we would recover
the Kolmogorov spectrum \cite{Kahniashvili:2020jgm}.

\begin{figure*}[t!]\begin{center}
\includegraphics[width=.49\textwidth]{f_s001EM_EGW.pdf}
\includegraphics[width=.49\textwidth]{f_s001_negEM_EGW.pdf}
\end{center}\caption[]{
Magnetic and GW energy (solid lines) and helicity (dashed lines)
spectra, similar to \Fig{pspecm_hel_initial}, for runs with an almost
non-helical forcing ($\sigma = \pm 0.01$) for short times ($1\leq t \leq 1.1$).
}\label{pspecm_hel_forc001}\end{figure*}

\begin{figure*}[t!]\begin{center}
\includegraphics[width=.49\textwidth]{f_s03EM_EGW.pdf}
\includegraphics[width=.49\textwidth]{f_s05EM_EGW.pdf}
\includegraphics[width=.49\textwidth]{f_s07EM_EGW.pdf}
\includegraphics[width=.49\textwidth]{f_s1_negEM_EGW.pdf}
\end{center}\caption[]{
Magnetic and GW energy (solid lines) and helicity (dashed lines)
spectra, similar to \Fig{pspecm_hel_initial}, for runs with forcing
for short times ($1\leq t \leq 1.1$), for different values of $\sigma$.
}\label{pspecm_hel_forc}\end{figure*}

The magnetic helicity spectrum around the peak $k=\kf$ is nearly saturated
for all values of $\sigma$ considered from 0.3 to 1.
For $\sigma = \pm 0.01$, 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, especially at high wave numbers.
We also see a systematic sign flip at higher $k$ both in $\HM$ and
in $\HGW$, which also occurs for the cases with $\sigma \geq 0.3$.
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.
The helicity spectrum $k \HM$ is steeper than the magnetic spectrum,
as opposite to the case of the runs with an initial magnetic field.
This is observed for all values of $\sigma$.
However, for the full helical case, the difference in slope is smaller, and it
becomes larger for smaller fractional helicity.

\subsection{Spectrum of polarization}

In \Fig{ppol_comp}, we plot the GW polarizations, $\PPGW (k)$, and
$\PPh(k)$, for the two types of simulations.
We see that $\mathcal{P}_{\rm GW}(k)$ reaches $\pm1$ at
$k\approx2k_*\approx1200$ when $\sigma=\pm1$.
This is larger than what was found in previous analytic predictions, see 
Ref.~\cite{Kisslinger:2015hua}, and recently used in Ref.~\cite{Ellis:2020uid}
in relation to detectability with LISA.
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
$\PPGW (k)=0$, especially in the case of an initial given magnetic field,
due to the deviation of the helical GW spectrum at low wave numbers.

\begin{figure*}[t!]\begin{center}
\includegraphics[width=.49\textwidth]{PGW_initial_sigma.pdf}
\includegraphics[width=.49\textwidth]{PGW_forcing_sigma.pdf}
\includegraphics[width=.49\textwidth]{Ph_initial_sigma.pdf}
\includegraphics[width=.49\textwidth]{Ph_forcing_sigma.pdf}
\end{center}\caption[]{
Spectrum of the polarization degree for the different values of $\sigma$,
for the types of simulations:
runs with a given initial field (panels (a) and (c)), and
runs with a forced magnetic field for a short time $1 \leq t \leq 1.1$
(panels (b) and (d)).
The solid line represents the saturated value of the polarization (i.e.,
averaged over times in the oscillatory regime), and the shaded regions
are the maximum and minimum values of the polarization over oscillations.
%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{Pol:2019yex},
%for $\sigma=1$, 0.5, and 0.1.
}\label{ppol_comp}\end{figure*}

The GW polarization integrated over all wave numbers, $\PPGW$, as a function of the
magnetic fractional helicity, $\PPM$, is shown in \Fig{pppol1}.
The value of the integrated magnetic fractional helicity $\PPM$ as a function
of time is shown in \Fig{PM_vs_t}.
In the runs with an initial given magnetic field, $\PPM$ stays constant
for a short time, $t - 1 \sim 10^{-2}$, which is similar to the time that the
GW spectrum takes to enter the stationary regime.
Hence, the value of $\PPM$ does not change while the GW spectrum is established.
However, for the case in which the magnetic field is driven for a time
$\tmax - 1 = 10^{-1}$, the fractional helicity $\PPM$ has changed more significantly
when the GW spectrum is established, such that the variation of helicity
affects the GW polarization.

\begin{figure}[t!]\begin{center}
\includegraphics[width=.51\textwidth]{PGW_vs_PM_sigmaM.pdf}
\end{center}\caption[]{
GW polarization versus magnetic polarization.
}\label{pppol1}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=.49\textwidth]{PM_vs_t_initial.pdf}
\includegraphics[width=.49\textwidth]{PM_vs_t_forcing.pdf}
\end{center}\caption[]{
Time evolution of the fractional helicity of the magnetic field, for the different
values of $\sigma$, for (a) the runs with an initial given magnetic field, and
(b) those with a forced magnetic field.
}\label{PM_vs_t}\end{figure}

%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.
%AR: I think that the tangent data fit works when the polarization of the
%AR: magnetic field is taken at the time \tmax.
%AR: However, this does not correspond to the initial \sigma value of the forcing
%AR: (see the new plots). Using the values of the initial sigma, this seems
%AR: linear relation.
%%%%%%%%%%%%%AR: text about Beltrami field commented out for now.
%AR: Once I am done with the appendix, I will review this text, since the
%AR: plots have been changed
%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.
%%AR: removed tildes
%analytical solution for both $\PPGW$ and
%$\PPM$; see Appendix~A.
%The resulting dependence of $\PPGW$ on
%$\PPM$ 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}$.
%%%%%%%%%%%%%AR: text on Beltrami field to be updated!!

\subsection{Analytical prediction of the spectrum of polarization}

The GW degree of circular polarization of signals produced by primordial
MHD turbulence has been estimated in Refs.~\cite{Kahniashvili:2005qi,
Kisslinger:2015hua}.
In this section, we use their model to predict the spectrum of GW
polarization $\PPh$ and compare to the numerical results, to explore
the validity of their model and the impact of the assumptions made.
In first place, we define the unequal time correlation (UTC) function of the
magnetic field as
\begin{equation}
\bra{\tilde B_i (\kk', t_1) \tilde B^*_j (\kk, t_2)} = (2 \pi)^6 \delta^3
(\kk - \kk') \left(P_{ij} \frac{\FM (k, t_1, t_2)}{4 \pi k^2} +
i \varepsilon_{ijl} k_l \frac{\GM(k, t_1, t_2)}{8 \pi k^2}\right),
\end{equation}
where we recover the equal time correlation (ETC) of the magnetic energy 
density and helicity, defined in \Eq{equaltime}, when $t_1=t_2=t$.
Following Ref.~\cite{Kahniashvili:2005qi}, which assumes stationary 
turbulence, we express the spectral functions of the UTC of the freely 
decaying turbulent source as a function of only the time
difference, modelled as $\FM(k, t_1, t_2) =
\FM(k, t_1, t_1 + \delta t) \approx \EM(k, t_1) D_1 (\delta t)$.
Similarly, the helical contribution to the UTC spectrum is
$\GM(k, t_1, t_2) = \GM(k, t_1, t_1 + \delta t)
\approx \HM(k, t_1) D_2 (\delta t).$
The functions $D_1(\delta t)$, $D_2(\delta t)$ are monotonically decreasing 
functions.
To characterize the two types of turbulence considered in the present work,
$t_1$ is considered to be the time when the turbulence free decay starts.
This corresponds to the initial time ($t_1=1$) for the cases
with an initial magnetic field, and to the time at which the forcing term
is switched off ($t_1=\tmax$) otherwise.
Hence, $D_1(0)=D_2(0)=1$.

Assuming that the duration of the turbulence sourcing is short
(i.e., $\tau = t_{\rm fi} - 1 \ll 1$), such that the expansion of the Universe
can be neglected, and after averaging over the time oscillations of the
source and the solution to the GW equation, the functions $S_h(k,t_1)$,
$A_h(k,t_1)$ can be obtained as
\cite{Kahniashvili:2005qi, Kosowsky:2001xp},
\begin{align}
    S_h(k)= A \frac{\tau}{k^2} \int \dd \ln p_1 \int 
    \dd \ln p_2  \bar \Theta \Big[ & (1 + \gamma^2) (1 + \beta^2)
    \EM(p_1) \EM (p_2) \nonumber \\
    & + 4 \gamma \beta \HM(p_1) \HM (p_2) \Big],
    \label{Sh_an}\\
    A_h(k)=  2 A \frac{\tau}{k^2} \int \dd \ln p_1 \int 
    \dd \ln p_2  \bar \Theta \Big[ & (1 + \gamma^2) \beta
    \EM(p_1) \HM (p_2) \nonumber \\
    & + (1 + \beta^2) \gamma \HM(p_1) \EM (p_2) \Big],
    \label{Ah_an}
\end{align}
where $\bar \Theta = \theta(p_1 + p_2 -k) \theta (p_1 + k - p_2)
\theta (p_2 + k - p_1)$, being $\theta$ the Heaviside step function,
$\gamma = (k^2 + p_1^2 - p_2^2)/(2 k p_1)$,
$\beta=(k^2 + p_2^2 - p_1^2)/(2 k p_2)$,
and $A$ is a constant that we omit, since we are interested
in the spectral shapes, and we will use \Eqs{Sh_an}{Ah_an}
to compute the polarization, which does not depend on
$A$\footnote{\EEqs{Sh_an}{Ah_an} correspond to
Eqs.~(10) and (11) of Ref.~\cite{Kahniashvili:2005qi}, in which
$P_{\rm S} (k) = 2\pi^2 \EM(k)/k^2$,
$P_{\rm A} (k) = \pi^2 \HM (k)/k$,
$H(k) = 2 \pi^2 S_h(k)/k^2$,
${\cal H} (k) = 2 \pi^2 A_h(k)/k^2$,
and their value $A$ contains $\tau/k^2$ and it differs by constant
coefficients due to the normalization we use in \Eq{GW}.
References~\cite{Kisslinger:2015hua,Ellis:2020uid} use the notation 
${\cal I}_{\rm S}(k)=P_{\rm S} (k)$, and ${\cal I}_{\rm A}=
P_{\rm A} (k)$.}.
Previous analytical assumptions \cite{Kahniashvili:2005qi,
Kisslinger:2015hua,Ellis:2020uid} consider two types of turbulence:
\begin{itemize}
    \item  Helical Kolmogorov (HK) turbulence driven by magnetic
    energy dissipation at small scales, resulting in spectral powers
    $n_{\rm S}=-5/3$, and $n_{\rm A}=-8/3$\footnote{We refer here
    to the spectral slopes of $\EM(k) \propto k^{n_{\rm S}}$, and
    $\half k \HM(k) \propto k^{n_{\rm A}}$, while Ref.~\cite{Kahniashvili:2005qi} use spectral slopes of
    $P_{\rm S}$, and $P_{\rm A}$, which are divided by $k^2$.},
    \item Turbulence determined by helical transfer (HT) and helicity
    dissipation at small scales, which results in $n_{\rm S} =
    n_{\rm A} = -7/3$, based in Ref.~\cite{MC96},
\end{itemize}
and they use power law spectra $\EM(k) \sim k^{n_{\rm S}}$,
$\half k \HM (k) \sim h k^{n_{\rm S}}$ in the range
$\kf < k < k_{\rm D}$, with $h$ being the fraction of helicity 
dissipation, and $k_{\rm D}$ is the dissipative scale,
given by the fluid viscosity, or the magnetic diffusivity.

We extend this analytic approach to consider a broken power
law with a subinertial Batchelor spectrum $k^4$ below $\kf$,
and modify the HT spectrum to $n_{\rm S}=n_{\rm A}=-5/3$,
corresponding to the spectral slopes of the stochastic magnetic
fields that we use in our numerical simulations (see 
\Fig{pspecm_hel_initial}) and that is based in previous MHD
simulations applied to cosmological phase transitions
\cite{Brandenburg:2017neh}.
\FFig{PGW_analytical} shows the resulting polarization degree
$\PPh(k)$ for the HT and HK types of turbulence.
The inclusion of the subinertial range leads to an increase on
polarization at wave numbers right below the peak, which is a more
realistic scenario, when compared to the numerical results,
especially for the case of HK turbulence.

\begin{figure}[t!]\begin{center}
\includegraphics[width=.44\textwidth]{P_analytical_Kol.pdf}
\includegraphics[width=.54\textwidth]{P_analytical_Mo.pdf}
\end{center}\caption[]{
Polarization spectra $\PPh (k)$ for HK turbulence with
$n_{\rm S} = -5/3$, $n_{\rm A} = -8/3$ in the subinertial
range (blue), and for HT turbulence using $n_{\rm S} = n_{\rm A} = 
-5/3$ (red), using a single power law in the inertial range (dashed),
and using a broken power law with a Batchelor spectrum in the
subinertial range (dotted), for fractional helicities
$h=0.1, 0.3, 0.5, 0.8, 1$.
We use $\kf=600$ for comparison with the numerical
simulations.
}\label{PGW_analytical}\end{figure}

The model described above assumes that the partial helicity at $k=\kf$ is $h$,
and describes the magnetic energy, and helical spectra using power laws.
The assumption that the slopes are the same along the inertial range is
accurate in the case with an initial magnetic field; see
\Figs{pspecm_hel_initial}{PGW_analytical_comp_HT}.
In the latter, we compare the analytical results obtained from
HT turbulence (with modified slopes and extended to broken power law)
with those of our numerical simulations.
Note that the extension of the turbulent spectra to $k < \kf$ allows to
get a more accurate position of the polarization spectral peak.
In general, we observe good agreement of the analytical results with the
polarization $\PPGW (k)$ computed from the simulations, while the spectrum
$\PPh(k)$ agrees well at low and intermediate $k$, but decays with respect
to the analytical models for large $k$.
However, in the case in which the magnetic field is forced for a short
duration of time, we observe in \Fig{pspecm_hel_forc} that the energy
and helicity spectra are similar around the spectral peak (where the magnetic
field is maximally helical),
while the helical spectrum starts decaying with a steeper slope
than the magnetic spectrum at some wave number $k_{\rm hel}$ which increases
with the helicity of the forcing.
This behavior is not captured by previous analytical models.
In \Fig{PGW_analytical_comp_HK} we compare the results of the polarization
spectra $\PPh(k)$, and $\PPGW(k)$ obtained from HK turbulence (with modified
slopes, and extended to broken power law) with those obtained from the
numerical simulations.
In addition, to take into account the deviations from the assumption
of constant slopes of the helical spectrum, we compute the polarization
using \Eqs{Sh_an}{Ah_an} to integrate over the numerically computed
spectra $\EM (k)$, $\HM (k)$ of the turbulence source.
We observe that the latter gives an accurate approximation of $\PPGW(k)$
to the numerical results, while $\PPh$ shows again some deviations at large
$k$.
The polarization spectrum shows features of both HK and HT types of turbulence
in the spectrum.
As we have mentioned, around the spectral peak, the numerical simulations
show that the magnetic field is maximally helical in all the considered
runs with $\sigma_{\rm M} > 0.3$ ($h=0.55$), which is better represented
by the HT spectrum (with $h=1$), while in the inertial range the different
slopes lead to a decrease of the polarization towards a HK spectrum, which
is more noticeable as we decrease the fractional helicity.
However, since the polarization spectrum is obtained via convolution over wave 
numbers, a result linearly combining both types of turbulence is not possible.

\begin{figure}[t!]\begin{center}
\includegraphics[width=.49\textwidth]{PGW_an_i_lows.pdf}
\includegraphics[width=.49\textwidth]{PGW_an_i_highs.pdf}
\end{center}\caption[]{
Polarization spectra $\PPGW (k)$ (red dots) and $\PPh (k)$ (blue dots)
obtained from the numerical simulations, compared to $\PPh (k)$ computed from
the analytical model (extended to a broken power law) using HT type of
turbulence (dotted thin), and from the analytical integrals;
see \Eqs{Sh_an}{Ah_an}), using the numerical spectrum of the turbulence (solid).
}\label{PGW_analytical_comp_HT}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=.4\textwidth]{PGW_an_f_s03.pdf}
\includegraphics[width=.4\textwidth]{PGW_an_f_s05.pdf}
\includegraphics[width=.4\textwidth]{PGW_an_f_s07.pdf}
\includegraphics[width=.4\textwidth]{PGW_an_f_s1.pdf}
\end{center}\caption[]{
Polarization spectra $\PPGW (k)$ and $\PPh (k)$ obtained from
the numerical simulations, compared to $\PPh (k)$ obtained from
the analytical using HK type of turbulence (solid),
and obtained from the analytical integral, but using the
numerical spectrum of the turbulence (dotted thin).
}\label{PGW_analytical_comp_HK}\end{figure}

We have confirmed using numerical simulations that in the case of an
initial magnetic field with a Kolmogorov spectrum for both the magnetic
energy density and helicity, the analytical model previously considered in 
Refs.~\cite{Kahniashvili:2005qi,Kisslinger:2015hua,Ellis:2020uid}
using HK turbulence gives an accurate prediction of the degree of circular
polarization.
However, when we consider the scenario in which the magnetic field is generated
via MHD forcing for a short amount of time, and non-linear interactions appear,
allowing to have different spectral slopes at different scales, previous
analytical estimates underpredict considerably the polarization degree
peak and fails to predict the appropriate shape, which is not given by
either the HK nor the HT models, but presents non-linearly combined features of
both.
We observe a helical inverse cascade that produces larger degree of circular
polarization at large scales.

\section{Gravitational radiation from the electroweak phase
transition}

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

\begin{figure}[t!]\begin{center}
\includegraphics[width=.49\textwidth]{OmGW_initial_detectors_OmM01.pdf}
\includegraphics[width=.49\textwidth]{OmGW_forcing_detectors_OmM01.pdf}
\includegraphics[width=.49\textwidth]{OmGW_initial_detectors_OmM005.pdf}
\includegraphics[width=.49\textwidth]{OmGW_forcing_detectors_OmM005.pdf}
\end{center}\caption[]{GW spectrum $h_0^2 \OmGW$ for
signals produced at the EWPT from magnetic fields
with $\OmM = 0.1$ for the driven case  (a) and
the case with initial condition (b),
and with $\OmM = 0.05$ for the driven case  (c) and
the case with initial condition (d).
}\label{OmGW_detectors}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=.49\textwidth]{XiGW_initial_detectors_OmM01.pdf}
\includegraphics[width=.49\textwidth]{XiGW_forcing_detectors_OmM01.pdf}
\includegraphics[width=.49\textwidth]{XiGW_initial_detectors_OmM005.pdf}
\includegraphics[width=.49\textwidth]{XiGW_forcing_detectors_OmM005.pdf}
\end{center}\caption[]{Helical GW
spectrum $h_0^2 \XiGW$ for
signals produced at the EWPT from magnetic fields
with $\OmM = 0.1$ for the driven case  (a) and
the case with initial condition (b),
and with $\OmM = 0.05$ for the driven case  (c) and
the case with initial condition (d).
}\label{XiGW_detectors}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=.6\textwidth]{OmGW_efficiency_vs_t_more2.pdf}
\end{center}\caption[]{
GW energy density $\OmGW$ compensated by $\Omega_i^2/\kf^2$, for $i=$M (magnetic),
K (kinetic), for the runs with initial
given magnetic field ('initial'), and for the runs with a forced magnetic field
('forcing (short)').
Added for comparison are the runs in Ref.~\cite{Kahniashvili:2020jgm}
('forcing (long)'), in which the magnetic field is forced for longer times
($\tmax = 1$), and the runs in Ref.~\cite{Pol:2019yex}, which contain runs
with initial given magnetic field, with forced magnetic field, at $\kf=10, 100,
1000$, and runs for acoustic turbulence ('acoustic').
}\label{OmGW_efficiency}\end{figure}

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}
%pspecm_scl5_hel_sig, pspecm_scl5_hel_sig
\includegraphics[width=.49\textwidth]{OmGW_initial_sigma.pdf}
\includegraphics[width=.49\textwidth]{OmGW_forcing_sigma.pdf}
\includegraphics[width=.49\textwidth]{XiGW_initial_sigma.pdf}
\includegraphics[width=.49\textwidth]{XiGW_forcing_sigma.pdf}
\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*}

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.

\section*{Acknowledgements}

%AB: adapted from earlier paper. Please check.
Support through the Swedish Research Council (Vetenskapsr{\aa}det),
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.

%\newpage
%AB: commented out for now; or is it required by JCAP?

\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
%TK agree
%\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}
%AB: commented out this section heading.

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}
%TK to be added, right?

%r e f
\begin{thebibliography}{}

\bibitem{Brandenburg:2017neh}
A.~Brandenburg, T.~Kahniashvili, S.~Mandal, A.~Roper Pol, A.~G.~Tevzadze and T.~Vachaspati,
{\em Evolution of hydromagnetic turbulence from the electroweak phase transition},
Phys. Rev. D \textbf{96} (2017) no.12, 123528.

\bibitem{Neronov:1900zz}
A. Neronov and I. Vovk\ysci{2010}{328}{73}
{Evidence for Strong Magnetic Fields from Fermi Observations of TeV Blazars}

\bibitem{Vachaspati:2020blt}
T.~Vachaspati,
{\em Progress on Cosmological Magnetic Fields},
arXiv:2010.10525 (2020).

\bibitem{Vachaspati:1991nm}
T.~Vachaspati,
{\em Magnetic fields from cosmological phase transitions},
Phys. Lett. B {\bf 265}, 258-261 (1991)

\bibitem{Cornwall:1997ms}
J.~M.~Cornwall,
{\em Speculations on primordial magnetic helicity},
Phys. Rev. D \textbf{56}, 6146 (1997)

\bibitem{Vachaspati:2001nb}
T.~Vachaspati,
{\em Estimate of the primordial magnetic field helicity},
Phys. Rev. Lett. {\bf 87}, 251302 (2001)

\bibitem{GarciaBellido:2002aj}
J.~Garc\'ia-Bellido, M.~Garc\'ia P\'erez and A.~Gonz\'alez-Arroyo,
{\em Symmetry breaking and false vacuum decay after hybrid inflation},
Phys. Rev. D{\bf 67}, 103501 (2003).

\bibitem{Sharma:2019jtb}
R.~Sharma, K.~Subramanian and T.~R.~Seshadri,
{\em Gravitational wave generation in a viable scenario of inflationary magnetogenesis},
%Phys. Rev. D \textbf{101} (2020) no.10, 103526.
%AB: I don't think we need to give the issue number here
Phys. Rev. D \textbf{101} (2020), 103526.

\bibitem{OF21}
S. Okano and T. Fujita\yjcap{2021}{03}{026}
{Chiral gravitational waves produced in a helical magnetogenesis model}

\bibitem{Kahniashvili:2005qi}
T. Kahniashvili, G. Gogoberidze, and B. Ratra\yprl{2005}{95}{151301}
{Polarized Cosmological Gravitational Waves from Primordial Helical Turbulence}

\bibitem{Kisslinger:2015hua}
L.~Kisslinger, T.~Kahniashvili,
{\em Polarized gravitational waves from cosmological phase transitions},
Phys. Rev. D {\bf 92}, 043006 (2015).

\bibitem{Kahniashvili:2020jgm}
T.~Kahniashvili, A.~Brandenburg, G.~Gogoberidze, S.~Mandal and A.~Roper Pol,
{\em Circular polarization of gravitational waves from early-Universe helical turbulence},
Phys. Rev. Res. \textbf{3}, no.1, 013193 (2021).
%doi:10.1103/PhysRevResearch.3.013193
%[arXiv:2011.05556 [astro-ph.CO]].

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

\bibitem{MC96}
S. S. Moiseev and O. G. Chkhetiani\yjetp{1996}{83}{192}
{Helical scaling in turbulence}

\bibitem{Ellis:2020uid}
J.~Ellis, M.~Fairbairn, M.~Lewicki, V.~Vaskonen and A.~Wickens,
{\em Detecting circular polarisation in the stochastic gravitational-wave background from a first-order cosmological phase transition},
JCAP {\bf 10} (2020) 032.

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

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

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

\bibitem{Orlando:2020oko}
G.~Orlando, M.~Pieroni and A.~Ricciardone,
{\em Measuring Parity Violation in the Stochastic Gravitational Wave Background with the LISA-Taiji network},
[arXiv:2011.07059 [astro-ph.CO]].

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

\bibitem{Pol:2019yex}
A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky\yprd{2020}{102}{083512}
{Numerical simulations of gravitational waves from early-universe turbulence}

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

%\bibitem{Sub16}
%K. Subramanian\yrpp{2016}{79}{076901}
%
%\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{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{Brandenburg:1996fc}
A. Brandenburg, K. Enqvist and P. Olesen\yprd{1996}{54}{1291}
{Large-scale magnetic fields from hydromagnetic turbulence in the very early universe}

\bibitem{Durrer:2013pga}
R. Durrer and A. Neronov\yanar{2013}{21}{62}
{Cosmological Magnetic Fields: Their Generation, Evolution and Observation}

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

\bibitem{Durrer:2003ja}
R. Durrer and C. Caprini\yjour{2003}{J.\ Cosmol.\ Astropart.\ Phys.}{0311}{010}
{Primordial Magnetic Fields and Causality}

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

%\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{Caprini:2003vc}
C. Caprini, R. Durrer, and T. Kahniashvili\yprd{2004}{69}{063006}
{The Cosmic Microwave Background and Helical Magnetic Fields:
The Tensor Mode}

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

\bibitem{Kosowsky:2001xp}
A.~Kosowsky, A.~Mack and T.~Kahniashvili,
``Gravitational radiation from cosmological turbulence,''
Phys. Rev. D \textbf{66}, 024030 (2002).

%\bibitem{CCFTV15}
%W. Chen, B. D. Chowdhury, F. Ferrer, H. Tashiro, T. Vachaspati\ymn{2015}{450}{3371}
%

%\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}
