# MHD in a cylindrical shearing box II: Intermittent Bursts and Substructures
in MRI Turbulence

Takeru K. Suzuki School of Arts & Sciences, The University of Tokyo, 3-8-1,
Komaba, Meguro, Tokyo 153-8902, Japan<EMAIL_ADDRESS>Department
of Astronomy, The University of Tokyo, 7-3-1, Hongo, Bunkyo, Tokyo, 113-0033,
Japan  Komaba Institute for Science, The University of Tokyo, 3-8-1 Komaba,
Meguro, Tokyo 153-8902, Japan

###### Abstract

By performing ideal magnetohydrodynamical (MHD) simulations with weak vertical
magnetic fields in unstratified cylindrical shearing boxes with modified
boundary treatment, we investigate MHD turbulence excited by magnetorotational
instability. The cylindrical simulation exhibits extremely large temporal
variation in the magnetic activity compared to the simulation in a normal
Cartesian shearing box, although the time-averaged field strengths are
comparable in the cylindrical and Cartesian setups. Detailed analysis of the
terms describing magnetic-energy evolution with “triangle diagrams”
surprisingly reveals that in the cylindrical simulation the compression of
toroidal magnetic field is unexpectedly as important as the winding due to
differential rotation in amplifying magnetic fields and triggering
intermittent magnetic bursts, which are not seen in the Cartesian simulation.
The importance of the compressible amplification is also true for a
cylindrical simulation with tiny curvature; the evolution of magnetic fields
in the nearly Cartesian shearing box simulation is fundamentally different
from that in the exact Cartesian counterpart. The radial gradient of epicycle
frequency, $\kappa$, which cannot be considered in the normal Cartesian
shearing box model, is the cause of this fundamental difference. An additional
consequence of the spatial variation of $\kappa$ is continuous and ubiquitous
formation of narrow high(low)-density and weak(strong)-field localized
structures; seeds of these ring-gap structures are created by the compressible
effect and subsequently amplified and maintained under the the marginally
unstable condition regarding “viscous-type” instability.

accretion disks — burst astrophysics — MHD — MHD simulations — protoplanetary
disks — turbulence

##  1 Introduction

The local shearing sheet/box model (Goldreich1965; Narayan1987; Hawley1995;
Latter2017) is a strong tool to study various physical properties of
differentially rotating systems. To date, local shearing box simulations have
been applied to a variety of astrophysical problems, such as
magnetohydrodynamical (MHD hereafter) turbulence (e.g., Matsumoto1995;
Stone1996; Kawazura2022) excited by magnetorotational instability (MRI
hereafter; Velikhov1959; Chandrasekhar1961; Balbus1991), dynamical and thermal
properties of protoplanetary disks (e.g., Kunz2013; Mori2017; Pucci2021) and
accretion disks around compact objects (e.g. Hirose2009; Dempsey2022), heating
accretion disk coronae (Io2014; Bambic2023arxiv), driving disk winds and
outflows (Suzuki2009ApJ; Suzuki2010ApJ; Bai2013; Lesur2013; Fromang2013),
acceleration of non-thermal particles (Hoshino2015; Kimura2016; Bacchini2022),
and amplification of magnetic fields in compact objects (Masada2012ApJ;
Guilet2022).

One of the drawbacks in the normal Cartesian shearing box model is that the
positive radial direction is not well defined because the symmetry with
respect to the $\pm x$ directions is assumed. As a result, the angular
momentum is not defined and the radial accretion flow cannot be properly
captured in the Cartesian shearing system. Therefore, mass accretion rate
cannot be directly measured in numerical simulations but is inferred from the
sum of Maxwell and Reynolds stresses. In addition, the basic physical
properties of magnetic-field amplification is considered to be qualitatively
different in linear shear flows in Cartesian coordinates and more realistic
differentially rotating flows in cylindrical coordinates (Ebrahimi2016MNRAS).
Although these problems are not inherent in global simulations (e.g.,
Armitage1998; Machida2000; Flock2013; Suzuki2014ApJ; Bethune2017; Takasao2018;
Zhu2020; Jacquemin-Ide2021), their numerical cost is generally more expensive.
Thus, if one can construct a modified model that takes into account global
effects in the local approach by breaking the $\pm x$ symmetry, it will be an
extremely efficient tool to investigate differentially rotating systems more
appropriately.

To this end, there have been several attempts that try to bridge the local and
global concepts over the past few decades (Brandenburg96; Klahr03;
Obergaulinger09). Along the same lines, Suzuki2019PASJ developed a cylindrical
shearing box model by extending the normal Cartesian box model; they applied
the conserved quantities of mass, momentum, and magnetic field in cylindrical
coordinates to periodic shearing conditions at both radial boundaries. They
demonstrated that quasi-time-steady inward mass accretion is induced by the
angular momentum flux that is outwardly transported by MRI-induced turbulence.

However, S19 performed the only single case with a thin vertical thickness of
one scale height. In addition, the treatment for the radial shearing
boundaries was not satisfactory; in particular, the difference in the
amplitude and frequency of epicyclic perturbations between both ends was
causing mismatched fluctuations that travel across the radial boundaries. In
this paper, we modify the treatment for the radial boundary condition and
conduct cylindrical MHD simulations with different ratios between the scale
height, $H_{0}$, and the radial location, $R_{0}$, in a larger vertical
domain. We examine how the amplification of magnetic fields depends on the
curvature effect arising from different aspect ratios, $H_{0}/R_{0}$, by
comparing the results to those in the exact Cartesian setup.

The structure of the paper is as follows. In Section 2, we summarize the setup
for our simulations in cylindrical shearing boxes. The shearing radial
boundary condition is modified from the original method in S19 by separating
the shearing variables into the mean and perturbation components, which is
described in Appendix A. Section 3 presents the main results of the
simulations, focusing on the transport among the different components of
magnetic fields. In Section 4 we first state the importance of the radial
gradient of epicycle frequency in the time evolution of magnetic fields and
discuss related topics; the detailed formulae for the linear perturbation
analyses on viscous-type instability is presented in Appendix B. We summarize
the paper in Section 5.

##  2 Cylindrical shearing box simulation

###  2.1 Basic setup

We solve ideal MHD equations in cylindrical coordinates, $(R,\phi,z)$,

$\frac{d\rho}{dt}+\rho\mbox{\boldmath${\nabla}$}\cdot\mbox{\boldmath${v}$}=0,$
(1)
$\rho\frac{d\mbox{\boldmath${v}$}}{dt}=-\mbox{\boldmath${\nabla}$}\left(p+\frac{B^{2}}{8\pi}\right)+\left(\frac{\mbox{\boldmath${B}$}}{4\pi}\cdot\mbox{\boldmath${\nabla}$}\right)\mbox{\boldmath${B}$}-\rho\frac{GM_{\star}}{R^{2}}\hat{R}$
(2) $\frac{\partial\mbox{\boldmath${B}$}}{\partial
t}=\mbox{\boldmath${\nabla\times(v\times B)}$},$ (3)

and

$\mbox{\boldmath${\nabla\cdot B}$}=0,$ (4)

where $\frac{d}{dt}$ and $\frac{\partial}{\partial t}$ denote Lagrangian and
Eulerian time derivatives, respectively; $\rho$, $p$, ${v}$ and ${B}$ are
density, gas pressure, velocity, and magnetic field, $G$ is the gravitational
constant, $M_{\star}$ is the mass of a central star, and the “hat” stands for
a unit vector.

We assume locally isothermal gas with an equation of state,

$p=\rho c_{\rm s}^{2},$ (5)

where $c_{\rm s}$ is isothermal sound speed that depends only on $R$ and does
not change with time. We employ the following radial dependence of temperature
($\propto c_{\rm s}^{2}$),

$c_{\rm s}^{2}=c_{{\rm s},0}^{2}\left(\frac{R}{R_{0}}\right)^{-1},$ (6)

where the subscript, $0$, indicates that the value is evaluated at the
reference radial location, $R=R_{0}$.

In the numerical simulations, we solve, instead of ${v}$ in the rest frame,
velocity, $\mbox{\boldmath${\delta v}$}=\mbox{\boldmath${v}$}-R\Omega_{\rm
eq}\hat{\phi}$, or

$(\delta v_{R},\delta v_{\phi},\delta v_{z})=(v_{R},v_{\phi}-R\Omega_{\rm
eq},v_{z}),$ (7)

if expressed in components, evaluated in the frame corotating with the
equilibrium rotation frequency, $\Omega_{\rm eq}$. If we ignore magnetic field
and assume $v_{R}=v_{z}=0$, the radial component of equation (2) yields

$-\frac{v_{\phi,{\rm
eq}}^{2}}{R}=-\frac{GM_{\star}}{R^{2}}-\frac{1}{\rho}\frac{\partial
p}{\partial R}$ (8)

under the equilibrium condition. Defining $\Omega_{\rm eq}\equiv v_{\phi,{\rm
eq}}/R$, from the radial force balance equation (8) and assuming the radial
profile of density, $\rho\propto R^{-\mu}$, we obtain

$\Omega_{\rm eq}=\Omega_{\rm K}\sqrt{1-2f_{\rm K}},$ (9)

where $\Omega_{\rm K}=\sqrt{GM_{\star}/R^{3}}$ is the Keplerian rotation
frequency and

$f_{\rm K}=-\frac{1}{\rho}\frac{\partial p}{\partial R}\bigg{/}2R\Omega_{\rm
K}^{2}=\frac{(1+\mu)c_{\rm s}^{2}}{2R^{2}\Omega_{\rm K}^{2}}$ (10)

is a sub-Keplerian parameter (ada76, S19).

Equations (1) – (5) are solved with the second-order Godunov + CMoCCT method
(Sano99; Evans88; Clarke96). The azimuthal advection due to $\Omega_{\rm eq}$
(equation 7) is handled with the “FARGO (Fast Advection in Rotating Gaseous
Objects)”-type algorithm (Masset2000A&AS; Benitez2016ApJS), which is an
improved treatment from S19 to loosen the Courant-Friedrichs-Lewy condition
(CFL1928) for the time update. The vertical component of the gravity is
ignored so that the simulations are performed in vertically unstratified
boxes, and thus the periodic boundary condition is applied to the $z$
boundaries as well as the $\phi$ boundaries.

We use the same simulation units as in S19: the physical variables are
normalized by the units of $R_{0}=1$, $\rho_{0}=1$, and $\Omega_{\rm K,0}=1$;
the magnetic field is normalized by $R_{0}\Omega_{\rm
K,0}\sqrt{4\pi\rho_{0}}$, which cancels the $\sqrt{4\pi}$ factor in the cgs-
Gauss units. In the following sections, we conventionally define
$2\pi/\Omega_{\rm K,0}$ as “one rotation time”, which is slightly shorter than
the actual time, $2\pi/\Omega_{\rm eq,0}$, it takes for one rotation at
$R=R_{0}$ in the sub-Keplerian condition. There are several conserved
quantities to verify the numerical treatment in the cylindrical shearing box;
see S19 for the details on the other numerical implementation.

###  2.2 Shearing radial boundary condition

The radial boundary condition we are adopting is an extension of the shearing
periodic condition in a local Cartesian box (Hawley1995). The basic framework
of the shearing radial boundary condition in a local cylindrical box is
explained in S19. In short, conserved quantities, instead of primitive
variables, are utilized for the boundary treatment by explicitly including the
effect of geometrical curvature. An improvement from the original prescription
in S19 is that we separate the shearing variables, $S$, composed from the
conserved quantities, into mean and perturbation components,

$S\mbox{\footnotesize($R_{\pm},\phi,z$)}=\langle
S\mbox{\footnotesize($R_{\mp},\phi\pm\Delta\Omega_{\rm eq}t,z$)}\rangle+\delta
S,$ (11)

where the subscripts, $-$ and $+$, stand for the inner and outer radial
boundaries, $\Delta\Omega_{\rm eq}=\Omega_{\rm eq,-}-\Omega_{\rm eq,+}(>0)$,
$\langle\cdots\rangle$ indicates the average over the $\phi-z$ plane, and

$\delta S=S-\langle S\rangle$ (12)

is the perturbation component (see also Section 2.5).

In the original treatment, $S$ was directly passed to the ghost cell outside
the simulation box from the corresponding cells in the simulation domain. In
this prescription, however, as the properties of the disturbances in general
differ at the inner and outer boundaries, the mismatched perturbations are
exchanged across both boundaries without any correction in an unphysical
fashion, which can be most clearly seen in the initial phase of MRI (e.g.,
Fig. 1 of S19). In order to suppress such contamination across the radial
boundaries, in $\delta S$ we include not only the perturbations from the
corresponding sheared cells but also the perturbations on the radially
adjacent cell inside the simulation domain. Additionally, $\delta S$ is
rescaled in order that the root-mean-squared (rms) amplitude, $\sqrt{\delta
S^{2}}$, in the ghost cells averaged over the $\phi-z$ surface is matched to
the surface averaged $\sqrt{\delta S^{2}}$ in the adjacent cells inside the
simulation domain. The detailed numerical implementation for the boundary
treatment is described in Appendix A.

###  2.3 Initial condition

We start the simulations with weak net vertical magnetic field. We adopt the
same radial dependences of the initial density and vertical field as those
employed in S19:

$\rho_{\rm init}=\rho_{\rm init,0}\left(\frac{R}{R_{0}}\right)^{-1}$ (13)

and

$B_{z,{\rm init}}=B_{z,0,{\rm init}}\left(\frac{R}{R_{0}}\right)^{-1}.$ (14)

Equations (6), (13), and (14) guarantee a constant initial plasma $\beta$, and
we set

$\beta_{z,{\rm init}}=\frac{8\pi\rho_{\rm init}c_{\rm s}^{2}}{B_{z,{\rm
init}}}=10^{4}.$ (15)

In order to trigger MRI, we add random velocity perturbations in $v_{R}$ and
$v_{\phi}$ to the equilibrium velocity distribution,
$(v_{R},v_{\phi},v_{z})=(0,R\Omega_{\rm eq},0)$.

###  2.4 Simulation domain & resolution

The simulation domain is a local cylindrical region that covers $(R_{-}\sim
R_{+},\phi_{-}\sim\phi_{+},z_{-}\sim z_{+})$. The azimuthal and vertical
spacings, $\Delta\phi$ and $\Delta z$, of each grid cell are constant. The
radial grid size, $\Delta R$, is proportional to $R$.

We consider two cases with thin-disk conditions of $c_{\rm
s,0}=0.1R_{0}\Omega_{\rm K,0}$ and $0.01R_{0}\Omega_{\rm K,0}$ at $R=R_{0}$.
We define the scale height, $H_{0}=c_{\rm s,0}/\Omega_{\rm K,0}$, at
$R=R_{0}$. The same box size per $H_{0}$ is adopted in these two cases,
$(L_{R},L_{\phi},L_{z})=(R_{+}-R_{-},R_{0}(\phi_{+}-\phi_{-}),z_{+}-z_{-})=(4H_{0},\frac{5\pi}{3}H_{0},4H_{0})$,
and $(L_{R},L_{\phi},L_{z})$ is resolved with grid points of
$(N_{R},N_{\phi},N_{z})=(256,320,256)$ (see Table 1 for the summary). We note
that the vertical box size is four times as large as that adopted in S19 to
cover the sufficient number of MRI channel-mode wavelengths in the saturated
state (see, e.g., Bodo2008; Johansen2009; Shi2016, for discussion on the
numerical domain size in Cartesian simulations).

Model | Simulation domain | Resolution
---|---|---
Cylindrical box | $R_{-}$ | $R_{+}$ | $L_{R}$ | $\phi_{-}$ | $\phi_{+}$ | $L_{\phi}$ | $z_{-}$ | $z_{+}$ | $L_{z}$ | $N_{R}$ | $N_{\phi}$ | $N_{z}$
$H_{0}/R_{0}=0.1$ | $0.82R_{0}$ | $1.22R_{0}$ | $4H_{0}$ | $-\pi/12$ | $\pi/12$ | $(5\pi/3)H_{0}$ | $-0.2R_{0}$ | $0.2R_{0}$ | $4H_{0}$ | $256$ | $320$ | $256$
$H_{0}/R_{0}=0.01$ | $0.980R_{0}$ | $1.020R_{0}$ | $4H_{0}$ | $-\pi/120$ | $\pi/120$ | $(5\pi/3)H_{0}$ | $-0.02R_{0}$ | $0.02R_{0}$ | $4H_{0}$ | $256$ | $320$ | $256$
Cartesian box | $x_{-}$ | $x_{+}$ | $L_{x}$ | $y_{-}$ | $y_{+}$ | $L_{y}$ | $z_{-}$ | $z_{+}$ | $L_{z}$ | $N_{x}$ | $N_{y}$ | $N_{z}$
$R_{0}\Rightarrow\infty$ | $-2H_{0}$ | $2H_{0}$ | $4H_{0}$ | $-(5\pi/6)H_{0}$ | $(5\pi/6)H_{0}$ | $(5\pi/3)H_{0}$ | $-2H_{0}$ | $2H_{0}$ | $4H_{0}$ | $256$ | $320$ | $256$
Table 1: Summary of the simulation domains and resolutions.

The numerical resolution $=64/H_{0}$ for the $R$ and $z$ directions is the
same as that in S19. It is $\approx 61/H_{0}$ in the $\phi$ direction, which
is slightly higher than $\approx 49/H_{0}$ in S19. For comparison, we also
carry out a numerical simulation in a local Cartesian box with the “same” box
size, $(L_{x},L_{y},L_{z})=(L_{R},L_{\phi},L_{z})$, per $H_{0}$ and the same
numerical resolution, $(N_{x},N_{y},N_{z})=(N_{R},N_{\phi},N_{z})$ (Table 1).

###  2.5 Spatial and temporal average

To analyze numerical data, we take various averages of simulated quantities.
We express $\langle A\rangle$ as the $\phi$\- and $z$\- integrated average of
a variable, A, at $R$:

$\langle
A\rangle\equiv\frac{\int_{z_{-}}^{z_{+}}\int_{\phi_{-}}^{\phi_{+}}d\phi
dzA}{\int_{z_{-}}^{z_{+}}\int_{\phi_{-}}^{\phi_{+}}d\phi dz}.$ (16)

We define

$[A]_{R_{1}}^{R_{2}}=\frac{\int_{R_{1}}^{R_{2}}RdR\langle
A\rangle}{\int_{R_{1}}^{R_{2}}RdR}$ (17)

as the volume-integrated average. For the average over the entire simulation
domain, we simply express

$[A]\equiv[A]_{R_{-}}^{R_{+}}.$ (18)

We write

$\overline{A}=\frac{\int_{t_{1}}^{t_{2}}dtA}{\int_{t_{1}}^{t_{2}}dt}$ (19)

as the average over time between $t_{1}$ and $t_{2}$. In all the simulated
cases, the magnetic field grows to the saturated state after $t\gtrsim 20$
rotations. Thus, we take the temporal average from $t_{1}=25$ rotations to the
end of the simulation at $t_{2}=t_{\rm final}=300$ rotations.

###  2.6 Evolution of magnetic energy

| Label | Unabbreviated expression
---|---|---
| $[\phi\Rightarrow_{{}_{\phi}}~{}\,R]$ | $\left[\frac{B_{R}}{4\pi R}\frac{\partial}{\partial\phi}(v_{R}B_{\phi})\right]/\left[\frac{B_{R}^{2}}{8\pi}\right]$
$\frac{\partial\ln\left[{B_{R}^{2}}\right]}{\partial t}$ | $[R\Rightarrow_{{}_{\phi}}~{}\,R]$ | $-\left[\frac{B_{R}}{4\pi R}\frac{\partial}{\partial\phi}(v_{\phi}B_{R})\right]/\left[\frac{B_{R}^{2}}{8\pi}\right]$
| $[R\Rightarrow_{{}_{z}}~{}\,R]$ | $-\left[\frac{B_{R}}{4\pi}\frac{\partial}{\partial z}(v_{z}B_{R})\right]/\left[\frac{B_{R}^{2}}{8\pi}\right]$
| $[z\Rightarrow_{{}_{z}}~{}\,R]$ | $\left[\frac{B_{R}}{4\pi}\frac{\partial}{\partial z}(v_{R}B_{z})\right]/\left[\frac{B_{R}^{2}}{8\pi}\right]$
| $[z\Rightarrow_{{}_{z}}~{}\,\phi]$ | $\left[\frac{B_{\phi}}{4\pi}\frac{\partial}{\partial z}(v_{\phi}B_{z})\right]/\left[\frac{B_{\phi}^{2}}{8\pi}\right]$
$\frac{\partial\ln\left[B_{\phi}^{2}\right]}{\partial t}$ | $[\phi\Rightarrow_{{}_{z}}~{}\,\phi]$ | $-\left[\frac{B_{\phi}}{4\pi}\frac{\partial}{\partial z}(v_{z}B_{\phi})\right]/\left[\frac{B_{\phi}^{2}}{8\pi}\right]$
| $[\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]$ | $-\left[\frac{B_{\phi}}{4\pi}\frac{\partial}{\partial R}(v_{R}B_{\phi})\right]/\left[\frac{B_{\phi}^{2}}{8\pi}\right]$
| $[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]$ | $\left[\frac{B_{\phi}}{4\pi}\frac{\partial}{\partial R}(v_{\phi}B_{R})\right]/\left[\frac{B_{\phi}^{2}}{8\pi}\right]$
| $[R\Rightarrow_{{}_{{}_{R}}}~{}\,z]$ | $\left[\frac{B_{z}}{4\pi R}\frac{\partial}{\partial R}(Rv_{z}B_{R})\right]/\left[\frac{B_{z}^{2}}{8\pi}\right]$
$\frac{\partial\ln\left[B_{z}^{2}\right]}{\partial t}$ | $[z\Rightarrow_{{}_{{}_{R}}}~{}\,z]$ | $-\left[\frac{B_{z}}{4\pi R}\frac{\partial}{\partial R}(Rv_{R}B_{z})\right]/\left[\frac{B_{z}^{2}}{8\pi}\right]$
| $[z\Rightarrow_{{}_{\phi}}~{}\,z]$ | $-\left[\frac{B_{z}}{4\pi R}\frac{\partial}{\partial\phi}(v_{\phi}B_{z})\right]/\left[\frac{B_{z}^{2}}{8\pi}\right]$
| $[\phi\Rightarrow_{{}_{\phi}}~{}\,z]$ | $\left[\frac{B_{z}}{4\pi R}\frac{\partial}{\partial\phi}(v_{z}B_{\phi})\right]/\left[\frac{B_{z}^{2}}{8\pi}\right]$
Table 2: Summary of the volume-integrated averages for the evolution of
magnetic energy. See equations (21) – (23) for the detail.

Taking the inner product of the induction equation (3) with
$\mbox{\boldmath${B}$}/4\pi$, we have the equation for the evolution of
magnetic energy:

$\frac{\partial}{\partial
t}\left(\frac{\mbox{\boldmath${B}$}^{2}}{8\pi}\right)=\frac{\mbox{\boldmath${B}$}}{4\pi}\mbox{\boldmath${\cdot\nabla\times(v\times
B)}$}.$ (20)

For numerical investigation, we rewrite equation (20) for $B_{R}^{2}$,
$B_{\phi}^{2}$, and $B_{z}^{2}$ separately in a logarithmic derivative form:

$\displaystyle\frac{\partial\ln B_{R}^{2}}{\partial t}$
$\displaystyle=\left(\frac{B_{R}}{4\pi
R}\frac{\partial}{\partial\phi}(v_{R}B_{\phi}-v_{\phi}B_{R})\right.$
$\displaystyle-\left.\frac{B_{R}}{4\pi}\frac{\partial}{\partial
z}(v_{z}B_{R}-v_{R}B_{z})\right)\left(\frac{B_{R}^{2}}{8\pi}\right)^{-1}$
$\displaystyle\equiv(\phi\Rightarrow_{{}_{\phi}}~{}\,R)+(R\Rightarrow_{{}_{\phi}}~{}\,R)+(R\Rightarrow_{{}_{z}}~{}\,R)+(z\Rightarrow_{{}_{z}}~{}\,R)$
(21) $\displaystyle\frac{\partial\ln B_{\phi}^{2}}{\partial t}$
$\displaystyle=\left(\frac{B_{\phi}}{4\pi}\frac{\partial}{\partial
z}(v_{\phi}B_{z}-v_{z}B_{\phi})\right.$
$\displaystyle-\left.\frac{B_{\phi}}{4\pi}\frac{\partial}{\partial
R}(v_{R}B_{\phi}-v_{\phi}B_{R})\right)\left(\frac{B_{\phi}^{2}}{8\pi}\right)^{-1}$
$\displaystyle\equiv(z\Rightarrow_{{}_{z}}~{}\,\phi)+(\phi\Rightarrow_{{}_{z}}~{}\,\phi)+(\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi)+(R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi)$
(22) $\displaystyle\frac{\partial\ln B_{z}^{2}}{\partial t}$
$\displaystyle=\left(\frac{B_{z}}{4\pi R}\frac{\partial}{\partial
R}(R(v_{z}B_{R}-v_{R}B_{z}))\right.$ $\displaystyle-\left.\frac{B_{z}}{4\pi
R}\frac{\partial}{\partial\phi}(v_{\phi}B_{z}-v_{z}B_{\phi})\right)\left(\frac{B_{z}^{2}}{8\pi}\right)^{-1}$
$\displaystyle\equiv(R\Rightarrow_{{}_{{}_{R}}}~{}\,z)+(z\Rightarrow_{{}_{{}_{R}}}~{}\,z)+(z\Rightarrow_{{}_{\phi}}~{}\,z)+(\phi\Rightarrow_{{}_{\phi}}~{}\,z).$
(23)

To examine in detail the evolution of magnetic field, we analyze the
$(i\Rightarrow_{{}_{k}}~{}\,j)$ terms on the right-hand side of equations (21)
– (23) integrated over the simulation domain (equations 17 and 18). We
summarize the volume integrated averages, $[i\Rightarrow_{{}_{k}}~{}\,j]$, in
Table 2.

Similar numerical investigation on the induction equation has been done to
understand dynamo-like magnetic evolution (Brandenburg1995ApJ; Davis2010ApJ).
An essential extension from these works is that we separate the electromotive
force, ${v\times B}$, into the shearing and compressible parts (Section 3.2).

##  3 Results

###  3.1 Time Evolution

Figure 1: 3D views of magnetic field lines (black lines) and density contour
(color) of the cylindrical simulations with $H_{0}/R_{0}=0.1$ (left) and
$0.01$ (middle) and the Cartesian simulation (right) at $t=200$ rotations.
Movies222Movies are available at
https://ea.c.u-tokyo.ac.jp/astro/Members/stakeru/research/movie/index.html are
available for the cylindrical case with $H_{0}/R_{0}=0.1$ (0-10 rotations and
70-110 rotations) and the Cartesian case (0-30 rotations).

We perform numerical simulations of the three cases in Table 1 until $t=300$
rotation time. Figure 2 presents three-dimensional (3D) snapshots of these
cases at $t=200$ rotation time (Movies of the cylindrical case with
$H_{0}/R_{0}=0.1$ and the Cartesian case are available.). The three panels
illustrate that the magnetic field is turbulent and dominated by the toroidal
($\phi$) component with a moderate level of density fluctuations. Table 3
shows various quantities regarding the time and volume averaged magnetic field
strength. While the three cases give similar values, by a close look at the
table, one can recognize that the two cylindrical cases give slightly larger
ratios of the poloidal ($R$ and $z$) to toroidal magnetic fields (6th and 7th
columns).

Model | $\overline{[\alpha_{\rm M}]}$ | $\overline{[B_{R}^{2}]}/8\pi\overline{[p]}$ | $\overline{[B_{\phi}^{2}]}/8\pi\overline{[p]}$ | $\overline{[B_{z}^{2}]}/8\pi\overline{[p]}$ | $\overline{[B_{R}^{2}]}/\overline{[B_{\phi}^{2}]}$ | $\overline{[B_{z}^{2}]}/\overline{[B_{\phi}^{2}]}$
---|---|---|---|---|---|---
$H_{0}/R_{0}=0.1$ | $6.22\times 10^{-2}$ | $2.78\times 10^{-2}$ | $1.11\times 10^{-1}$ | $1.25\times 10^{-2}$ | $0.250$ | $0.112$
$H_{0}/R_{0}=0.01$ | $5.60\times 10^{-2}$ | $2.51\times 10^{-2}$ | $1.00\times 10^{-1}$ | $1.14\times 10^{-2}$ | $0.250$ | $0.113$
Cartesian Box | $7.16\times 10^{-2}$ | $2.64\times 10^{-2}$ | $1.29\times 10^{-1}$ | $1.32\times 10^{-2}$ | $0.204$ | $0.102$
Table 3: Various dimensionless magnetic quantities averaged over $t=25-300$
rotations and the entire simulation domain. The 2nd column shows Maxwell
stress; the 3rd, 4th, and 5th columns present the radial, azimuthal and
vertical components of magnetic energy normalized by the gas pressure. The 6th
and 7th columns show the ratio of the poloidal components of magnetic energy
to the azimuthal component of magnetic energy.  Figure 2: Comparison of the
time evolutions of $[\alpha_{\rm M}]$ averaged over the simulation box for the
cylindrical cases with $H_{0}/R_{0}=0.1$ (red solid) and $0.01$ (blue dotted)
and the Cartesian case (black dashed).

Figure 2 presents the time evolution of the domain-averaged dimensionless
Maxwell stress,

$[\alpha_{\rm M}]=-\frac{[B_{R}B_{\phi}]}{4\pi[p]}.$ (24)

The cylindrical case with $H_{0}/R_{0}=0.1$ (red solid) exhibits very large
temporal variability; $[\alpha_{\rm M}]$ reaches nearly $0.3$ in the active
period between $80\lesssim t\lesssim 105$ rotations, while $[\alpha_{\rm
M}]<0.05$ is considerably small in the inactive phases. In contrast, the
Cartesian case (black dashed) gives a more time-steady evolution with
$0.05\lesssim[\alpha_{\rm M}]\lesssim 0.1$ for most of the simulation time.
The case with $H_{0}/R_{0}=0.01$ (blue dotted) shows intermediate behavior
between these two cases. In spite of the different evolutionary properties,
the time-averaged $\overline{[\alpha_{\rm M}]}$’s of these three cases are
similar (Table 3) because larger $[\alpha_{\rm M}]$ in the active phases and
smaller $[\alpha_{\rm M}]$ in the inactive phases are obtained in higher-
variability cases.

###  3.2 Evolution of magnetic field

Figure 3: Comparison of the volume-integrated averages of the variation rates
in magnetic energy (equations 21 – 23 and Table 2) of the cylindrical case
with $H_{0}/R_{0}=0.1$ (left) and the Cartesian case (right). What are shown
from top to bottom are $\frac{\partial\ln[B_{R}]^{2}}{\partial t}$ (or
$\frac{\partial\ln[B_{x}]^{2}}{\partial t}$),
$\frac{\partial\ln[B_{\phi}]^{2}}{\partial t}$ (or
$\frac{\partial\ln[B_{y}]^{2}}{\partial t}$), and
$\frac{\partial\ln[B_{z}]^{2}}{\partial t}$ in units of (rotation time)-1.
Smoothing with spline interpolation is applied to the plotted lines to display
the long-time evolution.

Figure 4: Relation among the three components of magnetic field. The top-left
triangle illustrates the physical quantities describing each arrow. The top-
right, bottom-left, and bottom-right triangles show the results of the
cylindrical cases with $H_{0}/R_{0}=0.1$ and $0.01$, and the Cartesian case,
respectively. The numerical values indicate the rate of the change of magnetic
energy, $\left\langle\frac{d\ln[B_{i}^{2}]}{dt}\right\rangle$, per rotation
time (See Table 2). The positive (negative) rates with $\geq$($\leq$) $1$
($-1$) are shown by red solid (blue dashed) arrows. We note that the sum of
the values on the arrows entering each $B_{i}$ is $=0$ when the time-steady
condition is achieved.

Figure 5: Left: The shearing term, $[i\Rightarrow_{{}_{i}}~{}\,j]$, stands for
the generation (or decay) of $B_{j}$ from $B_{i}$ by the $i$ derivative of
$v_{j}$. Right: The compressible term, $[i\Rightarrow_{{}_{j}}~{}\,i]$,
represents the convergent amplification (or divergent disperse) of $B_{i}$ by
the $j$ derivative of $v_{j}$.

We have shown that the cylindrical simulations exhibit considerably different
temporal magnetic activity from the Cartesian simulation. In order to
understand the difference, we examine the volume averaged variation rates of
the magnetic energy (equations 21 – 23 and Table 2 in Section 2.6) of the
cylindrical case with $H_{0}/R_{0}=0.1$ and the Cartesian case in Figure 3.
Although some terms, e.g., $[z\Rightarrow_{{}_{z}}~{}\,\phi]$ and
$[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]$ of the cylindrical case and
$[x\Rightarrow_{{}_{x}}~{}\,z]$ and $[y\Rightarrow_{{}_{y}}~{}\,z]$ in the
Cartesian case, exhibit relatively large time variability, most of the terms
show rather time-steady behavior; at least the sign of each term is almost
unchanged during the simulation time.

We illustrate all these $[i\Rightarrow_{{}_{k}}~{}\,j]$ terms averaged between
$t=25$ and $300$ rotations (equation (19)) in “triangle diagrams” (Figure 4),
where the result of the cylindrical case with $H_{0}/R_{0}=0.01$ is also
presented here. We note that the sum of numerical values from all four arrows
entering each $B_{i}$ is $\approx 0$, which indicates that the quasi-steady-
state condition is well satisfied if we take the temporal averages over
sufficiently long time.

The variation rates of the two cylindrical cases are similar to each other,
but they are very different from those in the Cartesian case (Figure 4). This
indicates that the physical properties of the magnetic evolution in the local
Cartesian box is fundamentally different even from those in the nearly
Cartesian box ($H_{0}/R_{0}=0.01$). We speculate that this is because the
Cartesian shearing box approximation cannot consider the radial variation of
epicycle frequency (see Section 4.1). Additionally, the magnitudes of the
variation rates are systematically larger in the cylindrical simulations. This
is expected to cause the larger temporal variability observed in the
cylindrical simulations (Figure 2).

Figure 6: An example for the amplification of a magnetic field by shearing
terms. The magnetic field is distributed with $-B_{R}B_{\phi}>0$, which is the
usual situation in the system with inner-fast differential rotation.
$B_{\phi}$ increases by radial differential rotation of $v_{\phi}$, which is
represented by $(R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi)$ (blue). $B_{R}$
increases by the inward (outward) motion of the inner (outer) region as a
result of the outward transport of angular momentum, which is represented by
$(\phi\Rightarrow_{{}_{\phi}}~{}\,R)$ (red).

We examine each $\overline{[i\Rightarrow_{{}_{k}}~{}\,j]}$ term in more
detail, particularly focusing on the similarities and differences between the
cylindrical and Cartesian simulations. Let us begin with the
$\overline{[i\Rightarrow_{{}_{i}}~{}\,j]}$ terms located on the sides of the
triangle, which stand for the transfer of magnetic fields between different
components of ${B}$; the physical meaning is the growth or decay of magnetic
fields by shearing motion (left panel of Figure 5). Although the numerical
value of each shearing term is considerably different for the cylindrical and
Cartesian simulations, the signs are the same. Thus, the physical properties
on the shearing amplification (and attenuation for the negative values) of the
magnetic fields are similar at least in a qualitative sense. Positive
$\overline{[z\Rightarrow_{{}_{z}}~{}\,R]}$ and
$\overline{[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]}$ reflect the standard
pathway triggered by MRI (e.g., Brandenburg1995ApJ; Balbus1998); $B_{R}$ is
generated from $B_{z}$ by the MRI, and eventually, $B_{\phi}$ is amplified
from the fluctuating $B_{R}$ by the winding due to the radial differential
rotation (Figure 6).

$\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,R]}$ is also positive because the
transfer of angular momentum from inner to outer rings causes the inner
(outer) ring to move inward (outward), which generates $B_{R}$ from $B_{\phi}$
(Figure 6). We would like to note that
$\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,R]}$ ($=+12.2$ and $+12.9$) in the
cylindrical simulations are much larger than
$\overline{[y\Rightarrow_{{}_{x}}~{}\,x]}(=+0.468)$ in the Cartesian
simulation. The difference is probably because the angular momentum cannot be
defined in the Cartesian setup owing to the assumed symmetry with respect to
the $\pm x$ directions (See Section 1). The leakage from the dominate
component of $B_{\phi}$ ($B_{y}$) to $B_{z}$ also gives positive
$\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,z]}$
($\overline{[y\Rightarrow_{{}_{\phi}}~{}\,z]}$) with the similar tendency
concerning the difference between the cylindrical and Cartesian cases.

Next, we inspect the compressible terms,
$\overline{[i\Rightarrow_{{}_{j}}~{}\,i]}$, which correspond to the round
arrows around each $B_{i}$ in Figure 4. The physical meaning is the
amplification (attenuation) of magnetic field by convergent (divergent) flows
(right panel of Figure 5). One can see that the signs of some
$\overline{[i\Rightarrow_{{}_{j}}~{}\,i]}$ terms are different in the
cylindrical and Cartesian simulations. Most remarkable difference is found in
the compressible amplification of $B_{\phi}$ along $R$;
$\overline{[\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]}$ in the cylindrical cases
is quite large ($=+6.42$ for $H_{0}/R_{0}=0.1$ and $+6.98$ for
$H_{0}/R_{0}=0.01$), which is in contrast to negative
$\overline{[y\Rightarrow_{{}_{z}}~{}\,y]}(=-1.05)$ in the Cartesian case.
Surprisingly, the amplification rate by
$\overline{[\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]}$ is slightly larger than
that of $\overline{[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]}$, the winding term
by radial differential rotation. In other words, the contribution from the
compressible flows dominates that from the shear flows in the amplification of
$B_{\phi}$ in the realistic cylindrical simulations, which is fundamentally
different from the result obtained from the Cartesian simulation. We again
infer that this is because of the radial variation of epicycle frequency,
which cannot be considered in the Cartesian setup.

$\overline{[R\Rightarrow_{{}_{\phi}}~{}\,R]}$ is also different between the
cylindrical and Cartesian simulations. Although
$\overline{[x\Rightarrow_{{}_{y}}~{}\,x]}$ is nearly 0 in the Cartesian
simulation, $\overline{[R\Rightarrow_{{}_{\phi}}~{}\,R]}$ takes a relatively
large negative value in the cylindrical simulations, which partially cancels
positive $\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,R]}$ due to the outward
transport of angular momentum, described above (Figure 6).

The signs of the compressible terms regarding $B_{z}$ are all opposite between
the cylindrical and Cartesian simulations, whereas the variation rates
themselves are not so large. The radial gradient of epicycle frequency is a
key for the positive $\overline{[z\Rightarrow_{{}_{{}_{R}}}~{}\,z]}$ in the
cylindrical simulations, which will be discussed in Section 4.1.

###  3.3 Onset of enhanced magnetic activity

Figure 7: Time evolution of Maxwell stress (top) and variation rates of
$[B_{\phi}^{2}]$ (middle and bottom) averaged over $R=0.95R_{0}\sim 1.05R_{0}$
in the rising and saturated phases of the most active phase in the cylindrical
simulation with $H_{0}/R_{0}=0.1$. The middle panel compares the radial
compression (red solid) and shear (blue dotted) terms. The bottom panel
presents the total compression (red solid) and shear (black dashed) terms.

We investigate in detail the mechanism of the large temporal magnetic
variability (Figure 2) observed in the cylindrical case with $H_{0}/R_{0}=0.1$
by utilizing $\overline{[i\Rightarrow_{{}_{k}}~{}\,j]}$ terms. In this
subsection, we scrutinize the rising phase of magnetic activity

Figure 7 shows the time evolution of $[\alpha_{\rm
M}]_{0.95R_{0}}^{1.05R_{0}}$ (top panel) and variation rates of the most
dominate component of the toroidal field,
$\partial\ln([B_{\phi}^{2}]_{0.95R_{0}}^{1.05R_{0}})/\partial t$, (middle and
bottom panels) during $t=65-97$ rotation time, which covers the rising and
saturated phases of the largest peak in $[\alpha_{\rm M}]$ (Figure 2). In this
figure, the volume integral is not taken in the entire domain
($=0.82R_{0}-1.22R_{0}$) but a narrower region between $R=0.95R_{0}$ and
$1.05R_{0}$ to focus on local magnetic activity. As shown in the top panel,
the maximum value of $[\alpha_{\rm M}]_{0.95R_{0}}^{1.05R_{0}}=0.340$ is
moderately larger than the maximum, $[\alpha_{\rm M}]=0.275$, averaged over
the whole domain, because the former is more sensitively affected by local
variations.

The middle panel of Figure 7 compares the variation rates of $B_{\phi}^{2}$
due to the radial compression,
$[\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]_{0.95R_{0}}^{1.05R_{0}}$ (red solid)
and the winding by radial shear,
$[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]_{0.95R_{0}}^{1.05R_{0}}$ (blue dotted).
While the radial shear term greatly fluctuates from negative to positive
values with time, the compressible term takes positive values in a more stable
manner. In particular, the compressible amplification plays an indispensable
role in the growth of $B_{\phi}^{2}$ in the rising phase of the magnetic
activity when the radial shear term is negative during $80\lesssim t\lesssim
84$ rotations.

The bottom panel of Figure 7 shows the sum of the compressible terms,
$\left\\{[\phi\Rightarrow_{{}_{z}}~{}\,\phi]+[\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]\right\\}_{0.95R_{0}}^{1.05R_{0}}$,
(red solid) and the shearing terms,
$\left\\{[z\Rightarrow_{{}_{z}}~{}\,\phi]+[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]\right\\}_{0.95R_{0}}^{1.05R_{0}}$,
(black dashed). As shown in Figures 3 & 4, the vertical compression term,
$[\phi\Rightarrow_{{}_{z}}~{}\,\phi]$, is negative, and thus, the variation
rate by the total compression is reduced from that by the only radial term
(red line in the middle panel). However, the effect of the total compression
(red line in the bottom panel) still keeps positive in the initial rising
phase of $t\lesssim 83$ rotation time; the compressible amplification plays an
essential role in triggering the bursty magnetic enhancement. This is
substantially different from the situation of the Cartesian simulation, in
which both $[y\Rightarrow_{{}_{x}}~{}\,y]$ and $[y\Rightarrow_{{}_{z}}~{}\,y]$
are always negative, namely the compressible terms work as the decay of
$B_{y}^{2}$ by divergent diffusion.

After $t\gtrsim 83$ rotations, the total compressible effect is negative as
the increased magnetic pressure countervails the compressible amplification.
After that, the shearing terms take on the role of the magnetic amplification
to the first peak at $t\approx 85$ rotations.

###  3.4 Termination of magnetic activity

Figure 8: Time evolution of magnetic energy (top) and the variation rate of
$B_{z}^{2}$ by shearing from $\phi$ to $z$ (bottom) averaged over
$R=0.95R_{0}\sim 1.05R_{0}$ in the declining phase of the same active phase
shown in Figure 7. What are shown in the top panel are the three components of
$[B_{i}^{2}]$ with $i=R$ (black dotted), $\phi$ (blue solid), and $z$ (red
dashed), normalized by the respective time-averaged values,
$\overline{[B_{i}^{2}]}$, during $t=25-300$ rotations.
$[B_{i}^{2}]/\overline{[}B_{i}^{2}]=1$ is plotted by the black thin solid
line. In the bottom panel, the time average,
$\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,z]}_{0.95R_{0}}^{1.05R_{0}}=+9.77$,
is also plotted by the black dotted line.

Next, we focus on the declining phase of the same active phase discussed in
Section 3.3. The top panel of Figure 8 presents the three components of
magnetic energy averaged over the same narrow region with
$R=0.95R_{0}-1.05R_{0}$, adopted for Figure 7. In order to see relative
enhancements, each $[B_{i}^{2}]$ is normalized by the time-averaged
$\overline{[B_{i}^{2}]}$. One can recognize that
$[B_{z}^{2}]/\overline{[B_{z}^{2}]}$ (red dashed line) reaches the highest
peaks at $t\approx 90.8$ and $105.4$ rotations among the three components just
before the drops of the magnetic energy. The times for these $z$-component
peaks occur slightly later than those for the $R$ and $\phi$ components.
Moreover, in the subsequent declining period,
$[B_{z}^{2}]/\overline{[B_{z}^{2}]}$ is larger than
$[B_{R}^{2}]/\overline{[B_{R}^{2}]}$ (black dotted line) and
$[B_{\phi}^{2}]/\overline{[B_{\phi}^{2}]}$ (blue solid line). In light of
these properties on the time evolution, we speculate that the magnetically
active states are turned off when the magnetic energy being exchanged between
the $R$ and $\phi$ components rapidly leaks to the $z$ component.

In order to verify this speculation, we inspect the transfer rate,
$[\phi\Rightarrow_{{}_{\phi}}~{}\,z]_{0.95R_{0}}^{1.05R_{0}}$, from
$B_{\phi}^{2}$ to $B_{z}^{2}$ in the bottom panel of Figure 8. It
monotonically increases in the rising phase of $t\leq 90$ rotations as
$[B_{\phi}^{2}]$ increases, and it stays at a high level, which is about
double of the time-averaged rate, slightly after the time with the maximum
$[B_{\phi}^{2}]$. Even in the descending phase, the transfer rate is
maintained at a level comparable to the time-averaged rate, which further
reduces $[B_{\phi}^{2}]$. A qualitatively similar tendency is seen before and
after the later peak at $t\approx 105$ rotations, although it is not so
obvious as for the highest peak because the field strength at the later peak
is smaller. We can conclude that the enhanced leakage of the magnetic energy
on the $R-\phi$ plane to the $z$ component terminates the high magnetic
activity.

##  4 Discussion

###  4.1 Radial dependence of epicycle frequency

Figure 9: Radial distance (horizontal) vs. time (vertical) diagram for radial
displacements by epicycle oscillations for $H_{0}/R_{0}=0.1$.

A significant difference of the cylindrical approach from the Cartesian one is
the presence of the radial variation in epicycle frequency. For the
equilibrium rotation profile, equation (9), the epicycle frequency is written
as

$\kappa=\sqrt{R\frac{\partial\Omega^{2}}{\partial R}+4\Omega^{2}}=\Omega_{\rm
eq}=\Omega_{\rm K}\sqrt{1-2f_{\rm K}},$ (25)

where we used $f_{\rm K}=H^{2}/R^{2}=H_{0}^{2}/R_{0}^{2}=$ const. (equation
10), which practically satisfied in our simulations. While in the Cartesian
box obviously $\kappa$ is spatially constant, in the cylindrical box $\kappa$
is radially dependent, $\kappa\propto R^{-3/2}$,

Figure 9 demonstrates radial displacements, $\xi_{R}=\xi_{0}\sin(\kappa t)$,
arising from epicycle oscillations in the case with $H_{0}/R_{0}=0.1$, giving
$\Omega_{\rm eq}=0.99\Omega_{\rm K}$, where the arbitrary amplitude is set to
be $\xi_{0}=0.01R_{0}$. The figure clearly illustrates that the displacements
in the neighboring radial locations gradually become out of phase with time
because $\frac{\partial\kappa}{\partial R}\neq 0$. Consequently, the phase
mixing induces converging and diverging regions by the radial component of the
oscillations, which inevitably contributes to the compressible amplification
and diffusion of magnetic field.

Figure 10: Conceptual cartoon for the role of the compressible term in the
amplification of $B_{\phi}$ (or $B_{z}$) (red lines). The initial uniform
field (left picture) is amplified in converging flows (shaded regions of right
picture) and diluted in diverging flows (white regions of right picture) by
the radial component of epicycle oscillations. The magnetic energy
$\propto[B_{\phi}^{2}]$ in the right picture is double of the initial value,
while the magnetic field strength $=[B_{\phi}]$ is conserved (see text).

Let us suppose a simple situation in which toroidal (or vertical) magnetic
field with constant strength, $B_{\phi,{\rm init}}$ (or $B_{z,{\rm init}}$),
is initially distributed (left of Figure 10). If we pick out the radial
compression term from equation (22), the corresponding term of the induction
equation is

$\frac{\partial B_{\phi}}{\partial t}=\cdots-\frac{\partial}{\partial
R}(v_{R}B_{\phi})+\cdots.$ (26)

This is essentially the continuity equation for $B_{\phi}$, whereas the
geometrical curvature term, $\frac{1}{R}\frac{\partial}{\partial R}(R\cdots)$,
is omitted if it is compared to the mass continuity equation (1). Therefore,
$B_{\phi}$ is compressed (rarefied) in the converging (diverging) regions
(right of Figure 10). If we compare the left and right pictures of Figure 10,
the volume integrated field strength is conserved because the number of the
field lines does not change. On the other hand, the magnetic energy becomes
double as the regions with strong field, $2B_{\phi,{\rm init}}$, occupy 50% of
the volume (the shaded regions in the right picture of Figure 10) so that
$[B_{\phi}^{2}]=[2B_{\phi,{\rm init}}^{2}]\times 0.5=2[B_{\phi,{\rm
init}}^{2}]$.

The consideration based on this simple model indicates that the radial
compression term should systematically increase $B_{\phi}^{2}$, which does not
occur in the Cartesian setup with the spatially constant $\kappa$. Although
the evolution of the magnetic field is more complex in realistic situations as
discussed in Section 3, this simple mechanism is expected to work in the
nonlinear saturated states, and therefore, the radial compression,
$\overline{[\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]}$, takes the relatively
large positive values in the cylindrical simulations, which is in clear
contrast to the negative $\overline{[y\Rightarrow_{{}_{x}}~{}\,y]}$ in the
Cartesian simulation (Figure 4). The same argument can be applied to the
vertical magnetic field. Indeed, Figure 4 shows the radial compression of
$B_{z}$, $\overline{[z\Rightarrow_{{}_{{}_{R}}}~{}\,z]}$, is positive in the
cylindrical cases, which is also in contrast to the negative corresponding
term, $\overline{[z\Rightarrow_{{}_{x}}~{}\,z]}$, in the Cartesian case.

The abovementioned argument indicates that the presence of
$\frac{\partial\kappa}{\partial R}\neq 0$ significantly changes the
fundamental properties of the amplification of magnetic fields. This is the
reason why the nearly Cartesian simulation with $H_{0}/R_{0}=0.01$ gives the
very different conversion rates of $[i\Rightarrow_{{}_{k}}~{}\,j]$ from the
exact Cartesian case but the similar values to the moderately cylindrical case
with $H_{0}/R_{0}=0.1$ as shown in the triangle diagrams of Figure 4.

The next question one might ask would be regarding the timescale for the onset
of the compressible amplification. When the oscillatory phases are deviated
between the radially “neighboring” regions each other by $\pi$, the
compressible amplification works most effectively as shown in Figure 10. We
can estimate the time, $t=\tau_{\rm comp}$, to reach the phase difference of
$\Delta\phi_{\kappa}=\pi$ from the uniform oscillation at $t=0$ via
$\Delta\kappa_{0}t=\Delta\phi_{\kappa}$, where the deviation of $\kappa_{0}$
is derived from $\Delta\kappa=|\left(\frac{\partial\kappa}{\partial
R}\right)_{0}|\Delta R_{0}$ by using the radial spacing between the
neighboring regions, $\Delta R_{0}$. Since
$|\left(\frac{\partial\kappa}{\partial R}\right)_{0}|=\frac{3\Omega_{\rm
eq,0}}{2R}$, we have

$\displaystyle\tau_{\rm comp}$
$\displaystyle=\frac{20\pi}{3}\left(\frac{H_{0}/R_{0}}{0.1}\right)^{-1}\left(\frac{\Delta\phi_{\kappa}}{\pi}\right)\left(\frac{\Delta
R}{H_{0}}\right)^{-1}$
$\displaystyle=\frac{10}{3}\left(\frac{H_{0}/R_{0}}{0.1}\right)^{-1}\left(\frac{\Delta\phi_{\kappa}}{\pi}\right)\left(\frac{\Delta
R}{H_{0}}\right)^{-1}[{\rm rotations}]$ (27)

where the radial spacing is normalized by $H_{0}$ for a typical scale of the
system. This equation shows that the epicycle oscillation becomes completely
out of phase ($\Delta\phi_{\kappa}=\pi$) at radial spacing of $\Delta R=H_{0}$
when $\tau_{\rm comp}\approx 3$ rotations for the cylindrical case with
$H_{0}/R_{0}=0.1$. For $H_{0}/R_{0}=0.01$, longer $\tau_{\rm comp}\approx 30$
rotations are required. However, we would like to note that the estimate based
on equation (27) is rather conservative because the effect of the radial
compression is regarded to be already effective well before the phase
difference reaches $\pi$. Thus, even in this nearly Cartesian case, the
compressible amplification is already working at $t\approx 10$ rotation time,
which is comparable to the timescale for the transition from the initial
linear stage to the nonlinear saturated phase (Figure 2).

Another characteristic feature expected from equation (27) is that, as time
goes on, the compressible amplification is going to work for smaller $\Delta
R$; smaller-scale localized magnetic concentrations can be formed at later
times. In order to capture these fine-scale structures, numerical simulations
require sufficiently high resolution. Thus, the saturation level and time-
variability of magnetic field may depend on numerical resolution, which will
be addressed in our future work.

Global simulations should automatically take the effect of
$\frac{\partial\kappa}{\partial R}\neq 0$ into account. Therefore, if one
applied the same analyses on the magnetic evolution in Section 3.2 to global
simulation data, the similar result regarding the importance of the
compressible effect should be obtained, although such an attempt has not been
tried within our knowledge. However, here, numerical resolution would matter
as discussed above.

###  4.2 Substructures and viscous-type instability

Figure 11: Face-on view of dimensionless density ($\rho/\rho_{0}$; left),
magnetization with respect of $B_{z}$ ($B_{z}^{2}/4\pi p$; middle), and
dimensionless Maxwell stress ($\alpha_{\rm M}$; right) on the midplane ($z=0$)
at $t=82.4$ rotations. Velocity fields are also shown by arrows. We note that
the velocity vectors are directed along the perturbation component, $(\delta
v_{R},\delta v_{\phi})=(v_{R},v_{\phi}-R\Omega_{\rm eq})$ (equation 7). Movie
during $t=70-110$ rotations is available.  Figure 12: Physical picture of the
viscous-type instability. See text for the detailed explanation.  Figure 13:
Scatter plots of dimensionless Maxwell stress, $\alpha_{\rm M}$, against
$B_{z}^{2}/4\pi\rho c_{\rm s}^{2}$ (left), $(\rho/\rho_{0})^{-1}$ (middle),
and $B_{z}^{2}/4\pi$ (right) at $t=82.4$ rotations. The data are taken from
the 2D plots in Figure 11. The slope of $y\propto x$ is overlaid in the left
and middle panels by the red solid line, and $y\propto x^{0.5}$ is in the left
panel by the blue dashed line.

Figure 10 implies that converging and diverging regions move with time. In the
nonlinear stage, oscillatory motions are excited randomly by turbulence and
those at different radial positions undergo phase shift each other. Hence, it
is expected that the locations of converging regions change with time in a
more or less stochastic manner. These randomly excited compressed regions
could be seeds of substructures by connecting to various types of instability
(e.g., Suriano2019MNRAS; Cui2021, see Lesur2022 for recent review).

Figure 11 shows a two-dimensional (2D) $r-\phi$ slice of density,
$\rho/\rho_{0}$ (left), “magnetization”, $B_{z}^{2}/4\pi p$ (middle), defined
as the inverse of $\beta_{z}$ (c.f., equation 15), and Maxwell stress,
$\alpha_{\rm M}$ (right) on $z=0$ at $t=82.4$ rotation time, which is in the
onset phase of the high magnetic activity analyzed in Section 3.3 (see Figures
2 & 7). One can recognize stripe structures with anti-correlated density and
magnetic field; low-density and strong-magnetic field regions are sandwiched
between denser and weaker-field regions with typical spacing of a fraction of
$H_{0}$. These rings and gaps are ubiquitously and continuously formed as
shown in the movie file. Such ring-and-gap structures are already seen in
early global unstratified ideal MHD simulations (Steinacker2002ApJ) and
discussed from a viewpoint of viscous-type instability (Lightman1974)
generalized by including the physical properties of the MRI (Hawley2001).

One of the key ingredients of the viscous-type instability in the MHD
framework is magnetic diffusion. Although the resistivity is in principle zero
under the ideal MHD approximation, numerical simulations inevitably undergo
numerical diffusion. From a physical point of view, reconnections of turbulent
magnetic fields would also provide effective magnetic diffusion and
dissipation even when the ideal MHD condition is satisfied (Lazarian1999ApJ;
Lazarian2020Phpl). We carry out linear perturbation analyses for viscous-type
instability with explicitly taking into account resistivity in the induction
equation (3). We basically follow the formulation introduced by Riols2019, but
restrict our analyses to the case without the mass loss and angular momentum
removal by MHD disk winds. An important modification from the original setting
in Riols2019 is that dimensionless viscosity, $\alpha_{\nu}$, is assumed to
depend separately on density and magnetic field:

$\alpha_{\nu}\propto\left(\frac{\rho_{0}}{\rho}\right)^{q_{\rho}}\left(\frac{B_{z}^{2}}{B_{z,0}^{2}}\right)^{q_{B}},$
(28)

where we can practically assume $\alpha_{\nu}\approx\alpha_{\rm M}$ to
interpret our simulation results. We apply plane-wave expansions to the
linearized equations of mass continuity (equation 1), radial and angular
momenta (equation 2), and vertical magnetic field (equation 3) with
resistivity. We finally obtain the criterion for the presence of an unstable
mode is

$q_{\rho}>1,$ (29)

where the detailed derivation is presented in Appendix B.

We note that in Riols2019 $\alpha_{\nu}$ is assumed to depend on
magnetization, $B_{z}^{2}/4\pi\rho c_{\rm s}^{2}$, namely $q=q_{B}=q_{\rho}$
is imposed, and they derived the instability condition as $q>1$. The relation
between time- and volume-averaged $\alpha_{\rm M}$ and $B_{z}^{2}/4\pi\rho
c_{\rm s}^{2}$ have been extensively examined with Cartesian shearing box
simulations (Salvesen2016; Scepi2018) and global simulations (Mishra2020).
There is a rough consensus of $q\approx 0.5-0.7$, which does not satisfy the
instability condition in Riols2019. On the other hand, our analyses in
Appendix B shows that $q_{B}$ does not qualitatively affect the stability
criterion but only quantitatively controls the growth (or decay) rate.

Figure 12 illustrates the physical picture of the instability. Let us suppose
a situation in which a higher-density, $\delta\rho>0$, region (shaded area) is
created by random perturbations. If $\delta(\rho\alpha_{\nu})<0$ there, the
outward transport of the angular momentum is suppressed in this denser ring
(red arrows). As a result, at the inner (outer) edge of the ring the angular
momentum increases (decreases) to give $\delta v_{\phi}>(<)$ $0$ (black
arrows), which causes the inner (outer) edge to move outward (inward), $\delta
v_{R}>(<)$ $0$ (blue arrows). Hence, the density of the denser ring further
increases. If we pick out the dependence on density333In the ideal MHD
condition, the dependence on $B_{z}$ should also be considered here as $\rho$
and $B_{z}$ behave “in phase” each other. However, the inclusion of
resistivity breaks this constraint and it is justified to focus only on the
dependence on $\rho$; see Appendix B for the detailed algebra. in equation
(28) and take the Taylor expansion,
$\delta(\rho\alpha_{\nu})\approx(1-q_{\rho})\alpha_{\nu,0}\delta\rho$.
Therefore, the condition for the instability, $\delta(\rho\alpha_{\nu})<0$,
corresponds to $q_{\rho}>1$ (equation 29).

Figure 13 shows scatter plots between $\alpha_{\rm M}$ (vertical axis) and
various quantities (horizontal axis) of each grid point displayed in Figure
11. In the left panel, the correlation with magnetization is shown. The
dependence is roughly consistent with the slope derived in the abovementioned
previous works, whereas the plots are largely scattered, reflecting the
scatter in $B_{z}^{2}$ (right panel), because the data are not averaged over
time or domain. The middle panel, which shows the dependence on the inverse of
density, exhibits relatively tighter correlation particularly in the
larger-$\alpha_{\rm M}$ and lower-density (upper right) side. Moreover, the
slope is slightly steeper than $q_{\rho}=1$ and is in the unstable regime
(equation 29). We examined scatter plots for $\rho$-$\alpha_{M}$ in other time
frames and found that $q_{\rho}\gtrsim 1$ is kept in most of the time.
Therefore, we can interpret that the substructures seen in our simulations
(Figure 11) are amplified with a secular timescale (see Appendix B for the
derivation) and maintained in the cylindrical disk that are in the marginally
unstable condition regarding the generalized viscous-type instability.

###  4.3 Intermittency

Cartesian shearing box simulations for MRI often show quasi periodicities with
$t\sim 10$ rotation time (Davis2010ApJ; Gressel2010; Guan2011ApJ), which is
related to the recurrent growth and breakup of large-scale channel-mode flows
(Gogichaishvili2018ApJ). While quasi-periodic magnetic activity is more
clearly seen in vertically stratified simulations, being frequently associated
with mass outflows (e.g., Suzuki2009ApJ; Suzuki2010ApJ; Wissing2022), similar
periodicity is also observed in unstratified simulations (Sano2001; Shi2016).
Our Cartesian case is also showing mild time variation in $[\alpha_{\rm M}]$
(black dashed line in Figure 2).

On the other hand, the duration $\sim 50$ rotations between high-magnetic-
activity phases seen in the cylindrical case with $H_{0}/R_{0}=0.1$ (red solid
line in Figure 2) is considerably longer than the typical period $\sim$ 10
rotations described above. We speculate that this weak periodicity is at least
partially controlled by the compressible effect arising from the radial
variation of $\kappa$ (Section 4.1). Utilizing equation (27), we estimate the
time to reach the phase difference, $\Delta\phi_{\kappa}=2\pi$, between
neighboring radial regions as a possible source of the periodicity. One of the
plausible choices for the radial spacing, $\Delta R$, could be the width of a
ring (or a gap). If we adopt $\Delta R\approx 0.2R_{0}$, referring to Figure
11, we obtain $t\approx 33$ rotation time, which is roughly consistent with
the observed cycle of the enhanced magnetic activity, although this argument
does not exclude other possibilities.

###  4.4 Accretion structure in active and inactive phases

Figure 14: Comparison of time-averaged radial profiles of various physical
quantities of the case with $H_{0}/R_{0}=0.1$ in the inactive (left) and
active (right) periods. The top panels show angular momentum fluxes carried by
Maxwell stress (black solid), Reynolds stress (red dotted), and mean accretion
flow (blue dashed). The middle panels compare rms $B_{R}$ (magenta dash-
dotted), $B_{\phi}$ (blue dashed), and $B_{z}$ (black solid). The initial
vertical field strength, $B_{z,{\rm init}}$, is multiplied by a factor of 10
(red dotted), in order to compare to the field strength in the saturated
state. The bottom panels show density (black solid) in comparison to the
initial distribution (red dotted) and the inverse of plasma $\beta=8\pi
p/B^{2}$ (blue dashed). We note that the vertical scales of the top and middle
panels are different for the left and right panels. See Section 2.1 for the
physical units on the vertical axes.

We compare the radial dependences of various physical quantities during
magnetically active and inactive phases in Figure 14. In order to see the
difference between the active and inactive periods, we take the temporal
averages for $[\alpha_{\rm M}]<0.13$ (left) and $[\alpha_{\rm M}]>0.13$
(right). The inactive (active) phase occupies 246.2 (28.8) rotations out of
the total time-averaged duration of $275(=300-25)$ rotations.

The top panels of Figure 14 compare the radial profiles of angular momentum
fluxes, following S19. When the time-steady condition is satisfied, an
equation for angular momentum integrated over the $\phi$-$z$ plane is

$\frac{\partial}{\partial R}\left[R^{2}\left(\overline{\langle\rho
v_{R}\rangle}R\Omega_{\rm eq}+\overline{\langle\rho v_{R}\delta
v_{\phi}\rangle}-\frac{\overline{\langle
B_{R}B_{\phi}\rangle}}{4\pi}\right)\right]=0,$ (30)

where the first (blue dashed), second (red dotted), and third (black solid)
terms indicate the angular momentum carried by the mean accretion flow, the
turbulent Reynolds stress, and the Maxwell stress, respectively. As shown in
the top left panel, the steady accretion structure is well realized in the
inactive phase; the inward angular momentum flux by the accretion flow is
balanced with the outward flux by the Maxwell stress with a small contribution
from the Reynolds stress. On the other hand, in the active phase (top right
panel) the steady accretion structure is broken because the increased magnetic
pressure pushes gas inward (outward) from the inner (outer) boundary.

The middle panels of Figure 14 show each component of the rms magnetic field
strength, $\sqrt{\overline{\langle B_{i}^{2}\rangle}}$, in comparison to the
initial $B_{z,{\rm init}}$. The left panel indicates that the vertical
magnetic field strength, $\sqrt{\overline{\langle B_{z}^{2}\rangle}}$,
averaged over the inactive periods increases about 10 times from the initial
value with almost keeping the initial profile, $\propto R^{-1}$. The radial
and azimuthal field strengths are larger (Table 3) but the radial dependences
are similar to that of the $z$ component. In the active phase (middle right
panel), while the field strength is larger as expected, the radial dependence
of each component and the ratios between different components are not so
different from those in the inactive phases.

In the bottom panel of Figure 14, density and the inverse of plasma $\beta$,

$\overline{\langle\beta\rangle}^{-1}=\frac{\overline{\langle
B^{2}\rangle}}{8\pi\overline{\langle\rho\rangle}},$ (31)

are presented. The density structure is moderately altered from the initial
condition, $R^{-1}$, as it is affected by temporal zonal flows; the density in
the inner region slightly decreases, while it increases in the outer region.
As a result, $\overline{\langle\beta\rangle}^{-1}$ exhibits a weak radial
dependence, which can be more easily seen in the active phase.

###  4.5 Boundary effect

The radial distributions averaged over the active phases show sharply pointed
and discontinuous profiles in $B^{2}$ (middle right panel of Figure 14) and
$\rho$ (bottom right panel) near the inner and outer boundaries because the
steady accretion structure is not achieved (top right panel). This implies
that the evolution of magnetic energy discussed in Section 3.2 may suffer the
treatment of the radial boundaries during the active phases. To inspect the
boundary effect, we compare all $[i\Rightarrow_{{}_{k}}~{}\,j]$ terms of the
case with $H_{0}/R_{0}=0.1$ ($0.01$) in a smaller radial region between
$R=0.90R_{0}$ ($0.990R_{0}$) and $1.11R_{0}$ ($1.010R_{0}$) avoiding the
boundary effect to those averaged over the entire simulation domain shown in
Figure 4.

Two winding terms, $\overline{[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]}$ and
$\overline{[z\Rightarrow_{{}_{z}}~{}\,\phi]}$, in the case with
$H_{0}/R_{0}=0.1$ are slightly affected; specifically, the averages in the
narrow region are $\overline{[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]}=4.12$ and
$\overline{[z\Rightarrow_{{}_{z}}~{}\,\phi]}=-4.41$, which are both reduced by
$\approx 20\%$ from the domain-averaged values. On the other hand, the
modifications in the other winding terms and all the compressible terms are
less than 5%. In the case with $H_{0}/R_{0}=0.01$, the modifications of all
the terms are less than 5%, because the magnetic activity is relatively weak
so that the unphysical boundary effect is almost ignorable. In summary, we can
conclude that the effect regarding the boundary treatment is limited and then
the discussion based on the numerical results so far is unaffected.

##  5 Summary

Continuing from S19, we studied fundamental MHD properties of accretion disks
by cylindrical shearing box simulations. We modified the treatment for the
radial boundary condition from the original prescription in S19. The key
improvement is to separate the boundary variables into the mean and
perturbation components and the amplitude of the latter is adjusted to match
the fluctuations at both boundaries (Section 2.2 and Appendix A). This
modified prescription enables us to reduce unmatched fluctuations traveling
across the boundaries.

The radial gradient of epicycle frequency, $\kappa$, causes the phase mixing
of random oscillatory motions. As a result, the radial compression is
significant in amplifying the azimuthal and vertical magnetic fields (Section
4.1). This is an important finding in this paper by inspecting the spatial
derivative terms on the right-hand-side of the equations for magnetic energy.
In contrast, the compressible effect works as diverging dilution of the
magnetic energy in the Cartesian box simulation because of the absence of the
radial variation of $\kappa$ (Section 3.2).

The compressible amplification plays a significant role in enhanced bursty
magnetic activity, which is more clearly seen in the case with larger
$H_{0}/R_{0}$ (Section 3.1). This is expected from the argument on the
timescale of the phase shift due to the non-uniform distribution of $\kappa$
(Section 4.1). We also speculate that the phase-shift timescale is related to
the weak periodicity in the intermittent magnetic bursts (Section 4.3).

The compressible amplification is also expected to create seeds of small-scale
substructures. Indeed, there are narrow ring-gap structures continuously and
ubiquitously formed in the simulations. These structures show the anti-
correlation between density and magnetic field strength; in particular, the
steep dependence of Maxwell stress on density, $\alpha_{\rm
M}\propto\rho^{-q_{\rho}}$, with $q_{\rho}\gtrsim 1$ is obtained (Section
4.2). We revisited the viscous-type instability by considering the dependence
of $\alpha_{\rm M}$ separately on density and vertical magnetic field, and
found that the instability condition is $q_{\rho}>1$ (Appendix B). Thus, we
interpreted that the ring-gap structures are maintained in the simulated disks
that are under marginally unstable conditions.

In this paper, as we focused on the effects of the disk curvature,
$H_{0}/R_{0}$, we did not conduct simulations with different initial vertical
field strengths, box sizes, or numerical resolutions. However, the dependences
on these parameters are obviously important, which will be addressed in our
future papers.

A key physics in our work is the radial gradient of $\kappa$; its nonlinear
outcomes are time-variability, intermittency, and localized substructures
observed in the simulations. This problem can be framed as nonlinear processes
in a system where the eigen-mode frequency varies spatially. When the physical
quantities are non-uniformly distributed, as being found in various
astrophysical systems, such as interstellar medium in star-forming regions and
the interior, atmosphere and magnetosphere of stars and compact objects,
eigen-mode frequencies of various oscillatory modes are also expected to be
spatially dependent. This study could have potential applications in such
systems, which is one of the future directions.

Numerical computations were carried out on Cray XC50 at Center for
Computational Astrophysics, National Astronomical Observatory of Japan, and
Yukawa-21 (Dell PowerEdge R840) at YITP, Kyoto University. This work is
supported by Grants-in-Aid for Scientific Research from the MEXT/JSPS of
Japan, 17H01105, 21H00033 and 22H01263 and by Program for Promoting Research
on the Supercomputer Fugaku by the RIKEN Center for Computational Science
(Toward a unified view of the universe: from large-scale structures to
planets; grant 20351188-PI J. Makino) from the MEXT of Japan.

##  Appendix A Modified shearing radial boundary condition

Figure 15: Shearing boundary condition in a local cylindrical box. The
simulation domain, which is enclosed by the solid black line, is between
$R=R_{-}$ and $R_{+}$. $R_{-,{\rm d}}$ ($R_{+,{\rm d}}$) is the radial
location of the innermost (outermost) cells in the simulation domain, which
are illustrated by gray shade. $R_{-,{\rm g}}$ ($R_{+,{\rm g}}$) is the radial
location of the inner (outer) ghost cells. The physical quantities of the
ghost cell at $R=R_{\pm,{\rm g}}$ (red) are determined from those of the
corresponding sheared cells at $R=R_{\mp,{\rm d}}$ (blue) and the radially
adjacent cell(s) at $R=R_{\pm,{\rm d}}$ (green). See text for the detail.

We describe the specific method to implement the shearing boundary condition,
equation (11). As we stated in Section 2.2, the modification from the original
prescription in S19 is that the shearing variables are separated into the mean
values and perturbations (equation 12). We rewrite equation (11) in a more
precise way to distinguish ghost (“g”) and domain (“d”) cells (Figure 15):

$\displaystyle S\mbox{\footnotesize($R_{\pm,{\rm g}},\phi$)}$
$\displaystyle=\langle S\mbox{\footnotesize($R_{\mp,{\rm d}},\phi-(\Omega_{\rm
eq,\pm,g}-\Omega_{\rm eq,\mp,d})t$)}\rangle+\delta S$
$\displaystyle\approx\langle S\mbox{\footnotesize($R_{\mp,{\rm
d}},\phi+\Delta\Omega_{\rm eq}t$)}\rangle+\delta S,$ (A1)

where we omitted “$z$” in the arguments as it is redundant and can be
understood without it. We employ the standard periodic boundary condition for
the mean components, $\langle S\rangle$: the shearing variables averaged over
the innermost (outermost) $\phi$-$z$ plane inside the simulation domain (gray
shaded regions in Figure 15) are passed to the outer (inner) ghost cells
outside the domain (regions surrounded by the dashed lines in Figure 15).

On the other hand, for the perturbation components, $\delta S$, in a ghost
cell at $R_{\pm,{\rm g}}$(red in Figure 15), we consider the contribution from
not only the corresponding sheared cells at $R_{\mp,{\rm d}}$ (blue) but also
the adjacent cells at $R_{\pm,{\rm d}}$ (green). Specifically, we adopt the
following parameterization:

$\displaystyle\delta S\mbox{\footnotesize($R_{\pm,{\rm g}},\phi$)}$
$\displaystyle=f_{\rm amp,\pm}\left[f_{\rm sh}K\delta
S\mbox{\footnotesize($R_{\mp,{\rm d}},\phi-(\Omega_{\rm eq,\pm,g}-\Omega_{\rm
eq,\mp,d})t$)}+(1-f_{\rm sh})\delta S\mbox{\footnotesize($R_{\pm,{\rm
d}},\phi-(\Omega_{\rm eq,\pm,g}-\Omega_{\rm eq,\pm,d})t$)}\right]$ (A2)
$\displaystyle\approx f_{\rm amp,\pm}\left[f_{\rm sh}K\delta
S\mbox{\footnotesize($R_{\mp,{\rm d}},\phi+\Delta\Omega_{\rm eq}t$)}+(1-f_{\rm
sh})\delta S\mbox{\footnotesize($R_{\pm,{\rm d}},\phi$)}\right],$ (A3)

where $f_{\rm amp,\pm}(\approx 1)$ is a parameter that controls the amplitude
of $\delta S$, $f_{\rm sh}$ determines the fractional contributions from the
sheared position at the opposite side of the $R$ domain (blue in Figure 15),
and

$K=\sqrt{\frac{\langle\delta S^{2}\mbox{\footnotesize($R_{\pm,{\rm
d}},\phi-(\Omega_{\rm eq,\pm,g}-\Omega_{\rm
eq,\pm,d})t$)}\rangle}{\langle\delta S^{2}\mbox{\footnotesize($R_{\mp,{\rm
d}},\phi-(\Omega_{\rm eq,\pm,g}-\Omega_{\rm
eq,\mp,d})t$)}\rangle}}\approx\sqrt{\frac{\langle\delta
S^{2}\mbox{\footnotesize($R_{\pm,{\rm d}},\phi$)}\rangle}{\langle\delta
S^{2}\mbox{\footnotesize($R_{\mp,{\rm d}},\phi\pm\Delta\Omega_{\rm
eq}t$)}\rangle}}.$ (A4)

Since in general the relative azimuthal position between the ghost and domain
cells, which changes with time, is not an exact integer multiple of the
azimuthal size of the fixed grid cell, $\delta S$’s in the first and second
terms on the right-hand side of equation (A2) need to be interpolated from two
grid cells. However, as we ignore the tiny difference between $\Omega_{\rm
eq,\pm,d}$ and $\Omega_{\rm eq,\pm,g}$ and use equation (A3) in the numerical
implementation, $\delta S$ in the second term is taken from the adjacent grid
cell in the $R$ direction.

We note that $f_{\rm sh}=1$ recovers the original treatment of S19, and
$f_{\rm sh}=0$ corresponds to the “non-gradient” boundary condition for
$\delta S$. Throughout the current paper, we adopt $f_{\rm sh}=0.5$, namely
the fluctuations at both the shearing location and the neighboring position
are equally mixed. We also note that, when $f_{\rm amp,\pm}=1$, the correction
by this $K$ ensures that the rms amplitude averaged over the ghost cells on
$R=R_{\pm,{\rm g}}$ matches that averaged over the neighboring cells at
$R=R_{\pm,{\rm d}}$: $\sqrt{\delta S^{2}\mbox{\footnotesize($R_{\pm,{\rm
g}},\phi$)}}=\sqrt{\delta S^{2}\mbox{\footnotesize($R_{\pm,{\rm d}},\phi$)}}$.
MRI grows from earlier times at smaller $R$. Without this amplitude
correction, it can be easily seen that fluctuations in the inner regions leak
out of $R=R_{-}$ and seeps from $R=R_{+}$ to the quiet outer region in the
initial growth phase, as exhibited in Figure 1 of S19.

Figure 16: 3D snapshot of the case with $\beta_{z.{\rm init}}=10^{3}$ and the
vertical box size, $L_{z}=H$, at $t=3.20$ rotation. Movie between 0 and 10
rotations is available.

Figure 16 shows the initial growing stage of MRI with the amplitude correction
by equation (A4). The simulation setup is the same as in S19: $\beta_{z,{\rm
init}}=10^{3}$ and the vertical domain size, $L_{z}=H$, is smaller than that
adopt in the present paper. The figure is demonstrating that the contamination
of the inner fluctuations passed to the outer boundary, $R=R_{+}$, is greatly
suppressed, compared to that found in S19444As the FARGO advection scheme
adopted in this paper suffers from less numerical diffusion than the normal
advection scheme used in S19, MRI sets in from slightly earlier time. Thus, we
are showing the snapshot at $t=3.20$ rotations, which is earlier than $t=3.45$
and $4.00$ rotations presented in S19..

$f_{\rm amp,\pm}$ is tuned to realize the steady accretion structure as shown
in the top left panel of Figure 14. Larger $f_{\rm amp,\pm}$ yields larger
perturbations in the ghost zone at $R=R_{\pm,{\rm g}}$, which results in
larger effective turbulent pressure there. Therefore, larger $f_{\rm amp,+}$
($f_{\rm amp,-}$) tends to induce inward (outward) flows. We are adopting
$f_{\rm amp,+}=1.02$ and $f_{\rm amp,-}=0.93$ for the cylindrical case with
$H_{0}/R_{0}=0.1$ and $f_{\rm amp,+}=0.94$ and $f_{\rm amp,-}=0.93$ for
$H_{0}/R_{0}=0.01$. A different choice of $f_{\rm amp,\pm}$ disturbs the
balanced accretion profile. However, we note that the time-averaged magnetic
profiles are not significantly affected, provided that $f_{\rm amp,\pm}\approx
1$ is employed.

##  Appendix B Linear perturbation analyses for viscous-type instability

We conduct linear perturbation analyses on the 2D $r-\phi$ plane under the
axisymmetric approximation. In the momentum equation (2), we do not explicitly
include magnetic fields but consider them through $\alpha$ prescription
(Shakura73). Then, the radial and angular momentum equations can be written as

$\frac{\partial}{\partial t}(\rho w_{R})=2\rho\Omega_{\rm
K}w_{\phi}-\frac{\partial}{\partial R}(\rho c_{\rm
s}^{2})+\rho\frac{w_{\phi}^{2}}{R}$ (B1)

and

$\frac{\partial}{\partial t}(\rho
w_{\phi}R)+\frac{1}{R}\frac{\partial}{\partial R}(R^{2}\rho c_{\rm
s}^{2}\alpha_{\nu})+\frac{1}{2}\rho\Omega_{\rm K}w_{R}R=0$ (B2)

in Eulerian forms, where $(w_{R},w_{\phi})=(v_{R},v_{\phi}-R\Omega_{\rm K})$
is the velocity deviation from Keplerian rotation and $\alpha_{\nu}$ is
dimensionless viscosity, equation (28). We note that, in equation (B1), the
inward gravity, $-\rho\frac{GM}{R^{2}}$, and the centrifugal force, $\rho
R\Omega_{\rm K}^{2}$, should be present but they are exactly canceled out each
other. We adopt the radial profiles of $c_{\rm s}^{2}\propto R^{-1}$ (equation
6), $\rho\propto R^{-1}$, (equation 13) and $\alpha_{\nu}\propto R^{1/2}$
(S19) to match our simulation setting. Then, from the unperturbed state
($\partial_{t}=0$) of equations (B1) and (B2), we obtain

$w_{\phi}=R\Omega_{\rm K}(\sqrt{1-f_{\rm K}}-1)\propto R^{-\frac{1}{2}},$ (B3)

which is consistent with equation (9), and

$w_{R}(=v_{R})=-\frac{\alpha_{\nu}c_{\rm s}^{2}}{R\Omega_{\rm
K}}=-\frac{\alpha_{\nu}H^{2}\Omega_{\rm K}}{R}\propto R^{0}.$ (B4)

$\alpha_{\nu}$ is assumed to depend separately on $\rho$ and $B_{z}$ (equation
28), and thus, we have

$\alpha_{\nu}\approx\alpha_{\nu,0}\left(\frac{R}{R_{0}}\right)^{\frac{1}{2}}\left(1-q_{\rho}\frac{\delta\rho}{\rho_{0}}+2q_{B}\frac{\delta
B_{z}}{B_{z,0}}\right)\equiv\alpha_{\nu,0}^{\prime}\left(1-q_{\rho}\frac{\delta\rho}{\rho_{0}}+2q_{B}\frac{\delta
B_{z}}{B_{z,0}}\right).$ (B5)

Here, density is calculated from the continuity equation (1) and vertical
magnetic field from the induction equation (3) with resistivity, $\eta$,
explicitly included:

$\frac{\partial B_{z}}{\partial t}=\frac{1}{R}\frac{\partial}{\partial
R}\left[R\left(-w_{R}B_{z}+\eta\frac{\partial B_{z}}{\partial
R}\right)\right],$ (B6)

where we assumed $v_{z}=0$. We also adopt the $\alpha$ prescription for $\eta$
with the same dependence on $\rho$ and $B_{z}$,

$\eta=\alpha_{\eta}\frac{c_{\rm s}^{2}}{\Omega_{\rm
K}}=\alpha_{\eta}H^{2}\Omega_{\rm
K}=\eta_{0}\left(\frac{R}{R_{0}}\right)\left(1-q_{\rho}\frac{\delta\rho}{\rho_{0}}+2q_{B}\frac{\delta
B_{z}}{B_{z,0}}\right).$ (B7)

We note that $\eta\propto R$ ensures that equation (B6) has an unperturbed-
state solution of $B_{z}\propto R^{-1}$, which is consistent with equation
(14) (see also S19).

We expand $\rho$, $w_{R}$, $w_{\phi}$, and $B_{z}$ with the radial dependences
that match our simulation setting:

$\displaystyle\rho$
$\displaystyle=(\rho_{0}+\delta\rho)\left(\frac{R}{R_{0}}\right)^{-1}$
$\displaystyle\equiv\rho_{0}^{\prime}+\delta\rho^{\prime},$ (B8)
$\displaystyle w_{R}$ $\displaystyle=w_{R,0}+\delta w_{R}$
$\displaystyle\equiv w_{R,0}^{\prime}+\delta w_{R}^{\prime},$ (B9)
$\displaystyle w_{\phi}$ $\displaystyle=(w_{\phi,0}+\delta
w_{\phi})\left(\frac{R}{R_{0}}\right)^{-\frac{1}{2}}$ $\displaystyle\equiv
w_{\phi,0}^{\prime}+\delta w_{\phi}^{\prime},$ (B10) $\displaystyle B_{z}$
$\displaystyle=(B_{z,0}+\delta B_{z})\left(\frac{R}{R_{0}}\right)^{-1}$
$\displaystyle\equiv B_{z,0}^{\prime}+\delta B_{z}^{\prime}.$ (B11)

We apply these expansions to equations (1), (B1), (B2), and (B6) and further
apply plane-wave expansion,

$\delta\propto\exp(i\omega t-ikR),$ (B12)

to the first-order variables. Then, the corresponding four equations are

$i(\omega-kw_{R,0}^{\prime})\delta\rho^{\prime}-ik\rho_{0}^{\prime}\delta
w_{R}^{\prime}=0,$ (B13)

$i\omega\rho_{0}^{\prime}\delta
w_{R}^{\prime}-2\rho_{0}^{\prime}\left[\Omega_{\rm
K}+\frac{w_{\phi,0}^{\prime}}{R}\right]\delta w_{\phi}^{\prime}+\left[i\omega
w_{R.0}^{\prime}-2\Omega_{\rm K}w_{\phi,0}^{\prime}-c_{\rm
s}^{2}\left(\frac{2}{R}+ik\right)-\frac{{w_{\phi,0}^{\prime}}^{2}}{R}\right]\delta\rho^{\prime}=0,$
(B14)

$i\omega\rho_{0}^{\prime}\delta w_{\phi}^{\prime}+\frac{1}{2}\Omega_{\rm
K}\rho_{0}^{\prime}\delta w_{R}^{\prime}+\left[i\omega
w_{\phi,0}^{\prime}+\frac{1}{2}\Omega_{\rm K}w_{R,0}^{\prime}+c_{\rm
s}^{2}\alpha_{\nu}^{\prime}(1-q_{\rho})\left(\frac{1}{2R}-ik\right)\right]\delta\rho^{\prime}+c_{\rm
s}^{2}\alpha_{\nu}^{\prime}\frac{\rho_{0}^{\prime}}{B_{z,0}^{\prime}}\left(\frac{1}{R}-2ik\right)q_{B}\delta
B_{z}^{\prime}=0,$ (B15)

and

$\left(i\omega-
ikw_{R,0}^{\prime}+\eta_{0}k^{2}+ik\eta_{0}(1-2q_{B})\frac{1}{R}\right)\delta
B_{z}^{\prime}-ikB_{z,0}^{\prime}\delta
w_{R}^{\prime}+ik\eta_{0}q_{\rho}\frac{1}{R}\frac{B_{z,0}^{\prime}}{\rho_{0}^{\prime}}\delta\rho^{\prime}=0.$
(B16)

From equations (B13) - (B16), we can derive a fourth-order equation with
respect to $\omega$. Here, we ignore the solutions of
$\omega^{2}\approx\Omega_{\rm K}^{2}+k^{2}c_{\rm s}^{2}$ describing sound
waves in shearing systems, and focus on secular modes by assuming
$\omega^{2}\ll\Omega_{\rm K}^{2}$. In addition, we also assume $H/R\ll 1$ and
ignore the terms $\propto 1/R,1/R^{2},\cdots$ arising from the curvature of
cylindrical coordinates. Then, we finally obtain

$\displaystyle\omega^{2}(1+k^{2}H^{2})$
$\displaystyle-\omega\left[(kw_{R}+ik^{2}\eta)(1+k^{2}H^{2})+2i(1-q_{\rho}+2q_{B})k^{2}H^{2}\Omega_{\rm
K}\alpha_{\nu}\right]$ $\displaystyle+2\Omega_{\rm
K}\alpha_{\nu}k^{2}H^{2}\left[(1-q_{\rho}+2q_{B})ikw_{R}-(1-q_{\rho})k^{2}\eta\right]=0,$
(B17)

where we omitted the subscripts, 0, and the superscripts, $\prime$, because we
ignored the $R$-dependent terms. Substituting equations (B4) and (B7) into
equation (B17), we can derive

$\omega=\frac{k^{2}H^{2}\Omega_{\rm
K}\alpha_{\nu}}{2(1+k^{2}H^{2})}\left[i\left\\{2(1-q_{\rho}+2q_{B})+\frac{1+k^{2}H^{2}}{{\rm
Pm}}\right\\}\pm\sqrt{-\left\\{2(1-q_{\rho}+2q_{B})+\frac{1+k^{2}H^{2}}{{\rm
Pm}}\right\\}^{2}+8(1-q_{\rho})\frac{1+k^{2}H^{2}}{{\rm Pm}}}\right],$ (B18)

where

${\rm Pm}\equiv\frac{\alpha_{\nu}}{\alpha_{\eta}}$ (B19)

is the magnetic Prandtl number and we once again disregarded the terms
$\propto 1/R$ that originally involved $w_{R}$. As shown in Figure 13,
$q_{\rho}\approx 1$, and then, in the square-root of equation (B18) the first
term should dominate the second term. Thus, we take the Taylor expansion of
the square-root term. Recalling the form of the plane-wave expansion (equation
B12), one can find a growing mode, if it is present, only for the negative
sign of equation (B18). The growth rate can be calculated as

$s=i\omega=\frac{2(q_{\rho}-1)k^{2}H^{2}}{2(1-q_{\rho}+2q_{B}){\rm
Pm}+1+k^{2}H^{2}}\Omega_{\rm
K}\alpha_{\nu}=\begin{cases}2(q_{\rho}-1)\frac{k^{2}H^{2}}{1+k^{2}H^{2}}\Omega_{\rm
K}\alpha_{\nu}&{\rm Pm}\ll 1\\\
\frac{(q_{\rho}-1)}{(1-q_{\rho}+2q_{B})}k^{2}H^{2}\Omega_{\rm
K}\alpha_{\eta}&{\rm Pm}\gg 1\\\ \end{cases}.$ (B20)

This result shows that the growth (or decay) rate is determined by the slower
process of viscous or resistive diffusion; it is an order of $\Omega_{\rm
K}\alpha_{\nu}$ for Pm $\ll 1$ and $\Omega_{\rm K}\alpha_{\eta}$ for Pm $\gg
1$. In the former case, an unstable mode is present if $q_{\rho}>1$ (equation
29), and the growth rate does not depends on $q_{B}$. In the latter case,
although it depends on $q_{B}$, the instability condition is again
$q_{\rho}>1$ because $(1-q_{\rho}+2q_{B})>0$ is probably satisfied in usual
situations (Figure 13).