%\documentclass[superscriptaddress,aps,pra,a4paper,nofootinbib,preprint]{revtex4}
%AB: preprint -> reprint (for now)
\documentclass[superscriptaddress,aps,prd,nofootinbib,reprint]{revtex4-2}
%\documentclass[superscriptaddress,aps,preprintnumbers,amsmath,amssymb,prl,nofootinbib]{revtex4}

%%%%%%%%%%%% packages 
\usepackage{amsfonts,amsmath,amssymb}
\usepackage{bm,tensor}
\usepackage{graphicx}
\usepackage{braket}
\usepackage{ifpdf}
\usepackage[utf8]{inputenc}
\usepackage{color}
\usepackage{ulem}
\usepackage{comment}
% the color option for hyperref based on jcappub.sty
\usepackage[colorlinks=true
,urlcolor=blue
,anchorcolor=blue
,citecolor=blue
,filecolor=blue
,linkcolor=blue
,menucolor=blue
,linktocpage=true
,pdfproducer=medialab
,pdfa=true
]{hyperref}
%\usepackage[inline]{showlabels}


%%%%%%%%%%%% new command 
\graphicspath{{./fig/}{./png/}}
\input macros
\newcommand{\red}[1]{\textcolor{red}{#1}}
\newcommand{\blue}[1]{\textcolor{blue}{#1}}
\newcommand{\SH}[1]{\textcolor[rgb]{0,0.6,0.1}{#1}}


\newcommand{\KK}{\textcolor{cyan}}
\newcommand{\ckm}[1]{{\color{blue}#1}}
\newcommand{\HN}[1]{\textcolor{magenta}{#1}}

\newcommand{\df}{\text{d}}
\newcommand{\p}{\partial}
\newcommand{\Mpl}{{M_\text{pl}}}
%\newcommand{\GeV}{{\rm GeV}}
%AB: added space
\newcommand{\GeV}{\,{\rm GeV}}


\begin{document}


%%%%%%%%%%%%     title
\title{
Supplemental Material to ``Chiral Magnetohydrodynamics with Zero Total Chirality''
}

\date{\today}
\preprint{NORDITA-2022-0xx}

\author{Axel Brandenburg}
\affiliation{Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden}
\affiliation{The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden}
\affiliation{School of Natural Sciences and Medicine, Ilia State University, 0194 Tbilisi, Georgia}
\affiliation{McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA}

\author{Valerie Domcke}
\affiliation{CERN}
\affiliation{EPFL}

\author{Kohei Kamada}
\affiliation{Research Center for the Early Universe (RESCEU), Graduate School of Science,\\ The University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo 113-0033, Japan}\

\author{Kyohei Mukaida}
\affiliation{KEK}

\author{Kai Schmitz}
\affiliation{CERN}

\author{Jennifer Schober}
\affiliation{Laboratoire d'Astrophysique, EPFL, CH-1290 Sauverny, Switzerland}


\begin{abstract}
In this 
\end{abstract}

\maketitle

\KK{ %\section{Preface}
In this Supplemental Material, we present additional plots.
}

\subsection{Evolution of characteristic scales} 

%\KK{Main results should be the one in notes11.pdf.
%The investigation at notes12.pdf is a supporting discussion.
%The observations on the $\sigma$ dependence in notes13.pdf should be discussed.}
%AB: added these notes for now.

To make contact with earlier work and to facilitate parameter studies, we
begin by simulating cases with an initial magnetic field, whose spectrum
peaks initially at a certain wave number $k_{\rm I0}$ that can be varied.
The form of the spectrum is given by a $k^{-5/3}$ inertial range and a
$k^4$ subinertial range.
It has random phases and is fully helical.
The chiral chemical potential is initially assumed to be uniformly
distributed in space.
However, even for an initially uniform chiral chemical potential, there
is a specific length scale associated with the value of $\mu_5$ through
the wave number of the most unstable mode of the CPI, $k_5=\mu_{50}/2$.
It is then of interest to vary the separation between $k_5$ and $k_{\rm I0}$.
We consider two sets of runs: Set~A where we change $k_{\rm I}$ and keep
$k_5$ fixed, and Set~B where we change $k_{50}$ and keep $k_{\rm I0}$ fixed.
It both cases, of course, since we want to satisfy
$\mu_5+\mu_{\rm M}=\const$, we need to adjust the
amplitude of the initial magnetic field correspondingly.

With the start of the simulation, the helical random magnetic field, which is present
initially, drives turbulent motions through the Lorentz force.
Those motions are in approximate equipartition with the magnetic field at high
wave numbers; see \Fig{rspec_select_ub_1024a_mu10}, where we
compare kinetic and magnetic energy spectra at different times.

We have performed simulations for different values of $k_{\rm I0}$
and $\mu_{50}$.
% and compared the resulting values of $\eta_1$, $\eta_2$, and $\eta_3$; see \Tab{Tlist}.
It turns out that, even though $\mu_{\rm M}+\mu_5=0$, there is
still inverse cascading.
%We see this clearly in \Fig{rspec_select_W512a_mu20}, where the position
We see this clearly in \Fig{rspec_select_1024a_mu10}, where the position
of the magnetic peak, $k_{\rm I}(\eta)$, gradually moves to smaller values.
At the same time, the height of the peak decreases with decreasing
$k_{\rm I}$, following an approximate power law $\propto k^\beta$.
%for most of our cases.
The newly injected magnetic helicity from the CPI leads to a growth of the magnetic
field at large wave numbers.
It manifests itself mostly through the build-up of negative magnetic
helicity at high wave numbers.
At some point, we also see a gradual propagation of the
secondary peak $k_{\rm II}$ toward smaller $k$.
In \Fig{pkmax_comp_I}, we plot the evolution of $k_{\rm I}(\eta)$ and $k_{\rm II}(\eta)$
for the same run (Run~I).
We see that there are distinct times $\eta_1$, $\eta_2$, and, $\eta_3$,
where $k_\pm$ begins to decrease first rapidly (at $\eta=\eta_1$) and
later (at $\eta=\eta_2$) more slowly, approximately like $\eta^{-1/3}$.
The decay of $k_{\rm II}$ closely follows that of $k_\pm$.
The decay of $k_{\rm I}$, on the other hand, does not show the rapid
decay phase that we see in $k_\pm$ and $k_{\rm II}$, but turns
directly into the approximate $\eta^{-1/3}$ decay at $\eta=\eta_1$.
%\KK{KK: [In this example, the large-scale spectrum given as an initial condition
%does not show explicit inverse cascade.]}
%AB: yes, to be avoided

We see that $\eta_1$ decreases quadratically with increasing $\mu_{50}$;
see \Fig{peta1_vs_mu50}.
In other words, $k_{\rm II}$ starts to decline more rapidly when
$\mu_{50}$ is large.
The onset of the decline of $k_{\rm I}$ does not strongly depend on the
value of $\mu_{50}$.
Likewise, the time $\eta_3$ when the decay of $k_{\rm II}$ slows down,
does not strongly depend on $\mu_{50}$.
On the other hand, when we vary $k_{\rm I0}$, we find that the time
$\eta_2$ of the decline of $k_{\rm I}$ decreases quadratically with
increasing values of $k_{\rm I0}$.
Thus, similarly to the dependence of the decline of $k_{\rm II}$, we
find that the time $\eta_2$ of the value of $k_{\rm I}$ depends on the
initial value $k_{\rm I0}$.

Another example is shown in \Fig{rspec_select_1024a_mu10}, where the
initial magnetic peak $k_{\rm I}$ is fully included within the range
of available wave numbers.
We now see that both peaks show inverse cascading.

The evolution of the peaks of the spectrum can be summarized as follows. 
%i) After the onset of the run, the CPI induces the grow of the negative helicity modes at 
%AB: modified
(i) After the start of the run, the CPI induces a growth of the negative helicity modes at 
%the constant secondary peak $k_\mathrm{sec}$ until $\eta=\eta_1$, 
%AB: modified
the secondary peak $k_\mathrm{II}$, which stays constant until $\eta=\eta_1$, 
and then starts to decrease with time in a power law fashion, $k_\mathrm{II} \propto t^{-q_{\rm IIa}}$. 
(ii) The original large-scale spectrum is unchanged until some time $\eta=\eta_2$ and 
starts to decrease via an inverse cascade with $k_\mathrm{I} \propto t^{-q_{\rm I}}$.
%(iii) At time $\eta=\eta_3$, the decay of the secondary peak becomes milder with a smaller 
%AB: modified
(iii) At time $\eta=\eta_3$, the decay of the secondary peak becomes slower with a smaller 
%index, $k_\mathrm{sec} \propto t^{-\beta}$, $(\beta<\alpha)$.}
%AB: modified
index, $k_\mathrm{II} \propto t^{-q_{\rm IIb}}$, with $q_{\rm IIa}>q_{\rm IIb}$.

%\KK{For the large scale spectrum of the positive helicity mode,} 
%it turns out that $\beta$ is around 2, although in some cases it is
%closer to 3/2.
%The evolution with $\beta=3/2$ is somewhat similar to that of nonhelical
%turbulence (BKT15), where $\beta=1$ was found.
%When $\beta=3/2$, and the system is on the self-similarity line \KK{($p+2q=2$?~\cite{}[Olesen?])}, we
%have $q=4/9$ and $p=10/9$.

%\begin{table}[htb]\caption{
%}\vspace{12pt}\centerline{\begin{tabular}{ccccccc}
%Case& $k_1$ &$\mu_{50}$& $\nu$ & $N$ \\
%\hline
%  O &  0.1  &     5    &$5\times10^{-4}$&  512 \\% 512b_mu5
%  I &  0.05 &    10    &$2\times10^{-4}$& 1024 \\% 1024a_mu10
% II &  0.1  &    20    &$2\times10^{-4}$& 1024 \\% 1024a_mu20
%III &  0.2  &    50    &$2\times10^{-4}$&  512 \\% 512a_mu50
% IV &  0.1  &    50    &$5\times10^{-4}$& 1024 \\% 1024c_mu50
%\label{Ttimescale}\end{tabular}}\end{table}

%\begin{table}[htb]\caption{
%}\vspace{12pt}\centerline{\begin{tabular}{ccccccc}
%Case& $k_1$ & $\kf$ &$\mu_{50}$& $\nu$ & $N$ \\
%\hline
%  I &  0.05 &  1  & 10 &$2\times10^{-4}$& 1024 \\% 1024a_mu10
%  A &  0.05 & 0.5 & 10 &$2\times10^{-4}$& 1024 \\% 1024a_mu10_k05_1024proc2
%  B &  0.05 & 0.2 & 10 &$2\times10^{-4}$& 1024 \\% 1024a_mu10_k02_1024proc
%  C &  0.02 & 0.2 & 10 &$2\times10^{-4}$& 1024 \\% 1024b_mu10_k02_1024proc
%\label{Ttimescale2}\end{tabular}}\end{table}

%\begin{figure}\begin{center}
%\includegraphics[width=\columnwidth]{p_comp_512e2_1em1_mu10}
%\end{center}\caption[]{
%Solid: 1024a and 1024a3, $(C_\sigma,C_\nu,D)=$(3, 3, 1e-5) to (3, 6, 1e-4).
%Dotted: (0, 6, 1e-4).
%}\label{p_comp_512e2_1em1_mu10}\end{figure}

%\begin{figure}\begin{center}
%\includegraphics[width=\columnwidth]{rspec_select_W512a_mu20}
%\end{center}\caption[]{
%Spectral evolution for $\eta=10$, 15, 22, 32, 46, 68, and 100.
%The value of $k_\pm$ decreases from $k=10$ at $\eta=10$
%to $k=0.4$ at $\eta=100$.
%The inset shows magnetic energy and helicity spectra at
%$\eta=3200$ when finite domain size effects may have played
%a role and the fractional magnetic helicity is now small
%at all wavenumbers.
%The arrows indicate the temporal evolution for the realistic
%case with $H_{\rm f}=10^{11}\GeV$.
%}\label{rspec_select_W512a_mu20}\end{figure}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{rspec_select_ub_1024a_mu10}
\end{center}\caption[]{
%JS: I think the y-label needs to be corrected and the kinetic energy spectrum should be defined in III.B
%Evolution of $k_\pm(\eta)$.
%AB: this is done in line 4 of II.B.
Magnetic energy (red lines) and kinetic energy (blue lines) spectra for
Run~I at times $\eta=216$, 317, 465, 682, 1000, 1468, 2155, 3163, 4642,
6813, 10000, 14,680, and 21,550.
Initially, the magnetic peak of positive magnetic helicity (red symbols)
is at $k=k_{\rm I}\approx1$, while the second peak of negative magnetic
helicity at $k=k_{\rm II}$ (blue symbols) appears only later.
}\label{rspec_select_ub_1024a_mu10}\end{figure}

\begin{figure}[t]\begin{center}
\includegraphics[width=\columnwidth]{rspec_select_1024a_mu10}
\end{center}\caption[]{
%Evolution of $k_\pm(\eta)$.
Magnetic energy (solid lines) and normalized helicity spectra
($k\HM(k)/2$, dotted lines with red and blue symbols for positive
and negative helicity spectra, respectively) for Run~I.
The times are $\eta=216$, 317, 465, 682, 1000, 1468, 2155, 3163, 4642,
6813, 10000, 14,680, and 21,550.
}\label{rspec_select_1024a_mu10}\end{figure}

%Initially, magnetic helicity grows at small scales only.
%Toward later times, progressively larger length scales are
%being affected; see \Fig{rspec_select_W512a_mu20}.
%At $\eta>100$, the finite size of the domain begins
%to affect the simulation.
%At much later times, the magnetic field is
%nonhelical over the entire spectrum.
%\Figs{pkmax_comp_lin1024a_mu10}{pkmax_comp_red_1024a_mu10} show
%the evolution of various characteristic wave numbers.
%%The results are summarized in \Tabs{Ttimescale}{Ttimescale2}.

%\begin{figure}\begin{center}
%\includegraphics[width=\columnwidth]{rspec_select_1024a_mu10}
%\end{center}\caption[]{
%Magnetic energy (solid lines) and normalized helicity spectra
%($k\HM(k)/2$, dotted lines with red and blue symbols for positive
%and negative helicity spectra, respectively) for Run~I.
%The times are $\eta=216$, 317, 465, 682, 1000, 1468, 2155, 3163, 4642,
%6813, 10000, 14,680, and 21,550.
%}\label{rspec_select_1024a_mu10}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pkmax_comp_I}
\end{center}\caption[]{
Comparison of $k_{\rm peak}$ (red), $k_\pm$ (orange),
and $k_{\rm sec}$ (blue) for Run~I.
The dashed-dotted line indicates $\eta^{-1/3}$.
}\label{pkmax_comp_I}\end{figure}

%\begin{figure}\begin{center}
%\includegraphics[width=\columnwidth]{pkmax_comp_red_1024a_mu10}
%\end{center}\caption[]{
%Other way around:
%Comparison of $k_{\rm peak}$ (red), $k_\pm$ (orange),
%and $k_{\rm sec}$ (blue) for Run~1024a3
%(solid lines) with Run~1024a4 (dashed lines).
%The dashed-dotted line indicates $\eta^{-2/3}$.
%}\label{pkmax_comp_red_1024a_mu10}\end{figure}

In \Fig{pEMxi_pq_runI}, we show the $pq$ diagram for Run~I.
We see that the solution evolves along the $\beta_{\rm M}=3/2$ line
during an extended time interval.
This is expected if there is conservation of the Saffman-like helicity
invariant \cite{Hosking+Schekochihin20}.
At late times, however, the system evolves again away from the
self-similarity line, which may be a consequence of the simulation
domain not being large enough.

\begin{figure}[t]\begin{center}
\includegraphics[width=\columnwidth]{pEMxi_pq_runI}
\end{center}\caption[]{
$pq$ diagram for Run~I.
%JS: Are the times given here the same as for the spectra above?
%JS: Maybe indicate what the solid diagonal line indicates. 
%AB: not yet, need to check
}\label{pEMxi_pq_runI}\end{figure}

\begin{figure}[t]\begin{center}
\includegraphics[width=\columnwidth]{pEMxi_pq_runIb}
\end{center}\caption[]{
$pq$ diagram for Run~Ib.
}\label{pEMxi_pq_runIb}\end{figure}

%------------------ end of notes 11 ----------------------------
%---------------------- notes 12 ----------------------------

\begin{table}[b]\caption{
}\vspace{12pt}\centerline{\begin{tabular}{lccccccc}
Run & $\kf$ & $\mu_{50}$ & $\eta_1$ & $\eta_2$ & $\eta_3$ & $\eta_5$ \\
\hline
IV  &  1 & 50 & 1.8&   75 &   50 &  50 \\  % 1024c_mu50 (2e-4)
II  &  1 & 20 &  9 &   75 &  120 & 200 \\  % 1024a_mu20 (2e-4)
I   &  1 & 10 & 30 &   80 &  160 & 200 \\  % 1024a_mu10 (2e-4)
O   &  1 &  5 & 35 &  200 &  145 & 200 \\  % 512b_mu5   (5e-4)
A   & 0.5& 10 & 30 &  360 &  360 & 300 \\  % 1024a_mu10_k05_1024proc2  (2e-4)
C   & 0.2& 10 & 30 & 2000 & 2000 &3000 \\  % 1024b_mu10_k02_1024proc   (2e-4)
\label{Ttimescale3}\end{tabular}}\end{table}

Note that for Run~IV, $\eta_2>\eta_3$.
while for Run~I, $\eta_2<\eta_3$.
In many other cases, $\eta_2\approx\eta_3$.
It is unclear whether there is really a meaningful and significant
difference between $\eta_2$ and $\eta_3$.
These are the times when both $k_{\rm I}$ and $k_{\rm II}$ tend
to decrease with the same power, i.e., approximately like
$\eta^{-1/3}$.

%/home/brandenb/tex/mhd/kohei/idl
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{peta1_vs_mu50}
\end{center}\caption[]{
Dependence of $\eta_1$ on $\mu_{50}$, showing a
$\eta_1\propto\mu_5^{-2}$ behavior for $\mu_5>10$
%JS: Which symbol is eta_2 and which is eta_3? 
%AB: added now in the figure itself.
(filled circles), while $\eta_2$ (plus signs) and
$\eta_3$ (asterisks) are nearly constant.
}\label{peta1_vs_mu50}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{peta2_vs_kf}
\end{center}\caption[]{
Dependence of $\eta_2$ on $\kf$, showing a
$\eta_2\propto\kf^{-2}$ behavior.
}\label{peta2_vs_kf}\end{figure}

%/home/brandenb/tex/mhd/kohei/idl
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{peta1_vs_mu50_comp}
\end{center}\caption[]{
Dependence of the compensated $\eta_1\mu_{50}^2/\sigma$
on $\mu_{50}$.
}\label{peta1_vs_mu50_comp}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{phmax_comp}
\end{center}\caption[]{
Comparison of $\mu_{\rm M}^\pm$ (red and blue lines, respectively) and $\mu_5$ (black lines)
for Run~I (solid lines), Run~II (dashed lines), and Run~A (dotted lines).
The orange dotted-dashed line shows the $\eta^{-2/3}$ decay law,
which matches the actual decay rather well..
}\label{phmax_comp}\end{figure}

We see from \Fig{phmax_comp} that $\mu_-$ is surprisingly small.
The imbalance between $|\muM^+|$ and $|\muM^-|$ suggests that there is
never a phase with a near-complete cancellation of $\muM^+$ and $\muM^-$.
We now want to address the question whether this is related to the value
of $\sigma^{-1}$ and $\nu$.
We therefore now perform simulations with different values; see
\Fig{phmax_comp_diff}.
It turns out that runs with smaller magnetic diffusivity
($\sigma^{-1}=10^{-4}$) result in an even larger imbalance, while those
with a larger diffusivity ($\sigma^{-1}=5\times10^{-4}$) have a smaller
imbalance.

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{phmax_comp_diff}
\end{center}\caption[]{
Comparison of $\mu_{\rm M}^\pm$ (red and blue lines, respectively) and
$\mu_5$ (black lines) for Run~II (solid lines, as in \Fig{phmax_comp}
with $\sigma^{-1}=2\times10^{-4}$), smaller diffusion (dashed
lines, $\sigma^{-1}=10^{-4}$), and larger diffusion (dotted lines,
$\sigma^{-1}=5\times10^{-4}$).
}\label{phmax_comp_diff}\end{figure}

\subsection{Evolution of $\aver{\mu_5}$ and ${\cal H}^\pm$}

%\KK{
The chirality of the system also evolves. Since the magnetic helicity spectrum 
can be divided into the large-scale positively helical part and small-scale negatively helical part, 
we introduce $\langle \mu_M^{+} \rangle$ for the large-scale spectrum 
and $\langle \mu_M^{-} \rangle$ for the small spectrum, which satisfy $\langle \mu_M^{+} \rangle-\langle \mu_M^{-} \rangle=\langle \mu_M\rangle$.
We find the typical evolution of $\langle \mu_5 \rangle$  and  $\langle \mu_M^{\pm} \rangle$ as follows. 
(i) $\langle \mu_5 \rangle$ and $\langle \mu_M^{+} \rangle$ stays constant until a time $\eta=\eta_4$
and exhibits a power law decay after that. 
(ii)  $\langle \mu_M^{-} \rangle$ grows in a power-law $\langle \mu_M^{-} \rangle \propto t^\epsilon$ 
until a time $\eta = \eta_5$ and decays as another power law  $\langle \mu_M^{-} \rangle \propto t^{-\delta}$.
We find that the absolute value of $\langle \mu_M^{-} \rangle$  is always much smaller than
those of $\langle \mu_5 \rangle$ and $\langle \mu_M^{+} \rangle$. See Fig.~\ref{phmax_comp}.
%}
Note that the decay of $\mu_5$ is not driven by the enhancement
of ${\cal H}^-$.


\KK{KK: [How much are we confident in the value of the scaling index? 
Can it be an intermediate value to the real scaling index? Or can it be the onset of the exponential decay?
It could depend on the quantities.]}

%------------------ end of notes 12 ----------------------------
%---------------------- notes 13 ----------------------------

%\subsection{Decay of $\mu_5$: continued to lower Reynolds numbers}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pcomp_mudiff}
\end{center}\caption[]{
Top: power law decay of $\mu_{50}$; the dashed-dotted lines
have slopes $-0.5$ (upper line) and $-0.98$ (lower line).
Bottom: plot of $-1/\mu_{50}$ for Runs~1 (blue),
2 (red), 3 (black), and 4 (orange).
}\label{pcomp_mudiff}\end{figure}

\begin{table}[b]\caption{
All runs have $\mu_{50}=-20$.
For Runs~1--4, the viscosity is the same, so the magnetic Prandtl number is different in all runs. 
}\vspace{12pt}\centerline{\begin{tabular}{ccccccccccc}
Run & $1/\sigma$ &  $\nu$  &  $\mathrm{Pm}$  & $s$ & $N^3$ \\
\hline
0a & $5\times10^{-3}$ & $2\times10^{-3}$ & $0.4$ & $\approx1$ & $448^3$ \\% 448a_mu20_5em3
0b & $5\times10^{-3}$ & $5\times10^{-3}$ & $1$ & $\approx1$ & $448^3$ \\% 448a_mu20_5em3_Pm1
1 & $5\times10^{-4}$ & $2\times10^{-4}$ & $0.4$ & $0.96$ & $1024^3$ \\% 1024a_mu20_5em4
2 & $2\times10^{-4}$ & $2\times10^{-4}$ & $1$ & $0.69$ & $1024^3$ \\% 1024a_mu20
3 & $1\times10^{-4}$ & $2\times10^{-4}$ & $2$ &$0.58$ & $1024^3$ \\% 1024a_mu20_1em4
4 & $5\times10^{-5}$ & $2\times10^{-4}$ & $4$ & $0.51$ & $2048^3$ \\% 2048c_mu20_5em5
\label{Ttimescale4}\end{tabular}}\end{table}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{fig/pcomp_mudiff_jenny}
\end{center}\caption[]{
Same as \Fig{pcomp_mudiff}, but with two more runs with larger $1/\sigma$. 
Also, in addition to $\aver{\mu_5}$ (solid lines), the rms value, $\mu_{5,\mathrm{rms}}$, is shown.
}\label{pcomp_mudiff_jenny}\end{figure}


\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pdecay_exp_vs_sigma}
\end{center}\caption[]{
Dependence of $s$ on $1/\sigma$ for Runs~1--4.
For larger $1/\sigma$, the value of $s$ seems to approach $1$. 
}\label{pdecay_exp_vs_sigma}\end{figure}

%------------------ end of notes 13 ----------------------------

\subsection{Energy fluxes and dissipation}

%\KK{Shall we give a discussion similar to Sec. 4 of notes3.pdf?}

In MHD, there are two important energy reservoirs, kinetic and magnetic
energies.
Both reservoirs can lose energy through dissipation into thermal energy,
and energy is also exchanged between the two reservoirs by work done
against the Lorentz force.

%In \Fig{pdiss_comp} we show the
We have studied the
evolution of various normalized energy fluxes: work done by
the Lorentz force, $W_{\rm L}=\aver{\uu\cdot(\JJ\times\BB)}$ (black),
magnetic forcing by the chiral chemical potential,
$W_\mu=\aver{\mu_5\JJ\cdot\BB)}/\sigma$ (orange),
magnetic energy dissipation, $\epsM=\aver{\JJ^2/\sigma}$ (red), and
kinetic energy dissipation, $\epsK=\aver{2\rho\nu\SSSS^2}$ (blue).
All fluxes are expressed as the ratio of the dissipation wavenumber
to the lowest wavenumber of the domain, i.e., as $k_\nu/k_1$, where
$k_\nu/k_1=(\epsK/\nu^3)^{1/4}$, and likewise for the other quantities.
This allows us to see to what extent the full spectral range is used in
the simulations.
In the present case, the maximum value is of $k_\nu/k_1$ are around 200,
which is about the largest value that can be achieved in simulations
with a Nyquist wavenumber ratio of 256, i.e., with 512 mesh points.

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pboxes1}
\end{center}\caption[]{
Energy flux balance for the model without spin flipping.
}\label{pboxes1}\end{figure}

In \Fig{pboxes1} we show the flows between different energy
reservoirs for the case without spin flipping for early,
intermediate, and late times ($\eta=10$, 100, and 1000).
We normalize these fluxes here by $\mu_{50}^4/\sigma^3$.
Since the system is not in a steady state, energy fluxes
can also be provided by the decay of magnetic and kinetic
energies.
This can be seen in the panels for magnetic energy,
$\EEM=\aver{\BB^2/2}$, where $-\dd\EEM/\dd t=38$
(in units of $\mu_{50}^4/\sigma^3$) at $\eta=10$.
This drives kinetic energy through $W_{\rm L}\approx7.8$.
The kinetic energy gains the amount $7.4$
(so $\EEK=\aver{\rho\uu^2/2}$ increases, here with
$-\dd\EEK/\dd t=-7.4$).
Both energy reservoirs also dissipate (here
$\epsM=1.2$ and $\epsK=0.4$ at $\eta=10$).

Interestingly, at early times, 
$W_\mu=\aver{\mu_5\JJ\cdot\BB)}/\sigma$ is negative.
This is because $\mu_5$ is negative, and $\JJ\cdot\BB$ is
still positive.
This effectively corresponds to a loss of magnetic energy.
This soon changes, because $\JJ\cdot\BB$ becomes negative,
so that then $W_\mu$ becomes positive, and the chiral chemical
potential drives magnetic energy with $W_\mu=54$ at $\eta=100$.
This then also leads to considerable Joule losses ($\epsM=44$)
and large work by the Lorentz force ($W_{\rm L}\approx12$).

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pboxes2}
\end{center}\caption[]{
Same as \Fig{pboxes1}, but for the model with spin flipping.
}\label{pboxes2}\end{figure}

This is very different in the case with spin flipping,
which strongly reduces the value of $\mu_5$, and so
$W_\mu=4$ at $\eta=100$; see \Fig{pboxes2}.
At late times ($\eta=1000$), the losses of magnetic energy are much
smaller in the case with spin flipping ($-\dd\EEM/\dd t=0.1$) than without
($-\dd\EEM/\dd t=0.3$).
We emphasize that the actual energy content is, however, much
larger with spin flipping, just because the losses are small.
This is a direct consequence of the presence of magnetic helicity.

%\section{Effect of chirality-flipping interaction}
%AB: we discussed to change this to a subsection instead
\subsection{Effect of chirality-flipping interaction}

All the simulations presented so far had $\Gamma=0$ and they resulted
in a final state where $\mu_5$ and $\mu_{\rm M}$ vanish at late times.
As discussed in the introduction, spin flipping could prematurely lead
to a vanishing $\mu_5$, which, because of total chirality conservation,
would imply that the decay of $\mu_{\rm M}$ would slow down and level
off at a value away from zero.
To study this quantitatively, we show in \Fig{p_comp2} the evolution
of $\mu_5$, $\mu_{\rm M}$, and $\mu_5+\mu_{\rm M}$ for two models
%with different $\Gamma_{\rm f}$ ($10^{-2}$ and $10^{-4}$), where
%AB: do we really want/need an index on \Gamma? Avoided it now.
with different $\Gamma$ ($10^{-2}$ and $10^{-4}$), where
$\lambda=2\times10^4$, and one more model where $\lambda=10^4$.

%First, we see that for $\Gamma_{\rm f}=10^{-2}$, $\mu_5$
%AB: do we really want/need an index on \Gamma? Avoided it now.
First, we see that for $\Gamma=10^{-2}$, $\mu_5$
begins to decrease rapidly to zero after $\eta_{\rm f}=100$.
This slows down the decay of $\mu_{\rm M}$, which then declines at a
much smaller rate.

Qualitatively similar behaviors are also seen for $\Gamma=10^{-4}$.
In both cases, we see that $\mu_5+\mu_{\rm M}$ evolves away from zero.
When $\lambda_5$ is decreased, the main difference is that the
magnetic field retains slightly larger values for longer.
This is because the magnetic field then couples less strongly to $\mu_5$
and tends to evolve more independently from total chirality conservation.

Next, it is also of interest to study a case where spin flipping
acts only for a certain time interval and is then turned off again
at $\eta=\eta_{\rm off}$.
This case is shown in \Fig{p_comp3},
where we show mainly the evolution after $\eta=\eta_{\rm off}$.
%We see that, when $\Gamma_{\rm f}=0$ after $\eta_{\rm off}$,
%AB: do we really want/need an index on \Gamma? Avoided it now.
We see that, when $\Gamma=0$ after $\eta_{\rm off}$,
the sum $\mu_5+\mu_{\rm M}$ is strictly constant and away from zero, while for
%$\Gamma_{\rm f}\neq0$, the sum continues to decrease slowly.
%AB: do we really want/need an index on \Gamma? Avoided it now.
$\Gamma\neq0$, the sum continues to decrease slowly.
This also means that $\mu_{\rm M}$ stops to decline rapidly at a larger value.
In the end, however, $\mu_{\rm M}$ must still decline to zero.

%\KK{The effect of the flipping rate is studied based on notes3.pdf. }
%------------------ notes 3 ----------------------------

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{p_comp2}
%\includegraphics[width=\columnwidth]{p_comp_1024a_mu10_Gam1em2_eflip1e2_lam2e4}
\end{center}\caption[]{
Evolution of $\mu_{\rm M}$ (red), $\mu_5$ (blue), and their sum (black) for models with
$\eta_{\rm f}=100$, $\sigma^{-1}=2\times10^{-4}$, and
%$\Gamma_{\rm f}=10^{-2}$ and $10^{-4}$ for $\lambda_5=2\times10^4$,
%AB: do we really want/need an index on \Gamma? Avoided it now.
$\Gamma=10^{-2}$ and $10^{-4}$ for $\lambda_5=2\times10^4$,
and the latter case also for $\lambda_5=10^4$.
}\label{p_comp2}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{p_comp2_comp}
%\includegraphics[width=\columnwidth]{p_comp_1024a_mu10_Gam1em2_eflip1e2_lam2e4}
\end{center}\caption[]{
Same as \Fig{p_comp2}, but compensated with $\sqrt{\eta/\sigma}$.
}\label{p_comp2_comp}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{p_comp3}
\end{center}\caption[]{
Effect of resetting spin flipping to zero after some time.
The solid line corresponds to the run with $\lambda_5=10^4$
%and $\Gamma_{\rm f}=10^{-4}$, i.e., the same as the solid
%AB: do we really want/need an index on \Gamma? Avoided it now.
and $\Gamma=10^{-4}$, i.e., the same as the solid
lines in \Fig{p_comp2}.
The dotted lines denote the same run without spin flipping.
}\label{p_comp3}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{p_comp5_full}
\includegraphics[width=\columnwidth]{p_comp5}
\end{center}\caption[]{
%Same as \Fig{p_comp3}, but for $\Gamma_{\rm f}=10^{-2}$ and $k_1=0.5$.
%AB: do we really want/need an index on \Gamma? Avoided it now.
Same as \Fig{p_comp3}, but for $\Gamma=10^{-2}$ and $k_1=0.5$.
}\label{p_comp5}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pscales_128a2_mu10_Gam1em4_eflip1e2_lam2e4_noff_k05}
\includegraphics[width=\columnwidth]{pscales_128a_mu10_Gam1em2_eflip1e2_lam2e4_off1e3_k01}
\end{center}\caption[]{
$\mu_5$ black, $\kf$ red, and $\iota$ blue.
Note the sign reversal of $\tilde{\mu}_5$ at around $\eta=450$.
In both panels, we have $\Gamma\neq0$ for $\eta>\eta_{\rm flip}=10^2$,
but in the lower panel we have $\Gamma=0$ for $\eta>\eta_{\rm off}=10^3$,
Note that then, $\tilde{\mu}_5$ (black line) is clearly different from
zero and positive at late times.
Note that in both cases, $\mu_{\rm M}$ (red line) stays away from zero.
}\label{pscales}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{p_comp6_full}
\includegraphics[width=\columnwidth]{p_comp6}
\end{center}\caption[]{
%Same as \Fig{p_comp3}, but for $\Gamma_{\rm f}=10^{-3}$ and $k_1=0.5$.
%AB: do we really want/need an index on \Gamma? Avoided it now.
Same as \Fig{p_comp3}, but for $\Gamma=10^{-3}$ and $k_1=0.5$.
}\label{p_comp6}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pscales_128a_mu10_Gam1em3_eflip1e2_lam2e4_noff_k05}
\includegraphics[width=\columnwidth]{pscales_128a_mu10_Gam1em3_eflip1e2_lam2e4_off1e3_k05}
\end{center}\caption[]{
$\mu_5$ black, $\kf$ red, and $\iota$ blue
for the model of \Fig{p_comp6}.
}\label{pscales}\end{figure}

%AB: moved from above to down here
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pkmax_comp_1024a_mu10_Gam1em2_eflip4e2_noff}
\end{center}\caption[]{
%Same as \Fig{pkmax_comp_I}, but with $\Gamma_{\rm f}=10^{-2}$ turned on
%AB: do we really want/need an index on \Gamma? Avoided it now.
Same as \Fig{pkmax_comp_I}, but with $\Gamma=10^{-2}$ turned on
after $\eta>400$.
Here, $k_1=0.05$, as opposed to $k_1=0.01$, which is used further below.
}\label{pkmax_comp_1024a_mu10_Gam1em2_eflip4e2_noff}\end{figure}

%AB: NEW
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pkmax_comp_1024a_mu10_k001c}
\end{center}\caption[]{
Same as \Fig{pkmax_comp_I}, but with $\Gamma=0$
\red{and larger domains with $k_1=0.01$ and 0.02.}
}\label{pkmax_comp_1024a_mu10_k001c}\end{figure}

%Fig.20
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pkmax_comp_flip}
\end{center}\caption[]{
%With spin flipping, $\lambda=10^4$ and $\Gamma_{\rm flip}=10^{-3}$
%AB: do we really want/need an index on \Gamma? Avoided it now.
With spin flipping, $\lambda=10^4$ and $\Gamma=10^{-3}$
%(as opposed to $\Gamma_{\rm flip}=10^{-2}$ that is used below.
%AB: do we really want/need an index on \Gamma? Avoided it now.
(as opposed to $\Gamma=10^{-2}$ that is used below.
$k_1=0.05$.
}\label{pkmax_comp_flip}\end{figure}

%Fig.21
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pkmax_comp_1024a_mu0_Gam1em2_lam0}
\end{center}\caption[]{
%With spin flipping, $\lambda=2\times10^4$ and $\Gamma_{\rm flip}=10^{-2}$.
%AB: do we really want/need an index on \Gamma? Avoided it now.
With spin flipping, $\lambda=2\times10^4$ and $\Gamma=10^{-2}$.
$k_1=0.05$.
}\label{pkmax_comp_1024a_mu0_Gam1em2_lam0}\end{figure}

%Fig.22
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pkmax_comp_1024d_mu5_Gam1em2_lam2e4}
\end{center}\caption[]{
%\red{With spin flipping, $\eta_{\rm flip}=100$ and $\Gamma_{\rm flip}=10^{-2}$,}
%AB: do we really want/need an index on \Gamma? Avoided it now.
\red{With spin flipping, $\eta_{\rm flip}=100$ and $\Gamma=10^{-2}$,}
and $\lambda=2\times10^4$. But smaller $k_1=0.01$, allowing more dynamical range toward small $k$.
Thick dashed line with $\mu_5=-5$ initially, and thin solid line with $\mu_5=0$.
This runs has been performed with additional slope-limited diffusion
being invoked.
}\label{pkmax_comp_1024d_mu5_Gam1em2_lam2e4}\end{figure}

The run shown in \Fig{pkmax_comp_1024d_mu5_Gam1em2_lam2e4}
raises the question whether we still get the slower $\eta^{-1/3}$ decay
(or now also the faster $\eta^{-2/3}$ decay) when allowing for $k_1=0.01$.
To address this question, we now perform a run without spin flipping and
otherwise the same parameters.
The result is shown in XX

%Fig.22b
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pkmax_comp_1024d_mu5_Gam0_lam2e4}
\end{center}\caption[]{
The run shown in \Fig{pkmax_comp_1024d_mu5_Gam1em2_lam2e4}
raises the question whether we still get the slower $\eta^{-1/3}$ decay
(or now also the faster $\eta^{-2/3}$ decay) when allowing for $k_1=0.01$.
To address this question, we now perform a run without spin flipping and
otherwise the same parameters.
}\label{pkmax_comp_1024d_mu5_Gam0_lam2e4}\end{figure}

%Fig.23 (but now Fig.24)
\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pkmax_comp_1024d_mu5_Gam1em2_lam2e4_flow}
\end{center}\caption[]{
Comparison with earlier (shorter c) run with initia flow field.
The thick line is a `d' run.
In all cases, we use $k_1=0.01$.
There is a $1/3$ slope, but the range is a bit too short still.
It is therefore not clear (and also unexpected) if the 2/3 slope had
anything to do with the inclusion of a velocity field.
}\label{pkmax_comp_1024d_mu5_Gam1em2_lam2e4_flow}\end{figure}

%1024a_mu0_Gam1em2_lam0    
%1024d_mu0_Gam1em2_lam2e4  mu5=0  
%1024d_mu5_Gam1em2_lam2e4  mu5=-5

%                                                 k1
%aug 18 10:02 1024d_mu0_Gam0_lam2e4              0.01 (just submitted)
%aug 18 09:59 1024a_mu0_Gam1em2_lam0             0.05  2022-08-18 (cvs problem)
%aug 14 07:50 1024d_mu0_Gam1em2_lam2e4           0.01 2/3 2022-08-13 (also with pkmax_comp_1024d_mu5_Gam1em2_lam2e4.pro)
%aug 12 10:21 1024d_mu5_Gam1em2_lam2e4           0.01 2/3 2022-08-12 (Fig.22, with pkmax_comp_1024d_mu5_Gam1em2_lam2e4.pro, t=23387)
%aug 11 15:47 1024c_mu5_Gam1em2_lam2e4           0.01 1/3 2022-08-11 (initial urms, not long enough? t=12412, compared with d run, initial helical flow)
%aug 11 15:46 1024c_mu0_Gam1em2_lam2e4           0.01 2/3 2022-08-10 (Fig.21, with pkmax_comp_1024a_mu0_Gam1em2_lam0.pro)
%aug 11 14:47 1024b_mu0_Gam1em2_lam0             0.01  -  2022-08-10 crashed, because of no nu-slope-limited
%aug  9 19:29 1024c_mu0_Gam1em2_lam0             0.01 2/3 2022-08-09 (thin line in Fig.22, with pkmax_comp_1024d_mu5_Gam1em2_lam2e4.pro and pkmax_comp_1024a_mu0_Gam1em2_lam0.pro)
%aug  8 07:27 1024d_mu10_Gam1em2_eflip4e2_noff   0.05  2022-07-27
%aug  8 07:25 1024a1_mu10                        0.05  2021-11-14
%aug  8 02:36 1024a_mu0_Gam1em2_lam2e4           0.05  2022-08-18 (cvs problem)
%aug  7 13:13 1024a_mu10                         0.05  2021-11-14 Run I
%jul 23 07:51 1024c_mu10_Gam1em2_eflip4e2_noff   0.05 1/3 (Fig.5, sim to Fig.20, with pkmax_comp_1024c_mu10_Gam1em2_eflip4e2_noff.pro)
%jul 23 03:35 1024b_mu10_Gam1em2_eflip4e2_noff
%jul 21 16:29 1024a_mu10_Gam1em2_eflip4e2_noff   0.05 1/3 (Fig.5, with pkmax_comp_1024c_mu10_Gam1em2_eflip4e2_noff.pro)
%jul 21 15:24 1024a_mu10_Gam1em2_eflip1e2_lam2e4
%jul 21 15:24 1024a_mu10_Gam1em3_eflip1e2
%jul 21 15:24 1024a_mu10_Gam1em3_eflip1e2_lam2e4
%jul 21 15:24 1024a_mu10_Gam1em3_eflip1e3
%jul 21 15:24 1024a_mu10_Gam1em2_eflip1e2

%-rw-r--r-- 1 brandenb brandenb 2,2K aug 18 09:49 pkmax_comp_1024d_mu5_Gam1em2_lam2e4.pro
%-rw-r--r-- 1 brandenb brandenb 2,1K aug 12 10:23 pkmax_comp_1024a_mu0_Gam1em2_lam0.pro
%-rw-r--r-- 1 brandenb brandenb 1,9K jul 27 05:06 pkmax_comp_1024c_mu10_Gam1em2_eflip4e2_noff.pro
%-rw-r--r-- 1 brandenb brandenb 1,9K jul 22 05:19 pkmax_comp_1024a_mu10_Gam1em2_eflip4e2_noff.pro
%-rw-r--r-- 1 brandenb brandenb 1,8K jul 22 05:19 pkmax_comp_p_comp6.pro

Spin flipping brings $\mu_5$ close to zero.
This process stops or slows down the decline of magnetic helicity, which
therefore remains positive.
At late times, $\mu_5$, which was originally negative, now becomes
positive and settles at a value of around $\mu_5\approx k_1$.

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{pcomp}
\end{center}\caption[]{
Evolution of $\EEM(\eta)$ (black) and $\HHM(\eta)$ (red)
for a model with (dotted lines) and without (solid lines)
spin flipping invoked.
}\label{pcomp}\end{figure}

The evolution of the magnetic energy densities, $\EEM(\eta)$,
are compared in \Fig{pcomp} for runs with and without spin flipping.
We also compare with the evolution of the magnetic helicity densities
$\HHM(\eta)$ in those two cases.

The sign of the final value of $\mu_5$ is determined by the initial
strength of $\mu_5$ relative to $\mu_{\rm M}$.
In the cases presented above, its sign was positive and thus determined
by the magnetic helicity.
If the initial magnetic field is weaker, the sign of the final value
of $\mu_5$ can stay negative.

Our runs show that spin flipping can lead to a significant increase of the fraction
of the magnetic helicity that can be preserved in spite of the fact that
the system has vanishing total chirality.
This also reduces the total energy density dissipation of the system.
In the absence of spin flipping, both magnetic helicity and chiral chemical
potential would approach zero, so there would be no magnetic helicity available
for successful baryogenesis.

To approach regime~I, we must
increase the value of $v_\lambda/v_\mu=\sigma\lambda^{-1/2}$.
We can either increase $\sigma$ or decrease $\lambda$.
We recall that the value of $\lambda$ was adjusted so as to have $\mu_{\rm
M}=\mu_5\approx k_{\rm Ny}$, and $k_{\rm Ny}$ was limited by numerical
capabilities and the properties of the inflationary initially condition.
Thus, to compensate for the decrease of $\lambda$, we must increase the
magnetic field strength, which, in the present runs, has already an
energy density of about 2\% of the radiation energy density.
Another possibility is to decrease $\sigma$.
This is possible, but requires simply large values of $k_{\rm Ny}$
and therefore also more mesh points.

%------------------ end of notes 3 ----------------------------

%\section{Non-zero total chirality}
%AB: we discussed to change this to a subsection instead
\subsection{Non-zero total chirality}

%\KK{The case for $\mu_5+\mu_M \not = 0$ is studied based on notes.pdf. }

\begin{table}[b]\caption{
Values of $(\mu_{50}+\mu_{\rm M0})/2$ for the combinations
of $\mu_{50}$ and $\mu_{\rm M0}$ considered in this paper.
}\vspace{9pt}\centerline{\begin{tabular}{rrrr}
& \multicolumn{3}{c}{$\mu_{\rm M0}$} \\
\cline{2-4}
$\mu_{50}$ & $-2$ & $-9$ & $-24.5$ \\
\hline
40 & $ 19$ & $ 15$ & $ 8$ \\
20 & $  9$ & $  5$ & $-2$ \\
10 & $  4$ & $ 0.5$& $-7$ \\
\label{Tlist}\end{tabular}}\end{table}

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pall}
\end{center}\caption[]{
$\mu_5$ (blue), $\half\lambda\aver{\AAA\cdot\BB}$ (red), and their sum
(black) for $\mu_5=40$ (top row), $20$ (middle row), and $10$ (bottom
row), for three different magnetic field strengths ($B_0=0.03$, $0.06$,
and $0.1$).
}\label{pall}\end{figure}

In \Fig{pall} we show the evolution of $\aver{\mu_5}$ and $\aver{\mu_{\rm M}}$.
To understand the significance of the nine panels of \Fig{pall}, it is
useful to keep an eye on \Tab{Tlist}, where the nine runs are organized
in the same order as the panels of \Fig{pall}.
In particular, we see that for seven runs, the sign of $\mu_{\rm tot}$
is positive (referred to below as runs of type~I), and for two runs, the
sign is negative (referred below as runs of type~II).

We clearly see that the two runs where $\mu_{\rm tot}>0$ (type~I) behave
qualitatively differently from those where $\mu_{\rm tot}<0$ (type~II).
All simulations of type~I have in common that the magnetic field initially
approaches zero while $\aver{\mu_5}$ approaches $\mu_{\rm tot}$.
This is the stage~I identified in the work of \cite{Hirono+15}.
At some later point, the chiral dynamo instability sets in, and a helical
magnetic field begins to grow until it reaches some maximum value,
while $\aver{\mu_5}$ diminishes to its minimal value $k_1$.
This is the marginal value for the dynamo action.

This situation is different for runs of type~II, where both $\aver{\mu_5}$
and $\aver{\mu_{\rm M}}$ show an initial approach towards the same value,
which turns out to be about $\mu_{\rm tot}/2$.
Only later, they begin to separate again in such a way that $\aver{\mu_5}$
approaches its marginal value, but with the opposite sign, namely $-k_1$.
Correspondingly, $\aver{\mu_{\rm M}}$ then approaches the value $\mu_{\rm tot}+k_1$.
This is summarized in \Tab{Tsum}.

\begin{table}[b]\caption{
Early and late values of $\aver{\mu_{\rm M}}$ and $\aver{\mu_5}$ for runs of
types~I and II.
}\vspace{9pt}\centerline{\begin{tabular}{ccccc}
& Type~I & Type~II & Type~III \\
& $|\muM|>\mu_5|$ & $|\muM|<\mu_5|$ & $|\muM|=\mu_5|$ \\
\hline
early $\aver{\mu_{\rm M}}$ & 0               & $\to\mu_{\rm tot}/2$ \\
early $\aver{\mu_5}$ & $\mu_{\rm tot}$ & $\to\mu_{\rm tot}/2$ \\
late $\aver{\mu_{\rm M}}$ & $\mu_{\rm tot}-k_1$ & $\mu_{\rm tot}+k_1$ \\
late $\aver{\mu_5}$ & $k_1$               & $-k_1$ \\
\label{Tsum}\end{tabular}}\end{table}

For stronger fields ($\mu_{\rm M0}=-24.5$), we see that there is always an
early relaxation phase for $\eta k_1^2 t\la10^{-2}$.
During that time, we see that in most of our cases, $\mu_5$ and $\mu_{\rm M}$
approach the intermediate values $\mu_{\rm tot}$ and $0$, respectively.
The two exceptions are those with $(\mu_{50},\mu_{\rm M0})=(20,-24.5)$
and $(10,-14.5)$, i.e., where $|\mu_{\rm M0}|>|\mu_{50}|$,
which is the regime considered by \cite{Schober+2020}.
In that case, both approach a common line in between the two,
at $\approx(\mu_{50}+\mu_{\rm M0})/2$, which is close to the values $-2$
and $-7$ seen in \Figsp{pspec_all}{f}{i}.

\begin{figure}[t!]\begin{center}
\includegraphics[width=\columnwidth]{pspec_all}
\end{center}\caption[]{
Spectra for the same 9 runs as in \Fig{pall}.
Blue colors indicate the first few times (the times are given in
constant logarithmic time intervals), the next few times are red,
and the last few times are in orange.
}\label{pspec_all}\end{figure}

In \Fig{pspec_all}, we show magnetic energy spectra.
To distinguish early, intermediate, and late times, we use correspondingly
the three colors blue, red, and orange.
Within each group of the same color, we use different line types (solid,
dotted, dashed, and dashed-dotted lines).

We see that the two runs of type~II are the only ones where the
spectral peak has approximately the same height at all times; see
\Figsp{pspec_all}{f}{i}.
In \Figp{pspec_all}{h}, we see that the spectral energy declines
monotonically with time, while for the remaining runs, the spectral
magnetic energy begins to grow again at late times.

In \Figp{pspec_all}{g} we see a typical case where the spectrum
undergoes inverse cascading, and with a drastic decay until the
trend reverses and the spectral peak begins to grow again.
This is similar to what was also shown by \cite{Schober+2020};
see their Figure~1(f).

%AB: in this section, we may still want to add a more systematic
%AB: investigation showing how the net final chirality depends on Gamma_flip.

\section{Application to the early universe}

%\subsection{Mapping QED-like theory to the standard model}
%AB: this title was rather long. How about this?
%\subsection{QED applied to the standard model}
%AB: we decided to change like so:
\subsection{From QED to the standard model}

%\ckm{
%Up to here, we focus on a QED-like theory.
%AB: focus -> have focused
Up to here, we have focused on a QED-like theory.
Now we would like to show its relation to the dynamics in the early Universe.
The standard model of particle physics involves three flavors of leptons $e_R$, $\ell_L$ and quarks $u_R$, $d_R$, $q_L$ plus Higgs $\Phi$, which are in total $16$ species.
On top of this, we have gauge interactions of $\mathrm{U}(1)_Y \times \mathrm{SU}(2)_L \times \mathrm{SU}(3)_C$.
%Apparently it is not clear why we can reduce this complicated system to the chiral MHD based on a QED-like theory.
%AB: changed using: clear -> obvious
It is not obvious why we can reduce this complicated system to the chiral MHD based on a QED-like theory.

What we are interested in here is the slow dynamics of a long wave length compared with interactions among particles.
The key idea for the reduction is to assume the equilibration for fast interactions and to keep only the slow variables.
The hypermagnetic field of $\mathrm{U}(1)_Y$ with a correlation length much larger than the mean free path of the particles 
stands out as a slow variable as the magnetic flux cannot be cut thanks to the absence of a monopole.
This feature does not hold for non-Abelian gauge fields because they are charged under their own gauge group.
We also need the chiral chemical potential since it is related to the magnetic field via the anomaly equation.
Aside from these two fundamental building blocks, we can coarse-grain the microscopic property of each particle as transport coefficients such as the diffusion constant and the electric conductivity, 
and macroscopic quantities like the pressure, the energy density, and the velocity field.
In this way, one may see that the system can be reduced to the chiral MHD as far as the slow and long wave dynamics is concerned.

Still, one might wonder why we can just focus on one particular chiral chemical potential as in Eq.~\eqref{dAdt} 
since we have $15$ species of chiral fermions in the standard model.
%AB: we should perhaps quote the many papers that include many more species
To illustrate this, let us focus on the temperature right above $10^5\, \mathrm{GeV}$, where the electron Yukawa interaction is not efficient compared to the cosmic expansion but other interactions are fast enough.
In this case, the chiral chemical potential for the right-handed electron, $\mu_{e}$, should be counted as a slow variable as it is directly related to the hypermagnetic field via the anomaly equation.
On the other hand, other chiral chemical potentials are subject to fast standard model interactions, which provides non-trivial $11$ constraints among them.
Recalling that the standard model has four conserved charges, the hypercharge and the flavored $B/3 - L_f$ with $f = 1,2,3$,
one may immediately see that the remaining $15$ chemical potentials can be expressed as a function of $\mu_e$ by solving $11 + 4$ constraints.
The $\mu_5$ appearing in Eq.~\eqref{dAdt} is originated from the generalized Ohm's law of
\begin{equation}
  %\bm{J} = \sigma \bm{E} + \frac{2 \alpha_Y}{\pi}\mu_5 \bm{B},
%AB: above, we used the symbol J_em. But am not sure that we need it all the time
%AB: we agreed to put subscript Y here
  \bm{J}_Y = \sigma \bm{E}_Y + \frac{2 \alpha_Y}{\pi}\mu_5 \bm{B}_Y,
\end{equation}
with $\alpha_Y$ being a fine structure constant of $\mathrm{U}(1)_Y$.
For the standard model $\mathrm{U}(1)_Y$, one may express this $\mu_5$ as a summation of chiral chemical potentials for the standard model fermions as
%AB: ref here?
\begin{equation}
  \mu_5 = \sum \epsilon_i y_i^2 \mu_i = \frac{711}{481} \mu_e,
\end{equation}
where $i$ runs through all the standard model fermions, $\epsilon_i = \pm$ for the right- and left-handed fermions respectively, and $y_i$ represents the hypercharge of $i$.
In the second equality, we have inserted the solution of $15$ constraints that we have mentioned.
Up to this order one coefficient of $711/481$, it is shown that one particular chiral chemical potential suffices to describe the system.
%}

\subsection{Initial conditions from axion inflation}

\KK{KK and also KS and KM shall give the discussion after we have fixed the scaling law.}
%See \App{NonZeroNetHelicity} for our results.

See \Fig{pkmax_comp_1024Hf1e11_xi3_k5em3a3} for a run with $\xi=3$,
and \Fig{rspec_select_2048Hf1e11_xi7_k2em2a} for a run with $\xi=7$.

%In \Figss{}{}, we show runs with $\xi=3$, keeping still $\mu_{\rm M}=26$.
%The two runs have $k$ ranges.
%The run with $5\times10^{-2}\leq k\leq51.2$ leads to a slower decay than that
%with $5\times10^{-3}\leq k\leq5.12$, where the growth at $k=k_5$ is not resolved.

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{rspec_select_1024Hf1e11_xi3_k5em2a2}
\end{center}\caption[]{
For Run~X for $\xi=3$.
}\label{rspec_select_1024Hf1e11_xi3_k5em2a2}\end{figure}

\begin{figure}\begin{center}
%JS: I think you select the wrong figure here
%AB: still not right?
%\includegraphics[width=\columnwidth]{pkmax_comp_1024Hf1e11_xi3_k5em3a3}
\includegraphics[width=\columnwidth]{pkmax_comp_1024Hf1e11_xi3_k5em2a2}
\end{center}\caption[]{
For Run~X.
}\label{pkmax_comp_1024Hf1e11_xi3_k5em3a3}\end{figure}

\begin{figure}\begin{center}
%JS: I think you select the wrong figure here
%AB: replaced by xi=7 run
%\includegraphics[width=\columnwidth]{rspec_select_1024Hf1e11_xi3_k5em2a2}
\includegraphics[width=\columnwidth]{rspec_select_2048Hf1e11_xi7_k2em2a}
\end{center}\caption[]{
For Run~Y for $\xi=7$.
$5\times10^{-3}\leq k\leq5.12$.
}\label{rspec_select_1024Hf1e11_xi3_k5em2a2}\end{figure}

\begin{figure}\begin{center}
%\includegraphics[width=\columnwidth]{pkmax_comp_1024Hf1e11_xi3_k5em3a3}
\includegraphics[width=\columnwidth]{pkmax_comp_2048Hf1e11_xi7_k2em2a}
%JS: not needed
%AB: don't know
%}
\end{center}\caption[]{
For Run~Y.
$5\times10^{-3}\leq k\leq5.12$.
}\label{pkmax_comp_2048Hf1e11_xi7_k2em2a}\end{figure}

%-rw-r--r-- 1 brandenb brandenb  23K feb  9 09:56 pkmax_comp_1024Hf1e11_xi3_k5em2a2.eps
%-rw-r--r-- 1 brandenb brandenb  21K feb  9 09:53 pkmax_comp_1024Hf1e11_xi3_k5em3a3.eps
%-rw-r--r-- 1 brandenb brandenb 173K feb  9 09:52 rspec_select_1024Hf1e11_xi3_k5em2a2.eps
%-rw-r--r-- 1 brandenb brandenb 207K feb  9 09:51 rspec_select_1024Hf1e11_xi3_k5em3a3.eps

%\section{Summary}
%AB: I think at the end of the paper, we should not just summarize, and maybe not summarize explicitly at all.
\section{Conclusions}

In this work, we have seen that there are two distinct behaviors of runs,
depending just on the sign of the total chirality for given signs of
fermion chirality and magnetic helicity.
Here, we only considered a case where $\aver{\mu_5}>0$ and $\aver{\mu_{\rm M}}<0$.
The borderline case between the two types corresponds to the case where
$\aver{\mu_5}+\aver{\mu_{\rm M}}=0$.

\vspace{2mm}
{\bf Data availability}---The source code used for the
simulations of this study, the {\sc Pencil Code},
is freely available from Ref.~\cite{JOSS}.
The simulation setups and the corresponding data
are freely available from Ref.~\cite{DATA}.

\medskip
\begin{acknowledgments}
Support through the grant 2019-04234 from the Swedish Research Council
(Vetenskapsr{\aa}det) and ...
are gratefully acknowledged.
We acknowledge the allocation of computing resources provided by the
Swedish National Infrastructure for Computing (SNIC)
at the PDC Center for High Performance Computing Stockholm and Link\"oping.

\end{acknowledgments}

\appendix

\section{Non-zero net helicity}
\label{NonZeroNetHelicity}

\subsection{Numerical experiments}

We fix our units by setting $k_1=\cs=\meanrho=1$.
We perform nine sets of simulations.
In all cases, we start with a uniform $\mu_5=\mu_{50}$ and a fully helical
magnetic field with a spectral peak at $k=k_0=60$; see \cite{Bran+17} for
similar setups of a helical magnetic field with a given initial spectrum.
The magnetic field strength is characterized by a typical field strength
$B_0$, so the magnetic helicity is then of the order of $B_0^2/k_0$, but
there are numerical factors that depend on details of the spectrum.
We perform three groups of runs with different initial $\mu_5=40$, $20$,
and $10$, and each with three different magnetic field strengths.
The corresponding values of $\mu_{\rm M0}$ are computed after the simulation
has run for a few time steps.
This is because in the first moments after the simulation has started,
$\mu_{\rm tot}$ shows slight adjustments due to discretization errors.
Furthermore, the magnetic helicity corresponds only approximately to
$B_0^2/k_0$ due to some numerical factors.

In \Tab{Tlist}, we summarize the initial chiralities, 
$\mu_{50}$ and $\mu_{\rm M0}$ and their average values,
$(\mu_{50}+\mu_{\rm M0})/2=\mu_{\rm tot}/2$.
The significance of these values and their signs
will become evident when discussing the results.

In all cases, we keep $\eta=10^{-4}$ and $\lambda=2\times10^4$.
This means $v_\mu=\eta\mu=0.004$, 0.002, and 0.001, and
$v_\lambda=\mu/\sqrt{\meanrho\lambda}=0.28$, 0.14, and 0.07,
so their ratio is always $v_\lambda/v_\mu=70$, which corresponds
to case~B of \cite{Bran+17}.
Time is presented in diffusion times as $\eta k_1^2 t$ on the abscissa.
In addition, it is sometimes presented logarithmically, so as to be able
to distinguish the early and late phases in the evolution.

\bibliography{ref}{}
%\bibliographystyle{aasjournal}
\bibliographystyle{apsrev4-1}
\end{document}

%AB: to move to elsewhere
When $v_\lambda/v_\mu\ge8$, there is a finite inertial range, whose length
(in wavenumber space) increases with increasing $v_\lambda/v_\mu$.
Large values of this ratio require large numerical resolution and become
eventually prohibitively expensive.
For the early Universe at the electroweak energy scale,
$v_\lambda/v_\mu\approx2400$ can be estimated \cite{Bran+17}.
while the largest value they considered
in their simulation was only about 700 (their Run~A), and even that led
to numerical problems near the end of the runs.
We therefore consider here a compromise with $v_\lambda/v_\mu=70$,
which corresponds to their Run~B.

Finally, we look at the $pq$ diagram in \Fig{ppq_all}.
By plotting spectra at different times versus $\kappa\equiv k\xiM(t)$,
the spectral peak is then automatically always near $\kappa=1$.
If the turbulence is selfsimilar, the spectra can be made to collapse
on top of each other by multiplying them with a monotonous function
of $t$, for example a certain power $\beta$ of $\xiM(t)$, so
$\phi(\kappa)=\xiM^\beta\EM(k\xiM(t),t)$.
The energy is then $\EEM\propto\xiM^{-(\beta+1)}\propto t^{-(\beta+1)q}
\propto t^{-p}$, where $p=(\beta+1)q$.
For a fully helical magnetic field, $\EEM\xiM=\const$, so $p=q$,
or $\beta=0$.
In addition, if the magnetic field is selfsimilar, the spectrum must
be invariant under rescaling $t\to\tau t'$ and $x\to \tau^q x'$.
Since the spectrum has dimensions $L^3 T^{-2}$, it must then obey
$E(\kappa)\to E'(\kappa) \tau^{3q-2} = \tau^{-\beta q}$, so
$\beta=2/q-3$, or $q=2/(3+\beta)$, so $q=2/3$ for $\beta=0$.
Alternative selfsimilar scalings can only be obtained when the relevant
spectra have other dimensions; see \cite{Bran20} for such an example in
the context of the Hall effect, where the relevant spectrum has dimensions
$L^5 T^{-2}$.

The runs of type~II show a striking evolution along the $\beta=0$ line,
as expected for helical runs.
%\red{(Note sure about the evolution toward $p=q=0$. Maybe something is wrong....)}

