# Nonlinear magnetic buoyancy instability of a dynamo magnetic field

Yasin Qazi,1 Anvar Shukurov,1 Devika Tharakkal1 and Frederick A. Gent1,2
1School of Mathematics, Statistics and Physics, Newcastle University,
Newcastle upon Tyne, NE1 7RU, UK
2Astroinformatics, Department of Computer Science, Aalto University, PO Box
15400, FI-00076 Espoo, Finland E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Previous studies of the Parker instability employ an imposed magnetic field to
investigate the instability; however, in galaxies and accretion discs, the
large-scale magnetic field is continuously replenished by the dynamo action
near the disc’s midplane. We consider a mean-field dynamo model for the
generation of a galactic magnetic field, in which an imposed $\alpha$-effect
drives the exponential growth of the magnetic field. The growth is quenched by
the vertical motions produced by magnetic buoyancy and the Parker instability
in a stratified gas layer; these vertical motions assist in the driving of a
mean-field dynamo. The dominant saturation mechanism is the buoyant velocity,
which lifts the field away from the disc plane. This buoyant magnetic field
expands into the corona of the galaxy. Moreover during the non-linear stages
of the dynamo, although the magnetic field is saturated, the dynamo and
buoyancy produce a stable oscillatory state. A simple model 1D model is
compared with direct 3D simulations of the equations of motion and confirms
that magnetic buoyancy alone could quench the dynamo and is responsible for
the observed oscillations. Similar to other buoyancy-driven dynamos the
magnetic fields remain strong out to the boundary of the dynamo region. The
discretization of the density into pockets of gas and the thinning of the disc
are not seen in this model even when this instability is relatively strong
compared to the dynamo, which suggests that the dynamo can support the
galactic magnetic field against the instability.

###### keywords:

instabilities – magnetic fields – MHD – dynamo – galaxies: magnetic fields –
ISM:structure

††pubyear: 2022††pagerange: Nonlinear magnetic buoyancy instability of a
dynamo magnetic field–F

##  1 Introduction

Figure 1: From Model R10h2 (see Table 2 below) the fieldlines for the total
$B$, mean $\langle B\rangle$ and fluctuating $b$ magnetic field are
overplotted onto an isosurface of the gas density near the midplane, contours
of gas density at the rear edges. Mean $\langle B\rangle$ is obtained by
smoothing over a Gaussian kernel of $50$ pc. [FAG: What kernel is used? I
added this here to for visual appeal, where the result is not so analytic.
Check correct model stated.]

Magnetic fields are present in many astrophysical systems, often exhibiting
highly random structure driven by turbulent flows. It is also very common to
find magnetic fields coherently organised on large scales, particularly
aligned to large scale shear flows in rotating objects such as spiral
galaxies, convecting stars or accretion disks. Although the origin of the seed
fields that become these large-scale fields is unclear, the processes which
amplify them are well understood. An electromagnetic dynamo is a process which
converts kinetic energy of a system to magnetic energy . Dynamos play an
essential part of the evolution magnetic fields in spiral galaxies.

[FAG: Replace dummy references in the following paragraph] Among the many
dynamo theories, mean-field dynamo theory (see, e.g., dummy) has been very
successful in modeling galactic magnetic fields (see, e.g., dummy) and offers
a framework in which to study and interpret their results Beck_1996. The mean-
field dynamo produces a magnetic field at scales larger than the motion of the
fluid and demonstrates how differential rotation ($\Omega$-effect) (dummy) and
helical motions ($\alpha$\- effect) (dummy) resulting from the Coriolis force
and ambient gas density stratification lead to a large scale magnetic field.
The $\alpha$-effect occurs due to the twisting of magnetic loops that can form
from magnetic buoyancy. An asymptotic solution to the $\alpha^{2}$ dynamo is
shown in Anvar1983, yielding an approximation for the growth rates of the
magnetic field and also investigating how the dynamos respond to various forms
of $\alpha$.

There are distinct phases to the magnetic dynamo. The first phase is the
kinematic stage in which the magnetic field grows exponentially. During this
_linear_ phase the Lorentz force is too weak to affect the motion of the
fluid, and there may be the development of instabilities such as the magnetic
buoyancy or Parker instability (1966ApJ_145_811P), from which motions may
arise that are able to drive the mean field dynamo. As the magnetic field
grows its Lorentz force influences the velocity field significantly and the
dynamo enters the _nonlinear_ or saturation phase, in which the magnitude of
the magnetic field settles to a statistically steady state.

The Parker instability is a form of the magnetic buoyancy instability, which
is enhanced by cosmic rays. It is a prime example of an effect that magnetic
fields and cosmic rays can have on the structure of the interstellar medium
(ISM). It describes the process in which a density perturbation, arising from
small-scale fluctuations in the magnetic field, increases the local magnetic
pressure and in turn decreases the local plasma density disturbing their local
pressure balance. The weightless magnetic field becomes buoyant and starts to
rise. Due to the density stratification the field expands, increasing its
buoyancy amplifying the instability in the equilibrium between gravitational
forces, gas compression and magnetic tension. The instability is expected to
occur in all magnetized astrophysical disks and is often thought to be
responsible for many astrophysical phenomena, such as disk stability
(LFSR2015), the formation of the thermal gas into discrete pockets
(1966ApJ_145_811P) and the magnetic dynamo (1998Han&Les).

[FAG: First discuss the background results relating to the linear instability
…] The linear state of the magnetic buoyancy instability, with and without
cosmic rays has been investigating both analytically and numerically with
reference to galactic parameters. [FAG: summarize the range of parmaters and
key results in both cases with references … See Devika’s intros…]

The magnetic buoyancy nonlinear instability has been examined by DT2022a, who
detail the evolution and saturation of the Parker instability in a nonrotating
box with imposed density, magnetic field and cosmic ray energy density. Using
snapshots during the nonlinear stage of our model displayed in Fig. 1 it can
be seen that the buoyancy instability injects large-scale vertical structure
into the field from the primarily horizontally aligned dynamo-induced magnetic
field. The nonlinear state of that system shows a complete redistribution of
the density and the magnetic field, which is systematically lost from the
system as the instability approaches steady state. A natural question to ask
then is can a magnetic dynamo sustain the magnetic field against this loss and
how does this affect the evolution of the instability?

Vertical motions produced by magnetic buoyancy and Parker instability in a
stratified gas layer become helical under the action of the Coriolis force in
the same manner as the turbulent motions which drive the $\alpha$-effect, so
it is natural to suggest that the instability can drive a mean-field dynamo
(1992PARKER; 1999MOSS; Thelen_aw_2000). The strength of the mean-field dynamo
is described by two dimensionless parameters $R_{\Omega}$ and $R_{\alpha}$,
which specify the strength of the $\omega$ and $\alpha$-effects respectively.
Typically $|R_{\Omega}|>R_{\alpha}$, this is referred to as an
$\alpha\omega$-dynamo. However mean-field dynamo action is possible under
solid body rotation in this case the magnetic field generated by the
$\alpha$-effect is then acted on again by the $\alpha$-effect, this results in
the $\alpha^{2}$-dynamo (1980mfmd_Krause_Radler; Shukurov2000GlobalMS). Since
we neglect the effect so differential rotation our simulation will be of the
Parker instability on a $\alpha^{2}$ dynamo generated magnetic field. In this
work we shall extend the investigation of the nonlinear state of the magentic
buoyancy instability by generating the required magnetic field via an imposed
$\alpha^{2}$-dynamo instability.

##  2 Methods

[FAG: We should clearly split the paper into methods, results and
dicussion/conclusion.]

###  2.1 Model design

####  2.1.1 Numerical setup

We evolve a section of the galactic disk within a Cartesian box, with $x$, $y$
and $z$ representing the radial, azimuthal and vertical directions,
respectively. The simulation domain extends 4 kpc in each horizontal direction
and 3 kpc vertically, centered around the galactic midplane. The grid
resolution is $256$ $\times$ $256$ $\times$ $384$ mesh points, yielding a grid
spacing along each edge $\approx 15.6\text{pc}$. The size of the domain is
sufficent to be larger than the characteristic vertical and horizontal scales
of the instability, and the resolution sufficient to obtain convergent
solutions to model the instability (see, e.g., DT2022a) and appendix C.

The stratification about the plane of the disk for all models is supported by
the imposition of a gravitational acceleration due to stellar and dark matter
halo following K&G1989MNRAS; Ferriere_1998 expressed as

$g=-\frac{a_{1}z}{\sqrt{z_{1}^{2}+z^{2}}}-a_{2}\frac{z}{z_{2}},$ (1)

with $a_{1}=4.4\times 10^{-14}\,{\rm km}\,{\rm s}^{-2}$, $a_{2}=1.7\times
10^{-14}\,{\rm km}\,{\rm s}^{-2}$, $z_{1}=200\,{\rm pc}$ and $z_{2}=1\,{\rm
kpc}$.

In this study we neglect the effects of rotation and shear and excite an
$\alpha^{2}$-dynamo by imposing an external $\alpha$ through the induction
equation with dependence on $z$ as

$\alpha(z)=\begin{cases}\displaystyle\alpha_{0}\sin\bigg{(}\frac{\pi}{h_{\alpha}}z\bigg{)},&\quad\lvert
z\rvert\leq\dfrac{h_{\alpha}}{2}\\\ \\\
\displaystyle\operatorname{sign}(z)\alpha_{0}\exp\bigg{[}-\bigg{(}\frac{2}{h_{\alpha}}z-\operatorname{sign}(z)\bigg{)}^{2}\bigg{]},&\quad\textrm{otherwise},\end{cases}$
(2)

where $\alpha_{0}$ specifies the turbulent speed and $h_{\alpha}$ a scale
height, to emulate a galactic dynamo, most active about the midplane and
vanishing smoothly towards the boundary. The eighteen models in this study
then span a range of values for these two parameters, as summarised in Table
1. [FAG: here, we should now introduce how the parameter $R_{\alpha}$ is
calculated from this parameter pair. Also clarify how $R_{\alpha}$ and
$R_{\Omega}$ are combined to determine the Dynamo number and how
$R_{\Omega}=0$ affects it.]

The system explored is isothermal so as to focus on the effects of purely
magnetic buoyancy instability and the $\alpha$-effect. The exclusion of
entropic reactions results in a stronger instability because in the adiabatic
case, gas sliding down the field lines heats up adiabaticaly and compresses,
the resulting pressure slows the rate of the gas sliding down the field lines.
For now, we neglect the Parker instability effects from cosmic rays. An
initial condition in isothermal hydrostatic equilibrium is thus achieved by
application of a gas density profile

$\rho(z)=\rho_{0}\;\exp\bigg{[}\frac{a_{1}}{c_{s}^{2}}\bigg{(}z_{1}-\sqrt{z_{1}^{2}+z^{2}}-\frac{a_{2}}{2a_{1}}\frac{z^{2}}{z_{2}}\bigg{)}\bigg{]},$
(3)

given a negligible initial seed magnetic field. The sound speed $c_{s}=15$ km
s-1 is chosen such that $50\%$ of hydrogen (HI) is ionised, giving a mean
molecular weight $\mu=\rho/(m_{g}n)=0.74$, assuming the gas to be comprised of
$70\%$ hydrogen and $30\%$ neutral helium with a negligible amount of heavy
metals. The initial temperature is therefore around $T_{0}=20,000$ K and a
scale height $H=300$ pc. [FAG: this differs from the 200 pc listed in Table
1., I think we can state it here as an approximation either $\sim 200$ or
$\sim 300$ and omit it from the Table, because it is an outcome of the
pressure balance rather than an input parameter.]

The seed magnetic field is applied as Gaussian random noise for the vector
potential component $A_{z}$ with a mean amplitude $\propto\sqrt{\rho}(z)$
maximal at the midplane of $10^{-6}\upmu\text{G}$. The magnetic field in the z
direction $B_{z}$ is therefore initially zero. The weak initial field avoids
any long transient behaviour in the growth of the field, which may occur from
an initially uniform field. [FAG: How is this consistent with observed fields
of a few $\upmu\text{G}$? The resultant saturated field may be.]

####  2.1.2 Basic equations

We solve the system of isothermal non-ideal compressible MHD partial
differential equations using the sixth order in space third order in time
finite difference Pencil Code (brandenburg2002; Pencil-JOSS) as follows:

$\displaystyle\frac{\text{D}\rho}{\text{D}t}$
$\displaystyle=-\rho\nabla\cdot{{\mbox{\boldmath$u$}}}+\nabla\cdot\zeta_{D}\nabla\rho\,,$
(4) $\displaystyle\frac{\text{D}{{\mbox{\boldmath$u$}}}}{\text{D}t}$
$\displaystyle=g\hat{{{\mbox{\boldmath$z$}}}}-\frac{\nabla
p}{\rho}+\frac{{{\mbox{\boldmath$j$}}}\times{{\mbox{\boldmath$B$}}}}{4\pi\rho}+\frac{1}{\rho}\nabla\cdot\rho\nu\mathbf{\tau}+\nabla\left(\zeta_{\nu}\nabla\cdot{{\mbox{\boldmath$u$}}}\right)$
$\displaystyle+\nabla\cdot\left(2\rho\nu_{6}{\mathbf{\tau}}^{(6)}\right)-{{\mbox{\boldmath$u$}}}\nabla\cdot\left(\zeta_{D}\nabla\rho\right),$
(5) $\displaystyle\frac{\partial{{\mbox{\boldmath$A$}}}}{\partial t}$
$\displaystyle=\alpha{{\mbox{\boldmath$B$}}}+{{\mbox{\boldmath$u$}}}\times{{\mbox{\boldmath$B$}}}+\eta\nabla\times{{\mbox{\boldmath$B$}}}+\eta_{6}\nabla^{6}{{\mbox{\boldmath$A$}}},$
(6)

where the material derivative is

$\frac{\text{D}}{\text{D}t}=\frac{\partial}{\partial
t}+{{\mbox{\boldmath$u$}}}\cdot\nabla,$ (7)

$\rho$ and $u$ are the gas density and velocity, $A$ is the magnetic vector
potential and the magnetic field
${{\mbox{\boldmath$B$}}}=\nabla\times{{\mbox{\boldmath$A$}}}$, such that
$\nabla\cdot{{\mbox{\boldmath$B$}}}=0$ by construction.
${{\mbox{\boldmath$j$}}}=\nabla\times c{{\mbox{\boldmath$B$}}}(4\pi)^{-1}$ is
the electric current density, with $c$ the speed of light in a vacuum. The gas
pressure $p=\rho c_{s}^{2}$, with $c_{s}$ the speed of sound.

The traceless rate of strain tensor $\tau$ is

$\centering\tau_{ij}=\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial
u_{j}}{\partial x_{i}}-\frac{2}{3}\delta_{ij}\frac{\partial u_{k}}{\partial
k}\@add@centering$ (8)

with an additional diffusion $\zeta_{\nu}$ proportional to the flow
convergence to resolve shock discontinuities. Artificial diffusion $\zeta_{D}$
is applied to the continuity equation (4) proportional to convergence with a
correction applying to conserve momentum in (5) (following GMKSH20). Hyper
diffusion to resolve grid-scale numerical instabilities proportional to sixth
order derivatives applies with terms including $\nu_{6}$ and $\eta_{6}$
(ABGS02; GMKS21). The key parameters used are shown in table 1

Table 1: Parameter values used in the simulations Quantity | Symbol | Value
---|---|---
Grid resolution | $\delta{{\mbox{\boldmath$x$}}}$ | 0.0156 kpc
Sound speed | $c_{s}$ | 15 km s-1
Kinematic viscosity [FAG: Pm$<1$] | $\nu$ | 0.008 kpc2 Gyr-1
Magnetic diffusivity [FAG: Pm$<1$] | $\eta$ | 0.03 kpc2 Gyr-1
Diffusion time TBC | $h_{\alpha}^{2}/\eta$ | $33.3\,{\rm Gyr}(h_{\alpha}/1\,{\rm kpc})^{2}$
Initial column density [FAG: ?] | $\rho(t=0)$ | 3.6e-18 g cm-3
Dynamo scale height | $h_{\alpha}$ | 0.2, 0.3, 0.6 kpc
Imposed turbulent velocity | $\alpha_{0}$ | 0.15 – 3.0 km s-1
Shock diffusion | $\zeta_{\nu},\zeta_{D}$ | 12 kpc2 $|\nabla\cdot{{\mbox{\boldmath$u$}}}_{-\\{\rm ve}|$
Hyper diffusion | $\nu_{6},\eta_{6}$ |  $10^{-10}$ kpc6 Gyr-1
[FAG: 3 values for $h_{\alpha}=0.2,0.3,0,6$] | $h_{\alpha}^{2}/\eta$ | $33.3\,{\rm Gyr}(h_{\alpha}/1\,{\rm kpc})^{2}$

####  2.1.3 Boundary Conditions

Boundary conditions in both horizontal directions are periodic with large
scale rotation and shear neglected in these models. The vertical extent of the
computational domain is sufficiently larger than the density scale height, so
we do not expect the results to be sensitive to the boundary conditions at
$z=\pm 1.5\,{\rm kpc}$. To exclude artificial magnetic energy leaking into
domain from the surfaces we apply a vertical field condition, with

$\frac{\partial A_{x}}{\partial z}=\frac{\partial A_{y}}{\partial
z}=A_{z}=0\quad\rm{at}\>|z|=1.5\,{\rm kpc},$ (9)

such that on the boundaries $B_{x}=B_{y}=0$ and $\partial B_{z}/\partial z=0$.

For the velocity field we apply a stress free condition on the surfaces with

$\frac{\partial u_{x}}{\partial z}=\frac{\partial u_{y}}{\partial
z}=0\quad\rm{at}\>|z|=1.5\,{\rm kpc}.$ (10)

To permit physical flows across these surfaces, without exciting numerical
instabilities, the boundary condition for $u_{z}$ follows Gent_SN_ISM_1,
damping locally the direction of flow. The magnitude of the flow outwardly
normal to the surface is constrained as $\leq|u_{z}|$ for each location on the
boundary.

The density gradient is conserved on the surfaces for $\ln(\rho)$ as

$\frac{\partial\ln(\rho)}{\partial z}=\pm\frac{1}{0.9}\,{\rm
kpc}^{-1}\quad\text{at}\>|z|=1.5\,{\rm kpc}\,,$ (11)

which assumes that the average density scale height of 900 pc (Lockman84)
applies to the local density on the boundary, which is numerically stable
while retaining a physically motivated density gradient near the surfaces, and
sufficiently far from the thin disk scale heights so as not to significantly
influence the results.

###  2.2 Mean-field methods

[FAG: add description of 1D model solution and discuss how it fits in results?
Explain the averaging conventions used to identify the mean field - horizontal
averages used in the analysis, kernel smoothing used to in visualisation, and
arbitrary used in 1D model. I haven’t yet commented on this section as I do
not yet appreciate how it fits the narrative.] Equation 6 includes a
generation term $\alpha(z)$ which is the mean helicity and controls the
strength of the dynamo. Nondimensionalization of the induction equation
reveals $\alpha(z)$ has units of velocity. We define
$B=B_{0}B^{\prime}$,$\alpha=\alpha_{0}\alpha^{\prime}$, $t=\tau t^{\prime}$,
$(x,y,z)=L_{0}(x^{\prime},y^{\prime},z^{\prime})$ and $u=L_{0}u^{\prime}/\tau$
then,

$\displaystyle\dfrac{B_{0}}{\tau}\frac{\partial B^{\prime}}{\partial t}$
$\displaystyle=\dfrac{1}{L_{0}}\nabla^{\prime}\times\left(\dfrac{L_{0}}{\tau}u^{\prime}\times
B_{0}B^{\prime}\right)\dfrac{\nabla^{\prime}}{L_{0}}\times(\alpha_{0}B_{0}\alpha^{\prime}B^{\prime})+\eta\dfrac{B_{0}}{L_{0}^{2}}\nabla^{2}B^{\prime}$
(12) $\displaystyle\frac{\partial B^{\prime}}{\partial t}$
$\displaystyle=\nabla^{\prime}\times(u^{\prime}\times
B^{\prime})+\dfrac{\alpha_{0}\tau}{L_{0}}\nabla^{\prime}\times(\alpha^{\prime}B^{\prime})$
(13)
$\displaystyle+\dfrac{\alpha_{0}\tau}{L_{0}}\nabla^{\prime}\times(\alpha^{\prime}B^{\prime})+\dfrac{\eta\tau}{L_{0}}\nabla^{\prime
2}B^{\prime}$ (14)

To have this be a dimensionless equation $\alpha_{0}=L_{0}/\tau$ The form of
$\alpha$ is related to the mean helicity $\mathcal{H}$ of the random flow and
is given by the correlation function of velocity and vorticity
$\alpha=-\frac{1}{3}\tau\langle{{\mbox{\boldmath$u$}}}\cdot(\nabla\times{{\mbox{\boldmath$u$}}})\rangle$
with correlation time $\tau$. The meaning of the generation term is explained
in (1980mfmd_Krause_Radler) namely that since the electric field and as a
result the mean electric current are parallel to the magnetic field i.e
${{\mbox{\boldmath$j$}}}=-\sigma\alpha{{\mbox{\boldmath$B$}}}+...$ where
$\sigma$ is the conductivity of the fluid. This current means that screw-like
vortices with a non-zero mean helicity cause the generation of loops in the
magnetic field which are perpendicular to the original magnetic field these
large magnetic loop are clear to see in figure 1. Chapter 7 of
(shukurov_subramanian_2021) shows that interstellar turbulence is helical and
that the mean helicity is an odd function in z this is because a volume
element of the interstellar gas, moving away from the galactic mid-plane $z=0$
s, expands due to the ambient stratification and because of the conservation
of the angular momentum, the expanding gas rotates oppositely to the galactic
rotation. This systematic rotation coordinated with the cloud motion leads to
a non-vanishing correlation $\langle u\cdot\nabla\times u\rangle$, negative at
$z>0$ and positive at $z<0$. Meaning rotation and stratification break the
mirror symmetry of the interstellar turbulence to make it helical.

##  3 Results

[FAG: In the results section let us not include discussion, but stick to
outlining the data and statistical results.]

[FAG: I propose we order the results as follows: show the dynamo solutions
current Fig 8 and Table 2., then Fig 2. Next talk about the stages of the
buoyancy instability, TBC,. ]

###  3.1 The linear phase of the dynamo

####  3.1.1 Signatures of the dynamo: exponential growth

From the same initial condition in hydrostatic equilibrium simulations are
evolved for each combination of $h_{\alpha}$ and $\alpha_{0}$. In the absence
of fluid motions to support generation of the magnetic field it initially
decays, driving vertical flows as the pressure balance adjusts while the
system finds an eigenstate. When the magnetic field is weak and the
electromotive force (EMF), ${\cal
E}={{\mbox{\boldmath$u$}}}\times{{\mbox{\boldmath$B$}}}$, [FAG: use vec not
mathbf] is negligible, the induction equation 6 is linear in $B$ and as long
as the velocity field is steady the magnetic field grows exponentially in time
${{\mbox{\boldmath$B$}}}\approx{{\mbox{\boldmath$B$}}}_{0}\exp(\gamma_{D}t)$,
where the $D$ indicates action of dynamo. There are infinitely many of these
eigenmodes ${{\mbox{\boldmath$B$}}}_{0}$ but the growth rate of the dynamo is
determined by the eigenmode with the largest growth rate.

For larger values of $R_{\alpha}$ the period of this linear phase is short.
[FAG: where are $\tau_{\eta}$ and $\tau_{D}$ defined?] However, for smaller
values of $R_{\alpha}$ when $\tau_{\eta}\sim\tau_{D}$ this period can last
longer. In particular, increasing $h_{\alpha}$ even by small amount can
drastically lengthen this stage. Fig. 2 shows the time series for three
different values of $h_{\alpha}$, we see here that increasing $h_{\alpha}$
even by a small amount causes the time frame of our system to almost double.
Increasing $\alpha_{0}$ speeds up the linear phase. During the linear phase
the magnetic field is contained predominantly in the region
$-h_{\alpha}<z<h_{\alpha}$ where the $\alpha$ effect is strongest and the
kinetic energy density very rapidly reaches a statistically steady state
common to all models, which corresponds to a volume averaged r.m.s. velocity
at around [FAG: confirm from $E_{\rm kni}\simeq 1e-4$ in Fig. 2] $10^{-1}$ km
s-1.

Figure 2:  Evolution of the volume-averaged kinetic and magnetic energy
density for various models spanning a range of values for $h_{\alpha}$ and
$R_{\alpha}$ as identified in the captions. [FAG: The $y$-axis labels are not
log, (the scaling is log, but the quantiites are linear. General $E$ can be
used in the axis label or $e$, since they are listed in the legend. $e$ is
energy density vs E total energy? The units are energy density, so should
include $cm^{-3}$. I am also highly skeptical of the magnetudes, which look
like code units. These should be multiplied by par.unit_energy_density and
will saturate around $1e-13$ to $1e-12$. It would also be helpful to include
the Model labels in each plot. General comment for all plots: fonts should be
at least as large as the text in the caption, and fontsize=float is usually an
option in all labels. ] Table 2:  Growth rates for the r.m.s magnetic field in
the linear phase of the dynamo ($\gamma_{\rm{D}}$), the growth rate of the
Parker instability ($\gamma_{\rm{P}}$) [FAG: Given we exclude CRs, maybe use
$\gamma_{B}$ for buoyancy. Also convention is to use italic subscripts unless
multiple letters or proper nowns] obtained from measuring the enhanced growth
rate in the mean field and the enhanced growth rate measured in the kinetic
energy ($\gamma_{\rm{k}}$) and the period of the non-linear oscillations for
simulations which reached that stage. These quantities rates are given for
various values of the dimensionless dynamo number Rα and the scale height of
the dynamo $h_{\alpha}$. The diffusivity $\eta=0.03$ kpc2/Gyr. For values of
R${}_{\alpha}<5$ we either see no growth and the field decays to values within
the numerical precision, or the field grows slowly on long timescales. The
value of $h_{\alpha}=0.2$ kpc corresponds to the scale height of the density.

Model Rα $h_{\alpha}$ $\alpha_{0}$ $\gamma_{\rm{D}}$ $\gamma_{\rm{P}}$
$\gamma_{\rm{k}}$ [kpc] [km s-1] [Gyr-1] [Gyr-1] [Gyr-1] R3h2 3 0.2 0.45
$-$3.8 - - R3h3 0.3 0.3 $-$2.5 - - R3h6 0.6 0.15 - - - R5h2 5 0.2 0.75 1.6 1.1
2.1 R5h3 0.3 0.5 0.4 - 0.6 R5h6 0.6 0.25 0.1 - - R7h2 7 0.2 1.05 4.4 7.9 9.8
R7h3 0.3 0.7 1.9 2.2 4.2 R7h6 0.6 0.35 0.4 0.5 0.7 R10h2 10 0.2 1.5 12.4 29.1
28.3 R10h3 0.3 1.0 5.5 11.3 10.6 R10h6 0.6 0.5 1.6 1.8 2.6 R15h2 15 0.2 2.25
32.2 86.9 129.6 R15h3 0.3 1.5 14.9 36.6 38.5 R15h6 0.6 0.75 3.7 6.36 8.4 R20h2
20 0.2 3.0 62.7 171.4 226.9 R20h3 0.3 2.0 28.5 74.2 82.7 R20h6 0.6 1.0 7.1
14.7 17.1

Figure 3: Dimensionless growth rates
$\displaystyle\hat{\gamma}=\frac{h_{\alpha}^{2}}{\eta}\gamma=\frac{R_{\alpha}^{2}}{4}-\frac{R_{\alpha}}{2}$
for different values of $h_{\alpha}$. The solid red line shows the analytic
result obtained in (Anvar1983). The other lines are the growth rates obtained
from different simulations with varying $h_{\alpha}$ and $R_{\alpha}$. [FAG:
Please also plot $R_{\alpha}=3$, which are negative]

Each model in Fig. 2 has a clearly exponential growth rate $\gamma_{D}$ for
the magnetic field. The growth rates are listed in Table 2 for all models
where such a linear phase, or decay, exists. These can be fitted to the
asymptotic approximation of maximal growth rate in (Anvar1983)

$\displaystyle\gamma_{\rm
max}=\frac{\eta}{h_{\alpha^{2}}}\left[\frac{R_{\alpha}^{2}}{4}-\frac{R_{\alpha}}{2}+O(1)\right].$
(15)

Normalising equation 15 by the diffusion time
$\hat{t}=t/(h_{\alpha}^{2}/\eta)$ yields dimensionless growth rates, which are
displayed against $R_{\alpha}$ in Fig. 3. There is good agreement between the
analytic and numerical results for all simulations. suggest that this agrees
with We find the critical $R_{\alpha,{\rm crit}}\approx 5$ above which dynamo
growth is excited. This is consistent with (Anvar1983) who conclude that [FAG:
I do not understand this sentence, especially the concept of the extremum at
an isolated point] for helicity of $\alpha(z)$ having an extremum at an
isolated point that the critical dynamo number will lie in the range
$4<R_{\alpha}<10$ and further supported by Section 4 of (Anvar_1987) where
similar critical dynamo numbers are found.

Equation 15 is an asymptotic approximation so at small values of $R_{\alpha}$
the approximation should break down however the analytic and numerical results
still have good agreement even at these small dynamo numbers. In figure 3 you
can see that even for this small sample of $R_{\alpha}$ values that the
analytic growth rates and the numerical growth rates converge.

Figure 4: 2D power spectra from Model XXXX for the $z$ component of the
magnetic field taken 385 pc above the midplane during different epochs. [FAG:
Could your plot these vertically with sharex=True, to fit column width?]

Although the exponential growth of the magnetic field is strong evidence of
dynamo action there are other phenomena which are evidence of dynamo action.

####  3.1.2 Signatures of the dynamo: spectra

Figure 4 shows the power spectra from Model XXXX for the $z$ components of the
magnetic field at $z=385$ pc at times indicated. When $\langle V\rangle=0$
[FAG: What is $V$?] there exists simple solutions to the mean field induction
equation if we assume further that $\alpha$ and $\eta$ are constant then the
solution is a force free field with
$\nabla\times\langle{{\mbox{\boldmath$B$}}}\rangle=\pm
k\langle{{\mbox{\boldmath$B$}}}\rangle$ with $k$ a constant. This reduces the
induction equation to

$\dfrac{\partial\langle{{\mbox{\boldmath$B$}}}\rangle}{\partial t}=(\pm\alpha
k-\beta k^{2})\langle{{\mbox{\boldmath$B$}}}\rangle$ (16)

which has exponentially growing solutions, with a real eigenvalue

$\gamma=\pm\alpha k-\beta k^{2}.$ (17)

[FAG: We keep referring to the reference state and have not at any point
indicated what is the reference state. The results are relevant to all models
and not to a particular case, so let us drop the reference state and talk in
general or mention the actual model name.] The value of $k$ which maximises
this function is $k_{\rm max}=\dfrac{\lvert\alpha\rvert}{2\beta}$ this
indicates that as in Model XXXX a strong dynamo will create a small scale
field.

The spectra from 0.5 Gyr during the linear phase of the dynamo the power
spectra in Figure 4 is being generated in the dynamo generation region $\lvert
z\rvert<h_{\alpha}$ on small scales of $k_{y}\approx k_{x}\approx 15$. A
systematic transfer of the energy is seen in the subsequent moves toward
larger scales at 1.5 and 2.5 Gyr. This shift to larger scales is sometimes
referred to as the inverse cascade (Brandenburg_2001), which describes the
transfer of energy to systematically increasing scales. In reality this
transfer is not from systematically increasing scales but instead directly
from small scales to the largest attainable in the system. This type of non-
local spectral energy transfer is associated with the $\alpha$-effect and
evidence of dynamo action.

####  3.1.3 Signatures of the dynamo: magnetic loops

Fig. 1 shows the stream lines of the total, mean and fluctuating field for
Model XXXX just after (1.5 Gyr) and later (2.5 Gyr) into the saturation phase
of the dynamo. Magnetic loops form with a scale of around $1$ kpc.

The appearance of these loops is due to the combined action of magnetic
buoyancy and the $\alpha$-effect. Buoyant fluid parcels ascending (or
descending) from the equatorial plane, in the limit of high magnetic Reynolds
numbers $R_{m}$, when the magnetic field is frozen into the flow the buoyant
fluid parcels distort the horizontal magnetic field into the shape of capital
Greek Omega forming these large scale loops. The local magnetic loops are
twisted together with the gas. Each of these loops has an associated electric
current orthogonal to the plane of the loop. A large scale horizontal current
forms from the horizontal field by small helical motions by the
$\alpha$-effect, which generates an orthogonal magnetic field. These buoyant
loops are further evidence supporting the presence of the dynamo process.

Fig. 11 shows radial-altitude snapshots of each component of the evolving
magnetic field for Model XXXX. These also exhibit the large-scale helical
structures about the scale of $1-2$kpc around saturation of the magnetic field
that are indicative of the $\alpha$-effect.

###  3.2 The linear phase of the buoyancy instability

####  3.2.1 Instability assisted dynamo

In Fig. 2 models with sufficiently fast dynamo approach equipartition with the
kinetic energy density, and hence a horizontally aligned strong magnetic
field, which can become sensitve to the magnetic buoyancy instability. In
these models a transitional phase emerges in which the magnetic buoyancy
instability drives a new exponential growth $\gamma_{B}$ in the magnetic and
$\gamma_{K}$ kinetic energies, which are also listed in Table 2, where they
exist. This associated growth of the magnetic field due to the instability is
well documented (1992PARKER; 1999MOSS; Thelen_aw_2000). For the larger dynamo
numbers of $R_{\alpha}=15,20$, $\gamma_{k}$ is almost double $\gamma_{P}$.

The dynamo velocity field can be driven by direct forcing or by an
instability, either an internal fluid instability such as convection or, more
importantly for this work, a magnetic instability such as magnetic Rayleigh-
Taylor or Parker instability. The distinct secondary growth in the velocity
field after the kinetic energy and the magnetic energy reach equipartition is
likely due to the redistribution of the gas density, which causes there to be
large inflows in the midplane where the most mass displacement takes place.

Vertical motions produced by magnetic buoyancy or Parker instability in a
stratified gas layer become helical under the action of the Coriolis force.
The instability could therefore assist in driving a mean-field dynamo
(1999MOSS; 1992PARKER) in the same manner as turbulent motions. Fig. 5 plots
the strength for the mean (solid) and fluctuating (dashed) magnetic field
within a scale height of the midplane $z<\lvert h_{\alpha}\rvert$ (blue) and
$z>\lvert h_{\alpha}\rvert$ (red). The mean field shows a clear secondary
growth rate, indicative of the buoyancy instability assisting the mean-field
dynamo, whereas the fluctuating field does not.

Applying Gaussian smoothing with scaling length tending to zero or horizontal
averaging of a mean field with nontrivial vertical structure, as is evident
here, then some of the mean-field energy will be recorded as turbulent field.
Where the turbulent field is weak relative to the mean field at $z>\lvert
h_{\alpha}\rvert$ the dashed red line shows a weak signature of secondary
growth from the mean field, which is absent for the dashed blue line at
$z<\lvert h_{\alpha}\rvert$ for the relatively stronger turbulent field.

Figure 5:  The evolving magnitude of the mean $\langle B\rangle$ and
fluctuating $b$ magnetic field averaged over $z<\lvert h_{\alpha}\rvert$
(blue) and $z>\lvert h_{\alpha}\rvert$ (red). [FAG: Have you multiplied the
data by par.unit_magnetic or used code units?]

####  3.2.2 Anisotropic scale of linear instability

For Model XXXX the linear stage of the instability is apparent in Fig. 5 at
around $t=1.5\,{\rm Gyr}$. In the magnetic power spectra at $t=1.5\,{\rm Gyr}$
in Fig. 4 there is a shift to larger structures in the radial direction,
perpendicular to the mean-field with a peak at $k_{x}\simeq 5$, while a
secondary peak appears isotropic in $k_{x}\simeq k_{y}\simeq 10$. This is a
clear shift of power from the smaller scales associated with the strong dynamo
toward the characteristic length scales associated with the buoyancy
instability of around 1-2 kpc and the horizontal scale of these unstable modes
are several times that of the scale height of the disc.

####  3.2.3 Parity of the linear instability

The magnetic field parallel to the disc plane is buoyant, because regions with
stronger field have a reduced gas density to maintain the total pressure
balance along the disc plane. The density displacement due to the magnetic
pressure alone within a flux tube containing a magnetic field $B$ is
$\Delta\rho/\rho=B^{2}/(8\pi\rho)=V_{A}^{2}/(2c_{0}^{2})$ where $V_{A}$ is the
Alfvén speed within the tube and $c_{0}$ is the sound speed.

Magnetic buoyancy leads to the magnetic Rayleigh-Taylor and Parker
instabilities. The Parker instability and the $\alpha$-effect both arise from
buoyancy effects, and so trying to disentangle the effects of each is
difficult. However, there are multiple distinct signatures of the instability.
The first and most clear example is the shift of the dominant modes in the
power spectrum from smaller wavelengths to the larger characteristic
wavelength of the instability, which is evident in the latest power spectra
displayed in Fig. 4.

Linearized analyses of the instability (Giz1993; 1998Kim_Hong) propose that
instabilities will possess a symmetry with respect to reflection in the mid-
plane. The instability is called symmetric if $u_{x},u_{y},B_{x},B_{y}$ are
symmetric and $u_{z},B_{z}$ are anti-symmetric with respect to the mid-plane,
and anti-symmetric when $u_{x},u_{y},B_{x},B_{y}$ are anti-symmetric and
$u_{z},B_{z}$ are symmetric with respect to the mid-plane. The method for
calculating the parities is outlined in appendix D. We find that during the
linear phase of the dynamo there is no preference for either odd or even modes
since the dynamo produces a random magnetic field, when the instability sets
in there is only a slight preference for even parity modes. However when
investigating the parity of the mean and fluctuating fields we find that the
mean field has a clear preference for even parity modes whereas the
fluctuating field has no preference. This is to be expected since the
instability mainly affects the mean field.

###  3.3 Nonlinear buoyancy instability

Saturation effects, such as the associated removal of the magnetic field by
vertical velocity will eventually halt the growth of the magnetic field, but
as mentioned earlier the instability is also viable mechanism for the
saturation of the magnetic field. Figure 2 shows the steady state of different
simulations, when both kinetic and magnetic energy have saturated. [FAG: We
should clarify what is of same order. $E_{B}$ vs $E_{K}$ within a run? $E_{B}$
vs $E_{B}$ for runs with common $R_{\alpha}$?] The steady state energy
densities for the magnetic and velocity fields are of the same order of
magnitude across the same value of $R_{\alpha}$. Increasing $h_{\alpha}$
causes a decrease in the saturation strength. Dynamo action can be inhibited
by magnetic tension that resists the twisting of the $\Omega$-shaped loops by
the Coriolis force to reduce the effective magnitude of $\alpha$. If the
Lorentz and Coriolis forces are balanced (when the bouyant velocity
$u_{B}\approx V_{A}$ ) then the instability would work as a saturation
mechanism for the magnetic field. Although the vertical velocity alone can
saturate the magnetic field through advection, the saturation itself suggests
the presence of the instability.

####  3.3.1 Density variation in the buoyancy instability

Figure 6:  Horizontally averaged density from snaphots at times indicated in
the legends, for models with $R_{\alpha}=10$ and $h_{\alpha}=0.2$ kpc (left)
,$h_{\alpha}=0.3$ kpc (middle) and $h_{\alpha}=0.6$ kpc (right). [FAG: Could
your plot these vertically with sharex=True, to fit column width?]

Fig. 6 shows the difference in the horizontal averages of the vertical density
distribution from the initial condition for simulation with $R_{\alpha}=10$
and $h_{\alpha}=0.2$, 0.3 and 0.6 kpc. Using an imposed magnetic field
(DT2022a) find the instability to cause the density scale height to reduce.
With a dynamo generated magnetic field the increased magnetic pressure near
the midplane resists the effect of the buoyancy instability to increase the
scale height of the gas. The profile of the vertical gas distribution can be
approximated by

$\rho=\rho_{1}\exp\left(-\left(\frac{z}{2h_{1}}\right)^{2}\right)+\rho_{2}\exp\left({-\frac{z}{h_{2}}}\right).$
(18)

A Kolomagorov-Smirnoff test to this fit for the gas density at various stages
of Model XXXX indicates the gas scale height decreases by 15% (between when
and when?) and that the width [FAG: full width half maximum?] of the disk
increases by 50%. [FAG: These are contradictory, if the scale height reduces,
how does the disk get wider?]

####  3.3.2 Nonlinear oscillations of the instability

Figure 7:  Time-altitude diagrams for Model XXXX of
$\langle\hat{B}_{x}\rangle_{xy}$ and $\langle\hat{B}_{y}\rangle_{xy}$
normalized to their time-dependent maximal magnitude.[FAG: please use
sharex=True and increase the fontsize a little]

Figure 7 shows the horizontal averages for the $B_{x}$ and $B_{y}$ from Model
XXXX. To reveal the topological variation at early times the quantities are
normalised by the evolving maximal magntitude. Large-scale structure in the
initial condition for the field dissipates within 300 Myr and the linear phase
of the dynamo with large $R_{\alpha}$ is characterised by temporal and spacial
small-scale, chaotic structure to beyond 1 Gyr, when an ordered horizontal
field emerges. Throughout the linear dynamo the ordered field is concentrated
within $|z|\lesssim 0.2$ kpc and is consistently dipolar and nonoscillatory.

With the onset of the buoyancy instability the scale height of the field
extends to about 1 kpc, the ordered field adopts partial quadripolar
charcteristics and exhibits periodic reversals. While the scale heights,
timing and period of the reverals varies as function of $R_{\alpha}$ and
$h_{\alpha}$, this is a general property of the models with sufficiently
strong mean fields.

The increasing scale and period in space of the instability reversals, in
contrast to the earlier $\alpha^{2}$-dynamo generated field, is also evident
in the later panels of Fig. 11.

###  3.4 The $\alpha$ effect

There is a phase shift in Model XXXX of around quarter of a period between the
$x$ and $y$ components of the magnetic field. The period of oscillation in
various models is dependant on $h_{\alpha}$, increasing the dynamo scale
height causes the the period of oscillation to increase. An important aspect
of the nonlinear $\alpha$-effect is the evolution of the magnetic helicity,
which can be apporximated by the averaged turbulent quantities, denoted by ′
as [FAG: Yasin, equations are not paragraphs and can be subclauses of
sentences, so should not have a line space between them and the text, unless
they start or end a paragraph.]

$\alpha=\alpha_{K}+\alpha_{M}=-\frac{\tau}{3}\left[\langle{{\mbox{\boldmath$u$}}}^{\prime}\cdot(\nabla\times{{\mbox{\boldmath$u$}}}^{\prime})\rangle-\frac{1}{4\pi\rho}\langle{{\mbox{\boldmath$B$}}}^{\prime}\cdot(\nabla\times{{\mbox{\boldmath$B$}}}^{\prime})\rangle\right]$
(19)

If there are no helicity fluxes the magnetic helicity would quench the dynamo
due to the invariance of magnetic helicity in a closed system. Therefore
$\alpha$ requires the contribution from the turbulent current helicity
$\alpha_{M}$, driving motions through the Lorentz forecs, as well as the
turbulent kinetic helicity $\alpha_{K}$. Normalisation by the correlation time
$\tau$, which to first approximation is assumed to be common to both magnetic
and kinetic fluctuations, yields $\alpha$ in dimensional form as speed.

[FAG: I haven’t caught them all, but in general results from citations are in
the present tense, unless since the results have proven to have changed. For
example, in 2006 Anvar finds …, because we still find this to be the case.]
Mean field dynamos produce a large scale field with magnetic helicity
corresponding to the mutual linkage of the poloidal and toroidal magnetic
fields. The dynamo must, then, also produce a small scale magnetic field with
opposite magnetic helicity. Unless the small scale magnetic field can be
transported out of the system it quenches the $\alpha$-effect and
simultaneously the mean field dynamo. Anvar_2006_helicity shows that without
advection removing the small scale field there is catastrophic quenching of
the mean field.

[FAG: Please revise the next paragrah to accurately explain the Figure, how it
compares to Devika’s and the meaning/formulation of the various subscripts.
Comment on how it varies or not between models.] In Fig. 8 the time-latitude
diagrams for the componenets $\alpha$-effect are plotted for Model XXXX. The
helicity during the linear stage of the dynamo is conentrated about the
midplane with sign consistent with vertical flows driving helical magnetic
fields of opposite sign either side of the midplane. Figure 13 of dummy
presents the kinetic helicity for model with solid body rotation, it shows the
kinetic helicity has an anomalous sign at large $z$ away from the midplane,
similarly in our case this the kinetic helicity generated has anomalous sign
away from the midplane.

Figure 8:  Time-latitude diagrams of composition of horizonatlly averaged
$\alpha$-effect for Model XXXX. [FAG: please arrange this figure as a single
column with sharex=True, colorbar labels on the side and aspect ratio of each
plot similar to Fig. 7. Use fontsize to ensure text are as large as caption
fonts. Label y-axis. ]

##  4 Mean-field $\alpha^{2}$-dynamo model

In order to understand the structure of the dynamo and the buoyancy
instability interactions, we seek to reconstruct the DNS solutions with a
simpler mean-field model depending only on vertical location $z$ and time $t$.
Our objectives are to replicate the relevant growth rates in the magnetic
field and the shape, amplitude and period of the reversals in the horizontal
field components. The reversals are common to all simulations that reach
saturation and exhibit this similar vertical expansion of the field. Model
R10h2 has the fastest growth and so we begin by fitting a mean-field model to
its parameters. Can we attribute this reversal to the nonlinear dynamo and the
response of the magnetic field to the magnetic buoyancy instability.

To explore this we employ the use of simple Runge-Kutta algorithm which solves
the induction equation for the horizontal components of the magnetic field,
since the vertical component of the magnetic field can be obtained from the
solenoidality condition, and the vertical momentum equation. The vertical
momentum equation must contain both pressure and gravitational forces in order
to include buoyancy. Likewise the induction equation must contain the EMF
which models the effect of the turbulent velocity on the magnetic field.
Although the system of equations we solve is isothermal we can draw parallels
to the Boussinesq approximation. The Boussinesq approximation is used to
simplify the equations of motion when variations in density are due primarily
to variations in temperature, in our situation the variation in density is
purely due to the magnetic field. Making considerations of a fluid parcel in
pressure balance with its surroundings then the variation in pressure is

Saturation effects such as the associated removal of the magnetic field by
vertical velocity will eventually halt field growth in a more realistic
system. [FAG: For an ideal gas I believe we require $R/\mu T_{0}$.]

$\displaystyle\delta\rho=\frac{2\langle{{\mbox{\boldmath$B$}}}\rangle{{\mbox{\boldmath$b$}}}+{{\mbox{\boldmath$b$}}}^{2}}{8\pi
R\mu^{-1}T_{0}}=\frac{{{\mbox{\boldmath$b$}}}^{2}}{8\pi R\mu^{-1}T_{0}}$ (20)

where $T_{0}$ is the mean temperature of the system, the gas constant
$R=k_{B}/(\mu m]$, $k_{B}$ is the Boltzmann constant and $m$ is the atomic
mass of the gas. The first term is dropped since the average of the
fluctuating field will be zero. The vertical momentum equation and horizontal
induction equation with the EMF are given by

$\displaystyle\frac{\partial u_{z}}{\partial t}$
$\displaystyle=-u_{z}\frac{\partial u_{z}}{\partial
z}-\frac{\partial}{\partial
z}\left(\frac{\left|{{\mbox{\boldmath$B$}}}\right|^{2}}{8\pi\rho}\right)+\frac{\delta\rho}{\rho_{0}}g_{z}+\nu\nabla^{2}u_{z}$
(21) $\displaystyle\frac{\partial B_{x}}{\partial t}$
$\displaystyle=-R_{\alpha}\frac{\partial}{\partial
z}(\hat{\alpha}(z)B_{y})+u_{z}\frac{\partial B_{x}}{\partial
z}+\frac{\partial^{2}B_{x}}{\partial z^{2}}$ (22) $\displaystyle\frac{\partial
B_{y}}{\partial t}$ $\displaystyle=R_{\alpha}\frac{\partial}{\partial
z}(\hat{\alpha}(z)B_{x})+u_{z}\frac{\partial B_{y}}{\partial
z}+\frac{\partial^{2}B_{y}}{\partial z^{2}}$ (23)

where $\rho_{0}$ is the mean density. These equations are solved with vacuum
boundary conditions,

$\displaystyle\left.\frac{\partial B_{x}}{\partial
z}\right|_{z=0}=\left.\frac{\partial B_{y}}{\partial
z}\right|_{z=0}=\left.\frac{\partial u_{z}}{\partial z}\right|_{|z|=z_{0}}=0,$
(24)
$\displaystyle\quad\left.B_{x}\right|_{|z|=z_{0}}=\left.B_{y}\right|_{|z|=z_{0}}=0$
(25)

[FAG: Is this correct, is the problem not that without $\alpha$ the buoyancy
will not saturate the dynamo? Please clarify the role of these terms.] For
this simple mean-field model the magnetic field would be removed by magnetic
buoyancy which would eventually result in saturation, but the effects of the
magnetic helicity are ignored. To address this we can employ an algebraic
quenching of the $\alpha$-effect

$\hat{\alpha}=\dfrac{\alpha(z)}{1+(B/B_{0})^{2}}$ (26)

In the 1D model we can see that the $\alpha$-effect causes the x and y
components of the magnetic field have opposite sign to each other resulting in
the two stable configurations seen before the saturation phase. Once we enter
the saturation/non-linear phase the dynamo has “switched off" and both the x
and y components acquire the same sign. When the buoyancy advects the field
away, the dynamo reactivates and begins to generate a magnetic field based on
the strongest components of the $\alpha$ term which is
$\dfrac{\partial\alpha}{\partial z}B$ since both components now have the same
sign one of the components will have to switch sign to return to the stable
growing configuration from the linear stage, in the case of the reference
state the first, to switch sign is $B_{x}$ the change is quick due to strong
dynamo action, this causes an overdampening resulting in an oscillatory steady
state.

There is a possibility for the magnetic field to grow with $\left|z\right|$ by
the vertical velocity $V_{B}$, due to this buoyancy and because the
$\alpha$-effect intensifies with height. The expansion of the field seems
inevitable as such a strong magnetic field cannot be kept within the galactic
disc and should expand into the corona. Simulations by (1999MOSS) describe
buoyancy driven dynamos and also find the field expands into the outer region
of the simulation domain and retains similar levels of strength both in the
disc midplane and near the boundaries. The expansion of the field follows an
exponential distribution.

$z\lvert_{\langle B\rangle_{xy}=0}\sim\exp(t/\tau_{z})$ (27)

$\tau_{z}\approx 5$ in the region $z>0$ for $B_{y}$ and $B_{x}$ and
$\tau_{z}\approx 4$ in the region $z>0$ for $B_{y}$ and $B_{x}$ the
corresponding timescales do not match with either the timescale of diffusion,
the dynamo or the instability. Larger values of $\alpha_{0}$ have larger
growth rates for the magnetic field and steeper spatial gradients causing the
expansion to intensify with larger values of $\alpha_{0}$ and become weaker
for smaller vales of $\alpha_{0}$. Increasing the scale height of the dynamo
$h_{\alpha}$, does not affect the rate of the expansion as long as the dynamo
number $R_{\alpha}$ remains the same. where $\alpha(z)$ is the same as in
equation 2 and $B_{0}$ is the equipartition value of the magnetic field, when
including $\alpha$ quenching the we can more precisely control the period of
the resulting non-linear oscillations of the magnetic field and match them
with the full simulation. The value of $B_{0}$ is free to choose however to
match the oscillation periods the saturation height becomes to low. To remedy
this we consider that the $x$ and $y$ components of the magnetic field have
their own saturation mechanisms depending on their individual field strengths
this results in the right hand column of figure 9 where the period and
morphology of the different models is easily seen. The treatment of having the
$\alpha$ quenching depend on individual components is also seen in
1993_LLKit_GRud where $\alpha$ quenching of a strong magnetic field is
considered by splitting the quenching effects over different components.
Moreover the saturation level of the models with and without the modified
quenching are shown in figure 10. The saturated values of the model and the
are of the same order of magnitude with the DS saturating around $3\mu G$ and
the model saturating around $1\mu G$.

Figure 9, shows the oscilations in the full simulation the Runge-Kutta model
and the Runge-Kutta model with quenching for $B_{x}(z,t)$ and $B_{y}(z,t)$ are
normalized to $1$. We observe that the same oscillationsssss are present in
both the full and the 1D simulations

Figure 9: (left) Temporal evolution of $\langle\hat{B}_{x}\rangle$ and
$\langle\hat{B}_{y}\rangle$. showing the clear oscillation of the magnetic
field and subsequent expansion into the outer regions of the simulation
domain. (middle) 1D Runge-Kutta model without quenching, the period of
oscillation for this model is much shorter than the . (right) 1D Runge-Kutta
model with $\alpha$ quenching, the oscillation period matches that found in
the . The period of the non-linear oscillations depends on the the parameter
$B_{0}$. All of these quantities have been normalized to unity, and the
parameters used are those for the reference state. [FAG: Excellent work. Pleae
split this into two plots, each a column for $B_{x}$ and $B_{y}$, with
sharex=True, and aspect ratio of each plot similar to Fig. 7. Use fontsize to
ensure text are as large as caption fonts. ] Figure 10: The evolution of the
r.m.s magnetic field for the , the Runge-Kutta model (1DRa10ha200) and the
quenched Runge-Kutta model (1DQRa10ha200) are using the same parameters as the
reference state. The quenched model saturates at a value of around $1$ $\mu G$
and the saturates around a value of $3$ $\mu G$.

##  5 Conclusions

[FAG: This belongs to the discussion and conclusions.] The STF model
(zeldovich1983magnetic) for exponential growth of magnetic fields explains how
magnetic flux tubes are stretched and twisted forming a loop which is folded
onto itself as to retain the flux but increase the field strength, this
picture ends with the diffusive effects acting on the folded tube to create a
singular un-knotted flux tube with a stronger magnetic field. With lower
values of $R_{\alpha}$ the STF mechanism is weaker meaning the tangling of the
field takes place over a longer period.

[FAG: Move to discussion/conclusions] (Fletcher_2011) found that the dynamo
number $R_{\alpha}$ is limited by the turbulent velocity scale,
$R_{\alpha}<3h/l_{0}=15$, where $h=0.5$ kpc is the disk scale height and
$l_{0}=0.1$ kpc is the length scale of turbulence.

[FAG: Move to the discussion or the methods] In section $4$ of (Anvar_1987),
it is shown that when the maximum helicity is close to the boundary the
critical value of $R_{\alpha}$ will increase, moreover we observe the growth
rates for similar values of $R_{\alpha}$ become smaller. The choice of values
$h_{\alpha}$ was to allow the dynamics to occur on a reasonable time frame,
moreover $h_{\alpha}$ also controls the strength of the instability for larger
values of $h_{\alpha}$, the instability is weaker. During the linear stage of
the dynamo the velocity is at a statistically steady state.

[FAG: Move to disc/conc] Section 9.4 of (Brandenburg_2005) demonstrates that
the saturation level of the $\alpha^{2}$ and $\alpha\omega$ dynamos are
dependant on the dynamo strength and that the final saturation level of the
$\alpha^{2}$ dynamo can be estimated by using the following expression

$B_{\rm fin}=\left(\frac{C_{\alpha}}{C_{\alpha,{\rm
crit}}}-1\right)^{1/2}\left(1+R_{m}^{-1}\right)^{1/2}B_{\rm eq}$ (28)

with,

$C_{\alpha}=\alpha_{k}/(\eta k_{m}),\quad C_{\alpha,{\rm
crit}}=\alpha_{crit}/(\eta k_{m})$ (29)

where
$k_{m}^{2}=\langle\langle{{\mbox{\boldmath$j$}}}\rangle\cdot\langle{{\mbox{\boldmath$B$}}}\rangle\rangle/\langle\langle{{\mbox{\boldmath$A$}}}\rangle\cdot\langle{{\mbox{\boldmath$B$}}}\rangle\rangle$
is the characteristic wavenumber of the mean-field and $\langle.\rangle$ here
denotes some suitable averaging and $\alpha_{crit}=\eta_{T}\tilde{k}_{m}$.
This analytic formula seems to agree with our simulation with the final
saturation level of the dynamo being the same order of magnitude for many of
the different simulations, for instance with $h_{\alpha}=0.2$ kpc,
$R_{\alpha}=10$ and $R_{\alpha}=5$ the analytic results are $B_{fin}\approx
10^{-1}$ and $10^{0}$ and the final saturated values from the simulation are
$0.3516$ and $1.1995$ respectively.

We have used an imposed $\alpha$-effect to generate a magnetic field in which
the Parker instability arises naturally. We observed how the generated
magnetic field resists the characteristic sign of the Parker instability in
which the density scale height decreases and the gas forms into pockets, in
this system the increases of magnetic pressure from he dynamo resists the
reformation of the density and during the non-linear/saturation phase of the
dynamo the magnetic density scale height increases. The observed effects were
significantly different to what was seen in previous work with an imposed
magnetic field. Most notably is the arrival of non-linear oscillations in the
azimuthal components of the magnetic field, due to the combined effects of
magnetic buoyancy and the $\alpha$-effect. Using 1D model we showed how the
oscillations in the non-linear phase are due to the combined effects buoyancy
and the $\alpha$-effect leading to an oscillatory steady state for the
azimuthal magnetic field. Moreover these effects occur independent of the
Parker instability which in the transitional phase simply enhances the dynamo
action as was suggested by (1992PARKER).

## Acknowledgements

We would like to thank Dr Graeme Sarson for useful discussions about the
effect of quenching on the 1D model.

## Data Availability

The raw data for this work were obtained from numerical simulations using the
open source PENCIL-CODE available at https://github.com/pencil-code/pencil-
code.git. The derived data used for the analysis given in the paper is
available on request from the corresponding author.

###  .1 The phases of the instability

[FAG: moved this plot to appendix to be retained in this paper and referenced
in the main text. Other items in the appendix may be used or retained for
Yasin’s thesis?]

Figure 11: Different components for Model XXXX of $B$ advancing in time. After
the saturation of the magnetic field at $t\geq 1.5Gyr$ there is the appearance
of large helical structures in the $B_{x}$ and $B_{y}$ components.

####  .1.1 Linear phase

[FAG: Where do we show this and is it important in the linear phase of the
dynamo?] All three components of the magnetic field are equally subjected to
the $\alpha$-effect and yet ${{\mbox{\boldmath$B$}}}_{z}$ has a larger growth
rate than the ${{\mbox{\boldmath$B$}}}_{x}$ and ${{\mbox{\boldmath$B$}}}_{y}$
components whose rate of growth is similar, this may be attributed to the
magnetic field in the $z$ direction being modified by magnetic buoyancy, which
could account for the increased growth. However the difference in the growth
rates of the different components is negligible.

$\rho g_{z}=\dfrac{\partial}{\partial z}(P_{\rho}+P_{B})$ (30)

Equation 30 shows the balance of the magnetic field, density and gravitational
forces, and figure 6 shows the density variation for three different
simulations with the parameters varied around the reference state. Initially
the density has a Gaussian distribution and the density pressure in the
midplane is stronger than that in outer region the growth of the magnetic
field in the midplane is suppressed. amplification of the magnetic field due
to the dynamo leads to increased magnetic pressure. To retain pressure balance
along the layer the density is pushed out of regions with a strong local
magnetic field. This perturbation of the density leads to the onset of the
instability and also the saturation of the magnetic field. For smaller
${\alpha_{0}}$ dynamo action is much weaker which leads to a comparatively
weaker field consequently the density perturbation from the magnetic field is
weaker leading to a weaker instability. For larger values of $h_{\alpha}$ the
dynamo process produces a weaker magnetic field which again leads to a weaker
instability.

Figure 12: The vertical variation of the horizontally averaged force
components for the reference state all normalized to unity by their maximums.

##  Appendix A Calculating the magnetic Reynolds number

We obtain the magnetic Reynolds number of the flow from the second order
structure functions

$D(\ell)=\langle\lvert
f({{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$\ell$}}})-f({{\mbox{\boldmath$x$}}})\rvert\rangle_{\ell}$
(31)

where $x$ is a point in the $x-y$ plane and $\ell$ is a horizontal offset,
using the Taylor-microscale, as outlined in (Hollins_2017) and in section 13.1
of (MoninAndYaglom).

The chosen horizontal boundary conditions are periodic so the structure
functions are only calculated for the regions $0<\ell_{x}<L_{x}/2$
,$0<\ell_{y}<L_{y}/2$, where $L_{x}$ and $L_{y}$ are the half length of the
box in our case $2$ kpc. Using the structure function $D(\ell)$ we can
construct the autocorrelation function $C(\ell)$

$C(\ell)=1-\dfrac{D(\ell)}{2\sigma^{2}}.$ (32)

When $\ell\rightarrow 0$, then near the origin

$C(\ell)\approx 1-\left(\dfrac{\ell}{\lambda}\right)^{2}.$ (33)

$\lambda$ is the Taylor microscale and is given by

$\lambda\approx 3\ell_{0}\rm{Re}^{-1/2}$ (34)

To estimate the Reynolds number the spatial structure function is calculated
and fitted by a function of the form

$D(\ell)=2\sigma^{2}\left[1-\exp\left(-\dfrac{\ell^{2}}{2L_{0}^{2}}\right)\right],\quad\ell_{0}=\sqrt{\dfrac{\pi}{2}}L_{0}$
(35)

and then a second degree polynomial is fitted near the origin of $C(\ell)$ to
determine the value of $\lambda$ which can be inserted into 34 to determine
the value of the Reynolds number. If the procedure is preformed on the
velocity field we obtain the Reynolds number, but if the same procedure is
preformed on the magnetic field then we obtain the magnetic Reynolds number.

Figure 13: Calculating the magnetic Reynold number from the the fluctuating
magnetic field for various dynamo scale heights.

##  Appendix B Averaging methods

In similar studies of the galactic magnetic fields and their evolution a range
of different averaging methods are used to in order to study the behaviour of
the magnetic field, the main methods we employ are horizontal averaging and
smoothing by a Gaussian kernel. Both these averaging techniques have their own
merits, horizontal averaging is useful for analysis of large scale effects,
however horizontal averaging will destroy any significant spatial structure,
which is when Gaussian smoothing can be used to analyse the mean and
fluctuating field averaged over some appropriate scaling length, we follow the
method outlined in GENT_ISM_II and (Hollins_2017). Identifying the scale
separation in mean and random fields can be exceptionally difficult as often
there are no signs of scale separation in the Fourier spectrum with no local
minima visible, using the standard approach (e.g. horizontal averaging) the
system would be said to lack scale separation. However the filtering approach
to averaging has been used in GENT_ISM_II to identify clear scale separation
in the mean and random fields. The power spectrum of $B$ is
$M(k)=k^{-2}\langle|\mathcal{F}(k)|\rangle_{k}$ for spherical shells of
thickness $\delta k$ at radii $k=|{{\mbox{\boldmath$k$}}}|$, from
$\mathcal{F}({{\mbox{\boldmath$k$}}})=\int_{V}{{\mbox{\boldmath$B$}}}({{\mbox{\boldmath$x$}}})\exp(-2\mathrm{\pi
i}{{\mbox{\boldmath$k$}}}\cdot{{\mbox{\boldmath$x$}}})d^{3}{{\mbox{\boldmath$x$}}}$.
The power spectrum for the mean and random magnetic fields are found by using
Gaussian smoothing, the process is covered in detail in chapter $7.2$ of
(shukurov_subramanian_2021).

The magnetic field $B$ can be decomposed into two parts. The first part
${{\mbox{\boldmath$B$}}}_{l}$, the part averaged over the length-scale $l$,
and the fluctuation field ${{\mbox{\boldmath$b$}}}_{l}$.

${{\mbox{\boldmath$B$}}}={{\mbox{\boldmath$B$}}}_{l}+{{\mbox{\boldmath$b$}}}_{l},\quad{{\mbox{\boldmath$B$}}}_{l}=\langle{{\mbox{\boldmath$B$}}}\rangle_{l},\quad{{\mbox{\boldmath$b$}}}_{l}={{\mbox{\boldmath$B$}}}-{{\mbox{\boldmath$B$}}}_{l}$
(36)

The mean-field ${{\mbox{\boldmath$B$}}}_{l}$ is calculated using volume
averaging with a Gaussian kernel:

$\displaystyle\langle{{\mbox{\boldmath$B$}}}\rangle_{l}({{\mbox{\boldmath$x$}}})$
$\displaystyle=\int_{V}{{\mbox{\boldmath$B$}}}({{\mbox{\boldmath$x$}}}^{\prime})G_{l}({{\mbox{\boldmath$x$}}}-{{\mbox{\boldmath$x$}}}^{\prime})d^{3}{{\mbox{\boldmath$x$}}}^{\prime}$
(37) $\displaystyle G_{l}({{\mbox{\boldmath$x$}}})$
$\displaystyle=(2\mathrm{\pi}l^{2})^{-3/2}\exp{(-{{\mbox{\boldmath$x$}}}^{2}/2l^{2})}$
(38)

With $V$ as the volume of the computational domain. The convolution of
${{\mbox{\boldmath$B$}}}_{l}$ with the Gaussian kernel preserves the
solenoidality of the ${{\mbox{\boldmath$B$}}}_{l}$ and
${{\mbox{\boldmath$b$}}}_{l}$ fields and the three dimensional structure of
the field which is not the case for horizontal averaging. However since
repeated averaging changed the result
$\overline{\overline{{{\mbox{\boldmath$B$}}}}}\neq\overline{{{\mbox{\boldmath$B$}}}}$,
so Reynolds rules for averaging no longer apply and so the definitions for
averaged quantities should take the more general form proposed in
(germano_1992).

The length scale $l$ across which we average is a free parameter. For small
$l$ the mean-field and total field will be almost identical, and for large $l$
the mean-field and total field will be largely the same, so in order to choose
an appropriate choice of $l$ we aim to choose an $l$ such that the
intersection of ${{\mbox{\boldmath$B$}}}_{l}$ and
${{\mbox{\boldmath$b$}}}_{l}$ occur in between the energy-range scales for the
mean $2\mathrm{\pi}L_{B}^{-1}$ and random fields $2\mathrm{\pi}L_{b}^{-1}$,
where $L_{B}$ and $L_{b}$ are the integral scales given by

$L_{\phi}=\frac{3\pi\int^{\pi\Delta}_{2\pi/D}k^{-1}M(k)dk}{4\int^{\pi\Delta}_{2\pi/D}M(k)dk}$
(39)

Where $\Delta$ is the numerical grid separation and $D$ is the size of the
domain $\phi$ as some field, in our case either ${{\mbox{\boldmath$B$}}}_{l}$
or ${{\mbox{\boldmath$b$}}}_{l}$ MoninAndYaglom. Moreover the maxima for each
spectra should lie outside this range yet still relatively close to each of
the energy range scales. Using this method we have been able to isolate the
scale separation for different simulations at the kinematic and saturation
phases of the dynamo. Figure of 14 is the power spectrum for the magnetic
field for the reference state. During the kinematic phase of the dynamo the
mean-field is still growing as such the scale separation between the mean and
fluctuating field is smaller. The mean-field has a broad maximum at around $2$
pc and the fluctuation field has maximum at around $6-7$ kpc. In the
saturation phase of the magnetic field the power spectrum for the mean-field
shifts to larger scales due to the onset of the instability which enhances the
growth of the mean magnetic field, with the value at $1.8$ Gyrs corresponding
to a wavelength of around $1$ kpc, which is the characteristic length scale of
the instability.

Figure 14: Power spectrum $M(k)$ for simulations with the reference state.
(Top) Radially averaged power spectrum for the reference state during the
linear phase of the dynamo. (Bottom) During the saturation phase of the
dynamo.

##  Appendix C Dependence of Results on Box Width

The mean-field dynamo and fluctuation dynamo problems are often treated
separately based on the scales they act on, both temporally and spatially,
however they both describe a single dynamo system which is capable of
generating both large and small scale magnetic fields. To capture both large
and small scale effects special care must be placed on the resolution and box
size, and we must take care to see if the numerical domain of the simulation
is affecting the dynamics of the system. So we investigate the effects of
changing the width of the box, from
$(L_{x},L_{y})=(0.5\text{kpc},0.5\text{kpc})$ to
$(L_{x},L_{y})=(8\text{kpc},8\text{kpc})$. In all box sizes, the helical
structure of the system during the saturation phase remains unchanged and in
boxes with larger x and y extent more vortices to form. The wavelength and
wavenumber of the instability also remains similar with the initial simulation
showing wavelengths of $\lambda=0.2$kpc with a wave number of k $=5$ and the
box of width $8$kpc having a wavelength of $\lambda=0.2$kpc and a wavenumber
of k $=5$. The increased width of the box causes a shift in the spectra to
higher wavenumbers but, the shift is small and the wavelength of the
instability remains largely the same. For the smaller box sizes some of the
detail is being lost however as we increase the box size the profile of the
spectra begins to converge to where the $4$kpc and $8$kpc boxes have the same
profile for the their spectra although slightly shifted. Also of note is that
in figure 15 a) the maxima for all but the box with $L=0.5$ are in
approximately the same position around $log(k)=7$kpc-1, meaning that during
the saturation phase of the dynamo changing the box width is not changing the
dynamics too dramatically except for in the case for $L=0.5$ where the smaller
box size is affecting our results. During the kinematic phase of the dynamo
the shown in 15 b) where we see that the growth rates for all the spectra are
very similar and share a maxima at around $log(k)=30$kpc-1, which corresponds
to a wavelength of around $0.2$kpc, these modes are clearly seen in the
snapshots where there are variations of around $0.2$kpc visible in the
magnetic field. Although it is clear that some detail is still being lost for
a box size of $L\leq 2$kpc so we opt to use the the larger box of $4$kpc.

Figure 15: a) The evolution of the magnetic power spectrum $M(K)$ for boxes of
different widths (L). The spectra are shown for 5 different box sizes with
$L=0.5\text{kpc},1\text{kpc},2\text{kpc},4\text{kpc},8\text{kpc}$, at
$t=700$Myr. The instability sets in around $600$Myr, so these measurements
correspond to the saturation phase of the dynamo. b) The same spectra captured
above for the different boxes at an earlier time during the kinematic phase of
the dynamo.  Figure 16: Log-linear plot of the kinetic and magnetic energy
densities and also an exponential with twice the growth rate shown in table 2.
The growth of the field is linear initially when kinetic energy is small and
as the magnetic field gains strength, the Lorentz force becomes stronger
driving larger velocities this ends the linear phase of the dynamo and
eventually we enter the saturation phase of the dynamo.

##  Appendix D Calculating the parity of the magnetic field

The parity of the magnetic field is calculated using the method from
(LFSR2015), in which the magnetic field is broken down into symmetric
components

$\displaystyle{{\mbox{\boldmath$B$}}}_{x,y,z}^{(s)}=\frac{1}{2}\left({{\mbox{\boldmath$B$}}}(x,y,z)+{{\mbox{\boldmath$B$}}}(x,y,-z)\right)$
(40)

and anti-symmetric components,

$\displaystyle{{\mbox{\boldmath$B$}}}_{x,y,z}^{(a)}=\frac{1}{2}\left({{\mbox{\boldmath$B$}}}(x,y,z)-{{\mbox{\boldmath$B$}}}(x,y,-z)\right)$
(41)

The total energy associated with each magnetic field component is then
determined using
$E_{x,y,z}^{(s/a)}=\int\left({{\mbox{\boldmath$B$}}}_{x,y,z}^{(s/a)}\right)^{2}/8\pi
dV$ , the parity of the field is then given by,

$\displaystyle
P_{x,y,z}=\frac{E_{x,y,z}^{(s)}-E_{x,y,z}^{(a)}}{E_{x,y,z}^{(a)}+E_{x,y,z}^{(s)}}$
(42)

Figure 17: Parity of the $B_{x}$, $B_{y}$ and $B_{z}$ components of the
magnetic field for a dynamo scale height $h_{\alpha}=0.2$ kpc, with
$\alpha_{0}$=10 kpc Gyr-1.

The resulting parities for each field for reference state are shown in figure
17. The saturation of the magnetic field occurs at around $1.5$ Gyrs and after
that we see a slight preference in $B_{z}$ for symmetry and a slight
preference in $B_{x}$ and $B_{z}$ for anti-symmetry however these preferences
are weak and we can conclude the magnetic field has no preference for either
even or odd modes, however this slight preference for anti-symmetry is inline
with the preferred field configuration for the magnetic field under the
influence of Parker instability. The instability however is a large scale
effect and will mainly affect the mean magnetic field and so the same
procedure can be done on the mean and fluctuating fields and is shown in
figure 18, the fluctuating field is asymmetric and shows no preference for
either even or odd modes this is to be expected due tot he random nature of
the fluctuating field. Initially again the mean-field is asymmetric, however
when the instability sets in and saturates the magnetic field then the mean-
field shows a preference even modes.

Figure 18: Same as figure 17 but with each component broken up into the mean
and fluctuating parts.

##  Appendix E Mechanism of the oscillations

Small-scale inhomogeneities due to the random fluctuating field at the
nonlinear phase of the dynamo cause small-scale buoyancy forces correlated
with the mean magnetic field. Captured in the EMF of the mean-field induction
equation, the buoyant magnetic field rises and expands into the outer regions
of the disk. As the field expands outward it is replaced by a field of equal
strength but crucially opposite sign. The amplitude of the resulting wave
follows a steady state evolution. Despite these oscillations the rapidly
changing magnetic gradients ensure the magnetic diffusivity inhibits growth of
the magnetic field so that the field strength remains steady.

These oscillations arise through the combined action of magnetic buoyancy and
the $\alpha$-effect. The dynamo generated magnetic field once saturated
becomes buoyant and the magnetic field is buoyantly advected away from the
midplane. Where buoyancy results in loss of total local pressure there is an
increase in local velocity and subsequently also kinetic helicity. Local
kinetic helicity increases eventually surpass the magnetic helicity, thus
“switching on” the $\alpha$-effect to produce large scale magnetic field that
replaces the field lost to buoyancy (1999MOSS; Parker1979). The newly active
$\alpha$-effect generates magnetic field of opposite direction.

The $\alpha$-effect contribution to evolving $B_{x}$ is therefore given by

$\dfrac{\partial\alpha}{\partial z}B_{y}+\alpha\dfrac{\partial B_{y}}{z}.$
(43)

$\alpha$ is sinusoidal in $z<\left|h_{\alpha}\right|$, so zero at the midplane
$z=0$, where the reversals originate. So these oscillations are driven
primarily by the $B_{y}\partial_{z}\alpha$. This is confirmed by the last two
rows of Fig. 19, with similar conclusions for the generaton of $B_{y}$ in the
second column of Fig. 19.

Figure 19:  [FAG: I cannot read these colorbar labels, put them at the side,
which is the pyplot default] Horizontal averages of the nonlinear period of
the reference state for the $x$ and $y$ components of (Row 1) the magnetic
field, (Row 2) the generation term in the induction equation, (Row 3),
$\dfrac{\partial\alpha}{\partial z}B_{x}$, $\dfrac{\partial\alpha}{\partial
z}B_{y}$, (Row 4) $\dfrac{\partial B_{x}}{\partial z}\alpha$ and
$\dfrac{\partial B_{y}}{\partial z}\alpha$.

Figure 20 shows the absolute values of the kinetic and magnetic helicity
averages over the subvolume $\lvert z\rvert<0.1$kpc.

Figure 20: Kinetic and magnetic helicity in the region $|z|<0.1$kpc. When the
magnetic helicity is above kinetic helicity the $\alpha$ effect “switches
off”.

##  Appendix F Oscillatory solutions to the induction equation

Oscillatory solutions of the induction equation with the $\alpha$-effect are
seen in (Anvar_1987). In this section we show that the induction equation in
the non-linear regime permits these types of oscillations. We follow the same
method outlined in (Anvar_1983), starting with the induction equation with the
non-linear term
$\left|{{\mbox{\boldmath$B$}}}\right|^{2}\nabla{{\mbox{\boldmath$B$}}}$ which
models the back reaction of the EMF on the induction equation.

$\frac{\partial{{\mbox{\boldmath$B$}}}}{\partial
t}=\nabla\times(\alpha{{\mbox{\boldmath$B$}}})+\left|{{\mbox{\boldmath$B$}}}\right|^{2}\nabla{{\mbox{\boldmath$B$}}}+\nabla^{2}{{\mbox{\boldmath$B$}}}$
(44)

Using the symmetry relation $B_{r}=-iB$ to reduce the equations in the
azimuthal direction down to just one equation in $B$.

$\frac{\partial B}{\partial t}=\nabla\times(\alpha
B)+\left|B\right|^{2}\frac{\partial B}{\partial
z}+\frac{\partial^{2}B}{\partial z^{2}}$ (45)

An solution which oscillates with period $\omega$ and has a phase shift
dependant on $\alpha$ is employed

$B=\psi(z,t)\exp\left(\frac{i}{2}\int\alpha dz+i\omega t\right)$ (46)

Subbing this into the induction equation above, splitting into real and
imaginary parts and then rearranging each equation to make the term with the
highest derivative the subject gives us a timestepping equation for both
$\psi^{\prime}$ and $\psi^{\prime\prime}$

$\displaystyle\psi^{\prime}$
$\displaystyle=\frac{\omega+R_{\alpha}\alpha^{\prime}-\frac{\alpha^{\prime}}{2}+\frac{|\psi|^{2}\alpha}{2}}{\alpha-
R_{\alpha}\alpha}$ (47) $\displaystyle\psi^{\prime\prime}$
$\displaystyle=|\psi|^{2}\psi^{\prime}-\psi\left(\frac{R_{\alpha}\alpha^{2}}{2}-\frac{\alpha^{2}}{4}\right)$
(48)

These two equations are then solved numerically using the shooting method with
using quadrupole boundary conditions to replicate the structure in of the
field in both the 1D and DNS simulations.

$\displaystyle\frac{\partial B}{\partial
z}\bigg{\rvert}_{z=0}=0,\quad\frac{\partial B}{\partial
z}\bigg{\rvert}_{z=z_{0}}=0$ (49)

This leaves $\omega$ as a free parameter, matching it together with the DNS we
can obtain a solution with a similar period to the reference state, shown in
figure 21.

Figure 21: Plot of the solution to the non-linear equations 47 and 48. Solved
using the shooting method in the upper half z plane.